diff --git a/global_3d.c b/global_3d.c
index 44c265b..60ec2b8 100644
--- a/global_3d.c
+++ b/global_3d.c
@@ -27,6 +27,9 @@
 /* Choice of potential */
 
 #define POT_HARMONIC 1      /* harmonic oscillator */
+#define POT_COULOMB 2       /* Coulomb (1/r) potential */
+#define POT_PERIODIC 3      /* periodic potential */
+#define POT_DOUBLE_COULOMB 4      /* sum of Coulomb potentials located at focal points of ellipse */
 
 /* plot types used by rde */
 
@@ -39,10 +42,12 @@
 #define Z_ANGLE_GRADIENTX 231   /* direction of gradient of u */
 #define Z_NORM_GRADIENT_INTENSITY 24  /* gradient and intensity of polar angle */
 #define Z_VORTICITY 25  /* curl of polar angle */
+#define Z_VORTICITY_ABS 251  /* absolute value of curl of polar angle */
 
 /* for Schrodinger equation */
 #define Z_MODULE 30       /* module squared of first two fields */
 #define Z_ARGUMENT 31     /* argument of first two fields, with luminosity depending on module */
+#define Z_REALPART 32     /* first field, with luminosity depending on module */
 
 /* for RPSLZ equation */
 #define Z_MAXTYPE_RPSLZ 40      /* color of type with maximal density */
@@ -51,8 +56,8 @@
 
 /* macros to avoid unnecessary computations in 3D plots */
 
-#define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY))
-#define COMPUTE_THETAZ ((zplot == Z_POLAR)||(zplot == Z_NORM_GRADIENT)||(zplot == Z_ANGLE_GRADIENT)||(zplot == Z_NORM_GRADIENT_INTENSITY)||(zplot == Z_VORTICITY))
+#define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY)||(cplot == Z_VORTICITY_ABS))
+#define COMPUTE_THETAZ ((zplot == Z_POLAR)||(zplot == Z_NORM_GRADIENT)||(zplot == Z_ANGLE_GRADIENT)||(zplot == Z_NORM_GRADIENT_INTENSITY)||(zplot == Z_VORTICITY)||(zplot == Z_VORTICITY_ABS))
 
 #define COMPUTE_ENERGY ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)||(zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_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))
 
diff --git a/global_ljones.c b/global_ljones.c
index 42d3027..4b7747e 100644
--- a/global_ljones.c
+++ b/global_ljones.c
@@ -52,6 +52,9 @@
 #define S_POOL_TABLE 6      /* pool table with pockets */
 #define S_CENTRIFUGE_RND 7  /* segments forming centrifuge with more rounded bins */
 #define S_CENTRIFUGE_LEAKY 8  /* segments forming centrifuge with rounded bins and holes */
+#define S_CIRCLE_EXT 9      /* segments forming a repelling cicle */
+#define S_ROCKET_NOZZLE 10  /* segments forming a rocket bell-shpaed nozzle */
+#define S_TWO_CIRCLES_EXT 11  /* segments forming two repelling cicle */
 
 /* particle interaction */
 
@@ -86,6 +89,12 @@
 
 #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 */
+
+/* Gravity schedules */
+
+#define G_INCREASE_RELEASE 1    /* slow increase and instant release */
+#define G_INCREASE_DECREASE 2   /* slow increase an decrease */
 
 /* Plot types */
 
@@ -196,6 +205,7 @@ typedef struct
     double length;              /* length of segment */
     short int concave;          /* corner is concave, to add extra repelling force */
     short int cycle;            /* set to 1 if (x2, y2) is equal to (x1, y1) of next segment */
+    short int group;            /* group to which segment belongs (for several obstacles) */
     double angle1, angle2;      /* angles in which concave corners repel */
     short int active;           /* segment is active */
     double x01, x02, y01, y02;  /* initial values of extremities, in case of rotation/translation */
diff --git a/global_particles.c b/global_particles.c
index 5f67491..3e011ca 100644
--- a/global_particles.c
+++ b/global_particles.c
@@ -77,6 +77,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7;
 #define P_POLYGON 5       /* regular polygon, alternative for D_POLYGON */
 #define P_TOKA_PRIME 6    /* Tokarsky room made of 86 triangles */
 #define P_TREE 7          /* pine tree */
+#define P_TOKA_NONSELF 8  /* Tokarsky non-self-unilluminable room */
 
 /* Color palettes */
 
diff --git a/global_pdes.c b/global_pdes.c
index 07738bd..916bd06 100644
--- a/global_pdes.c
+++ b/global_pdes.c
@@ -53,9 +53,12 @@
 #define D_STAR 42               /* star shape */
 #define D_FRESNEL 43            /* Fresnel lens */
 #define D_NOISEPANEL 44         /* zigzag noise insulating panel */
+#define D_NOISEPANEL_RECT 441   /* comparison between zigzag noise insulating panel and flat walls */
 #define D_DOUBLE_FRESNEL 45     /* two facing Fresnel lenses */
 #define D_QRD 46                /* quadratic resonance diffuser */
+#define D_QRD_ASYM 461          /* asymmetric quadratic resonance diffuser */
 #define D_CIRCLE_SEGMENT 47     /* lens-shaped circular segment */
+#define D_GROOVE 48             /* groove array supposed to induce polaritons */
 
 #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) */
@@ -96,7 +99,8 @@
 #define D_VONKOCH_HEATED 27 /* von Koch snowflake in larger circle */
 
 /* Variable index of refraction */
-#define IOR_MANDELBROT 1    /* index of refraction depends on escape time in Mandelbrot set */
+#define IOR_MANDELBROT 1      /* index of refraction depends on escape time in Mandelbrot set (log) */
+#define IOR_MANDELBROT_LIN 100    /* index of refraction depends on escape time in Mandelbrot set (linear) */
 
 /* Boundary conditions */
 
@@ -154,6 +158,8 @@
 
 #define COL_TURBO_CYCLIC 101    /* TURBO color palette (by Anton Mikhailov) corrected to be cyclic, beta */
 
+#define NPWIDTH 0.02      /* width of noise panel separation */
+
 
 typedef struct
 {
diff --git a/heat.c b/heat.c
index 816d119..298544c 100644
--- a/heat.c
+++ b/heat.c
@@ -35,7 +35,7 @@
 #include <tiffio.h>     /* Sam Leffler's libtiff library. */
 #include <omp.h>
 
-#define MOVIE 0         /* set to 1 to generate movie */
+#define MOVIE 1         /* set to 1 to generate movie */
 
 /* General geometrical parameters */
 
@@ -671,7 +671,7 @@ void animation()
     glutSwapBuffers();
     
     draw_wave(phi, xy_in, 1.0, 0);
-    if (DRAW_BILLIARD) draw_billiard();
+    if (DRAW_BILLIARD) draw_billiard(0, 1.0);
 //     print_Julia_parameters(i);
     
 //     print_level(MDEPTH);
@@ -698,7 +698,7 @@ void animation()
         
         for (j=0; j<NVID; j++) evolve_wave(phi, phi_tmp, xy_in);
 
-        if (DRAW_BILLIARD) draw_billiard();
+        if (DRAW_BILLIARD) draw_billiard(0, 1.0);
         
 //         print_level(MDEPTH);
 //         print_Julia_parameters(i);
diff --git a/lennardjones.c b/lennardjones.c
index de3d945..fc27464 100644
--- a/lennardjones.c
+++ b/lennardjones.c
@@ -39,11 +39,10 @@
 #define MOVIE 0         /* set to 1 to generate movie */
 #define DOUBLE_MOVIE 0  /* set to 1 to produce movies for wave height and energy simultaneously */
 
-#define TIME_LAPSE 0     /* set to 1 to add a time-lapse movie at the end */
+#define TIME_LAPSE 1     /* set to 1 to add a time-lapse movie at the end */
                          /* so far incompatible with double movie */
 #define TIME_LAPSE_FACTOR 3    /* factor of time-lapse movie */
 
-
 /* General geometrical parameters */
 
 #define WINWIDTH 	1280  /* window width */
@@ -54,15 +53,15 @@
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
-#define INITXMIN -0.7
-#define INITXMAX 0.7	/* x interval for initial condition */
-#define INITYMIN -0.5
-#define INITYMAX 0.5	/* y interval for initial condition */
+#define INITXMIN -1.95
+#define INITXMAX 1.95	/* x interval for initial condition */
+#define INITYMIN -1.05
+#define INITYMAX 1.05	/* y interval for initial condition */
 
 #define BCXMIN -2.0
-#define BCXMAX 5.0	/* x interval for boundary condition */
-#define BCYMIN -1.6
-#define BCYMAX 1.6	/* y interval for boundary condition */
+#define BCXMAX 2.0	/* x interval for boundary condition */
+#define BCYMIN -1.125
+#define BCYMAX 1.125	/* y interval for boundary condition */
 
 #define OBSXMIN -2.0
 #define OBSXMAX 2.0     /* x interval for motion of obstacle */
@@ -70,10 +69,10 @@
 #define CIRCLE_PATTERN 8   /* pattern of circles, see list in global_ljones.c */
 
 #define ADD_FIXED_OBSTACLES 0   /* set to 1 do add fixed circular obstacles */
-#define OBSTACLE_PATTERN 2  /* pattern of obstacles, see list in global_ljones.c */
+#define OBSTACLE_PATTERN 3  /* 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 5   /* 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 9   /* pattern of repelling segments, see list in global_ljones.c */
 
 #define TWO_TYPES 0         /* set to 1 to have two types of particles */
 #define TPYE_PROPORTION 0.7 /* proportion of particles of first type */
@@ -89,14 +88,14 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 100        /* number of points for Poisson C_RAND_POISSON arrangement */
-#define PDISC_DISTANCE 3.33  /* minimal distance in Poisson disc process, controls density of particles */
+#define PDISC_DISTANCE 4.0  /* minimal distance in Poisson disc process, controls density of particles */
 #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.7	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
 #define MU 0.012  	    /* parameter controlling radius of particles */
-#define MU_B 0.018         /* parameter controlling radius of particles of second type */
-#define NPOLY 18             /* number of sides of polygon */
+#define MU_B 0.018          /* parameter controlling radius of particles of second type */
+#define NPOLY 20            /* number of sides of polygon */
 #define APOLY 0.666666666   /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 4            /* depth of computation of Menger gasket */
 #define MRATIO 3            /* ratio defining Menger gasket */
@@ -107,6 +106,7 @@
 #define NGRIDY 24           /* 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 X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -115,10 +115,11 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2000      /* number of frames of movie */
-#define NVID 500         /* number of iterations between images displayed on screen */
-#define NSEG 250         /* number of segments of boundary */
-#define INITIAL_TIME 50    /* time after which to start saving frames */
+#define NSTEPS 3000      /* number of frames of movie */
+#define NVID 80         /* number of iterations between images displayed on screen */
+#define NSEG 150         /* number of segments of boundary */
+#define INITIAL_TIME 100     /* time after which to start saving frames */
+#define OBSTACLE_INITIAL_TIME 50     /* 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 */
@@ -132,14 +133,16 @@
 
 /* Boundary conditions, see list in global_ljones.c */
 
-#define BOUNDARY_COND 3
+#define BOUNDARY_COND 0
 
 /* Plot type, see list in global_ljones.c  */
 
 #define PLOT 0
 #define PLOT_B 8        /* plot type for second movie */
 
+#define DRAW_BONDS 0    /* set to 1 to draw bonds between neighbours */
 #define COLOR_BONDS 1   /* set to 1 to color bonds according to length */
+#define FILL_TRIANGLES 0    /* set to 1 to fill triangles between neighbours */
 
 /* Color schemes */
 
@@ -162,48 +165,49 @@
 
 /* particle properties */
 
-#define ENERGY_HUE_MIN 330.0      /* color of original particle */
-#define ENERGY_HUE_MAX 50.0        /* color of saturated particle */
+#define ENERGY_HUE_MIN 330.0        /* color of original particle */
+#define ENERGY_HUE_MAX 50.0         /* color of saturated particle */
 #define PARTICLE_HUE_MIN 359.0      /* color of original particle */
 #define PARTICLE_HUE_MAX 0.0        /* color of saturated particle */
-#define PARTICLE_EMAX 2.0e2           /* energy of particle with hottest color */
-#define HUE_TYPE0 45.0     /* hue of particles of type 0 */
-#define HUE_TYPE1 300.0     /* hue of particles of type 1 */
-#define HUE_TYPE2 300.0      /* hue of particles of type 2 */
-#define HUE_TYPE3 300.0     /* hue of particles of type 3 */
+#define PARTICLE_EMAX 1.0e3         /* energy of particle with hottest color */
+#define HUE_TYPE0 280.0     /* hue of particles of type 0 */
+#define HUE_TYPE1 70.0     /* hue of particles of type 1 */
+#define HUE_TYPE2 180.0      /* hue of particles of type 2 */
+#define HUE_TYPE3 210.0     /* hue of particles of type 3 */
 
 #define RANDOM_RADIUS 0     /* set to 1 for random circle radius */
-#define DT_PARTICLE 2.0e-6    /* time step for particle displacement */
+#define DT_PARTICLE 3.0e-6    /* time step for particle displacement */
 #define KREPEL 12.0          /* constant in repelling force between particles */
-#define EQUILIBRIUM_DIST 4.5    /* Lennard-Jones equilibrium distance */
+#define EQUILIBRIUM_DIST 5.0    /* Lennard-Jones equilibrium distance */
 #define EQUILIBRIUM_DIST_B 3.5  /* 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 1.0e-1     /* damping coefficient of particles */
-#define PARTICLE_MASS 4.0   /* mass of particle of radius MU */
+#define DAMPING 25.0    /* damping coefficient of particles */
+#define PARTICLE_MASS 1.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 10.0        /* initial velocity range */
+#define V_INITIAL 0.0        /* initial velocity range */
 #define OMEGA_INITIAL 10.0        /* initial angular velocity range */
 
 #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.02           /* initial inverse temperature */
+#define BETA 0.1            /* initial inverse temperature */
 #define MU_XI 0.01           /* friction constant in thermostat */
 #define KSPRING_BOUNDARY 1.0e11    /* confining harmonic potential outside simulation region */
 #define KSPRING_OBSTACLE 1.0e11    /* harmonic potential of obstacles */
-#define NBH_DIST_FACTOR 4.0       /* radius in which to count neighbours */
-#define GRAVITY 0.0            /* gravity acting on all particles */
+#define NBH_DIST_FACTOR 6.0       /* radius in which to count neighbours */
+#define GRAVITY 1000.0            /* gravity acting on all particles */
 #define INCREASE_GRAVITY 0     /* set to 1 to increase gravity during the simulation */
-#define GRAVITY_FACTOR 100.0   /* factor by which to increase gravity */
-#define GRAVITY_INITIAL_TIME 500    /* time at start of simulation with constant gravity */
-#define GRAVITY_RESTORE_TIME 1000    /* time at end of simulation with gravity restored to initial value */
+#define GRAVITY_SCHEDULE 2     /* type of gravity schedule, see list in global_ljones.c */
+#define GRAVITY_FACTOR 50.0    /* factor by which to increase gravity */
+#define GRAVITY_INITIAL_TIME 200    /* time at start of simulation with constant gravity */
+#define GRAVITY_RESTORE_TIME 700    /* time at end of simulation with gravity restored to initial value */
 
 #define ROTATION 0           /* set to 1 to include rotation of particles */
 #define COUPLE_ANGLE_TO_THERMOSTAT 0    /* set to 1 to couple angular degrees of freedom to thermostat */
 #define DIMENSION_FACTOR 1.0  /* scaling factor taking into account number of degrees of freedom */  
-#define KTORQUE 50.0          /* force constant in angular dynamics */
+#define KTORQUE 100.0          /* force constant in angular dynamics */
 #define KTORQUE_B 10.0        /* force constant in angular dynamics */
 #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 */
@@ -214,10 +218,10 @@
 #define QUADRUPOLE_RATIO 0.6  /* anisotropy in quadrupole potential */ 
 
 #define INCREASE_BETA 1  /* set to 1 to increase BETA during simulation */
-#define BETA_FACTOR 0.025   /* factor by which to change 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 FINAL_CONSTANT_PHASE 1000  /* final phase in which temperature is constant */
+#define FINAL_CONSTANT_PHASE 0  /* 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 */
@@ -238,9 +242,11 @@
 #define N_P_AVERAGE 100      /* size of pressure averaging window */
 #define N_T_AVERAGE 50       /* size of temperature averaging window */
 #define MAX_PRESSURE 3.0e10  /* pressure shown in "hottest" color */
-#define PARTIAL_THERMO_COUPLING 0   /* set to 1 to couple only particles to the right of obstacle to thermostat */
-#define PARTIAL_THERMO_REGION 1     /* region for partial thermostat coupling (see list in global_ljones.c) */
-#define PARTIAL_THERMO_SHIFT 0.5    /* distance from obstacle at the right of which particles are coupled to thermostat */
+#define PARTIAL_THERMO_COUPLING 1   /* set to 1 to couple only particles to the right of obstacle to thermostat */
+#define PARTIAL_THERMO_REGION 2     /* region for partial thermostat coupling (see list in global_ljones.c) */
+#define PARTIAL_THERMO_SHIFT 0.2    /* distance from obstacle at the right of which particles are coupled to thermostat */
+#define PARTIAL_THERMO_WIDTH 0.5    /* vertical size of partial thermostat coupling */
+#define PARTIAL_THERMO_HEIGHT 0.2   /* vertical size of partial thermostat coupling */
 
 #define INCREASE_KREPEL 0   /* set to 1 to increase KREPEL during simulation */
 #define KREPEL_FACTOR 1000.0   /* factor by which to change KREPEL during simulation */
@@ -250,8 +256,8 @@
 #define NPART_BOTTOM 100     /* number of particles at the bottom */
 
 #define ADD_PARTICLES 0    /* set to 1 to add particles */
-#define ADD_TIME 500       /* time at which to add first particle */
-#define ADD_PERIOD 250       /* time interval between adding further particles */
+#define ADD_TIME 0       /* time at which to add first particle */
+#define ADD_PERIOD 10000       /* time interval between adding further particles */
 #define FINAL_NOADD_PERIOD 200  /* 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 */
 
@@ -263,16 +269,21 @@
 
 #define ROTATE_BOUNDARY 0           /* set to 1 to rotate the repelling segments */
 #define SMOOTH_ROTATION 1           /* set to 1 to update segments at each time step (rather than at each movie frame) */
-#define PERIOD_ROTATE_BOUNDARY 2500 /* period of rotating boundary */
+#define PERIOD_ROTATE_BOUNDARY 1000  /* period of rotating boundary */
 #define ROTATE_INITIAL_TIME 0       /* initial time without rotation */
-#define ROTATE_FINAL_TIME 500       /* final time without rotation */
+#define ROTATE_FINAL_TIME 100       /* final time without rotation */
+#define ROTATE_CHANGE_TIME 0.33     /* relative duration of acceleration/deceleration phases */
+#define OMEGAMAX 100.0              /* maximal rotation speed */
 #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 MOVE_BOUNDARY 1        /* set to 1 to move repelling segments, due to force from particles */
-#define SEGMENTS_MASS 100.0     /* mass of collection of segments */
-#define DEACTIVATE_SEGMENT 1    /* set to 1 to deactivate last segment after a certain time */
+#define MOVE_BOUNDARY 0        /* set to 1 to move repelling segments, due to force from particles */
+#define SEGMENTS_MASS 5.0      /* mass of collection of segments */
+#define DEACTIVATE_SEGMENT 0    /* set to 1 to deactivate last segment after a certain time */
 #define SEGMENT_DEACTIVATION_TIME 1000   /* time at which to deactivate last segment */
+#define SEGMENTS_X0 0.0         /* initial position of segments */
+#define SEGMENTS_Y0 -0.75         /* initial position of segments */
 
 #define POSITION_DEPENDENT_TYPE 0   /* set to 1 to make particle type depend on initial position */
 #define POSITION_Y_DEPENDENCE 0     /* set to 1 for the separation between particles to be vertical */
@@ -280,7 +291,7 @@
 
 #define REACTION_DIFFUSION 0    /* set to 1 to simulate a chemical reaction (particles may change type) */
 #define RD_TYPES 3              /* number of types in reaction-diffusion equation */
-#define REACTION_PROB 0.03      /* probability controlling reaction term */ 
+#define REACTION_PROB 0.0045      /* probability controlling reaction term */ 
 
 #define PRINT_PARTICLE_NUMBER 0     /* set to 1 to print total number of particles */
 #define PRINT_LEFT 1        /* set to 1 to print certain parameters at the top left instead of right */
@@ -297,11 +308,11 @@
 
 #define FLOOR_FORCE 1      /* set to 1 to limit force on particle to FMAX */
 #define FMAX 1.0e12         /* maximal force */
-#define FLOOR_OMEGA 1      /* set to 1 to limit particle momentum to PMAX */
+#define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
 #define PMAX 1000.0        /* maximal force */
 
-#define HASHX 100   /* size of hashgrid in x direction */
-#define HASHY 40   /* size of hashgrid in y direction */
+#define HASHX 60   /* size of hashgrid in x direction */
+#define HASHY 30    /* 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 */
 
@@ -312,6 +323,7 @@
 
 #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO))
 #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE))
+#define TWO_OBSTACLES (SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)
 
 double xshift = 0.0;      /* x shift of shown window */
 double xspeed = 0.0;      /* x speed of obstacle */
@@ -320,10 +332,10 @@ double vylid = 0.0;       /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c
 double xwall = 0.0;       /* x coordinate of wall (for BC_RECTANGLE_WALL b.c.) */
 double vxwall = 0.0;      /* x speed of wall (for BC_RECTANGLE_WALL b.c.) */
 double angular_speed = 0.0;    /* angular speed of rotating segments */
-double xsegments = 0.0;        /* x coordinate of segments (for option MOVE_BOUNDARY) */
-double ysegments = 0.0;        /* y coordinate of segments (for option MOVE_BOUNDARY) */
-double vxsegments = 0.0;       /* vx coordinate of segments (for option MOVE_BOUNDARY) */
-double vysegments = 0.0;       /* vy coordinate of segments (for option MOVE_BOUNDARY) */
+double xsegments[2] = {SEGMENTS_X0, -SEGMENTS_X0};  /* x coordinate of segments (for option MOVE_BOUNDARY) */
+double ysegments[2] = {SEGMENTS_Y0, SEGMENTS_Y0};  /* y coordinate of segments (for option MOVE_BOUNDARY) */
+double vxsegments[2] = {0.0, 0.0};       /* vx coordinate of segments (for option MOVE_BOUNDARY) */
+double vysegments[2] = {0.0, 0.0};       /* vy coordinate of segments (for option MOVE_BOUNDARY) */
 int thermostat_on = 1;    /* thermostat switch used when VARY_THERMOSTAT is on */
 
 #define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on)))
@@ -440,31 +452,68 @@ double obstacle_schedule_smooth(int i, int j)
 
 double gravity_schedule(int i, int j)
 {
-    double time, gravity;
+    double time, gravity, x, y;
     
-    if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY);
-    else 
-    {
-        time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) 
-        + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME);
-        gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0));
-//         printf("i = %i, time = %.3lg, Gravity = %.3lg\n", i, time, gravity); 
-        return(gravity);
+    switch (GRAVITY_SCHEDULE){
+    
+        case (G_INCREASE_RELEASE):
+        {
+            if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY);
+            else 
+            {
+                time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) 
+                + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME);
+                gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0));
+                return(gravity);
+            }
+            break;
+        }
+        
+        case (G_INCREASE_DECREASE): 
+        {
+            if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY);
+            else 
+            {
+                time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) 
+                + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME);
+                x = 2.0 - cos(DPI*time);
+                y = 0.5*((GRAVITY_FACTOR - 1.0)*x + 3.0 - GRAVITY_FACTOR);
+                gravity = GRAVITY*y;
+                return(gravity);
+            }
+            break;
+        }
     }
 }
 
 double rotation_angle(double phase)
 {
-    double omegamax = 15.0;
-    
     /* case of rotating hourglass */
-    while (phase > DPI) phase -= DPI;
-    return(phase - 0.5*sin(2.0*phase));
+//     while (phase > DPI) phase -= DPI;
+//     return(phase - 0.5*sin(2.0*phase));
     
     /* case of centrifuge */
-//     while (phase > DPI) phase -= DPI;
-//     angular_speed = 0.5*omegamax*(1.0 - cos(phase));
-//     return(0.5*omegamax*(phase - sin(phase)));
+//     while (phase > 1.0) phase -= 1.0;
+//     phase *= DPI;
+//     angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase));
+//     return(0.5*OMEGAMAX*(phase - sin(phase)));
+    
+    /* case of centrifuge remaining at constant speed for a while */
+    if (phase < ROTATE_CHANGE_TIME)
+    {
+//         angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase*PI/ROTATE_CHANGE_TIME));
+        return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME)));
+    }
+    else if (phase < 1.0 - ROTATE_CHANGE_TIME)
+    {
+//         angular_speed = OMEGAMAX;
+        return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME));
+    }
+    else
+    {
+//         angular_speed = 0.5*OMEGAMAX*(1.0 + cos((phase - 1.0 + ROTATE_CHANGE_TIME)*PI/ROTATE_CHANGE_TIME));
+        return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME)));
+    }
 }
 
 double rotation_schedule(int i)
@@ -487,7 +536,7 @@ double rotation_schedule(int i)
 
 double rotation_schedule_smooth(int i, int j)
 {
-    double phase;
+    double phase, angle, phase1, angle1;
     static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME;
     
     if (i < imin) 
@@ -497,9 +546,22 @@ double rotation_schedule_smooth(int i, int j)
     }
     else
     {
-        if (i > imax) i = imax;
-        phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*((double)(i - imin) + (double)j/(double)NVID);
-        return(rotation_angle(phase));
+        if (i > imax) 
+        {
+            angle = rotation_angle(1.0);
+            angular_speed = 0.0;
+        }
+        else 
+        {
+            phase = (1.0/(double)(imax - imin))*((double)(i - imin) + (double)j/(double)NVID);
+            angle = rotation_angle(phase);
+            
+            phase1 = (1.0/(double)(imax - imin))*((double)(i + 1 - imin) + (double)j/(double)NVID);
+            angle1 = rotation_angle(phase1);
+            angular_speed = 25.0*(angle1 - angle);
+        }
+        
+        return(angle);
     }
 }
 
@@ -660,37 +722,68 @@ void evolve_wall(double fboundary)
 
 void evolve_segments(t_segment segment[NMAXSEGMENTS])
 {
-    int i, nactive = 0;
-    double fx = 0.0, fy = 0.0;
+    int i, nactive = 0, group;
+    double fx[2] = {0.0, 0.0}, fy[2] = {0.0, 0.0}, x, y, padding = 3.0*MU, mass2 = SEGMENTS_MASS*TWO_CIRCLES_RADIUS_RATIO;
     
+    for (group=0; group<2; group++)
+    {
+        fx[group] = 0.0;
+        fy[group] = 0.0;
+    }
     for (i=0; i<nsegments; i++) if (segment[i].active)
     {
-        fx += segment[i].fx;
-        fy += segment[i].fy;
+        group = segment[i].group;
+        fx[group] += segment[i].fx;
+        fy[group] += segment[i].fy;
         nactive++;
+        if (BOUNDARY_COND == BC_SCREEN) /* add force from simulation boundary */
+        {
+            x = 0.5*(segment[i].x1 + segment[i].x2);
+            y = 0.5*(segment[i].y1 + segment[i].y2);
+            if (x < XMIN + padding) fx[group] += KSPRING_BOUNDARY*(XMIN + padding - x);
+            else if (x > XMAX - padding) fx[group] -= KSPRING_BOUNDARY*(x - XMAX + padding);
+            if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y);
+            else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding);
+        }
+        if (group == 0) fy[group] -= GRAVITY*SEGMENTS_MASS;
+        else fy[group] -= GRAVITY*mass2;
     }
-    if (nactive > 0)
+    if (nactive > 0) for (group=0; group<2; group++)
     {
-        fx = fx/(double)nactive;
-        fy = fy/(double)nactive;
+        fx[group] = fx[group]/(double)nactive;
+        fy[group] = fy[group]/(double)nactive;
     }
-    if (FLOOR_FORCE)
+    if (FLOOR_FORCE) 
     {   
-        if (fx > FMAX) fx = FMAX;
-        else if (fx < -FMAX) fx = -FMAX;
-        if (fy > FMAX) fy = FMAX;
-        else if (fy < -FMAX) fy = -FMAX;
+        if (fx[0] > FMAX) fx[0] = FMAX;
+        else if (fx[0] < -FMAX) fx[0] = -FMAX;
+        if (fy[0] > FMAX) fy[0] = FMAX;
+        else if (fy[0] < -FMAX) fy[0] = -FMAX;
+    }
+    vxsegments[0] += fx[0]*DT_PARTICLE/SEGMENTS_MASS;
+    vysegments[0] += fy[0]*DT_PARTICLE/SEGMENTS_MASS;
+    xsegments[0] += vxsegments[0]*DT_PARTICLE;
+    ysegments[0] += vysegments[0]*DT_PARTICLE;
+    if (TWO_OBSTACLES)
+    {
+        if (FLOOR_FORCE)
+        {
+            if (fx[1] > FMAX) fx[1] = FMAX;
+            else if (fx[1] < -FMAX) fx[1] = -FMAX;
+            if (fy[1] > FMAX) fy[1] = FMAX;
+            else if (fy[1] < -FMAX) fy[1] = -FMAX;
+        }
+        vxsegments[1] += fx[1]*DT_PARTICLE/mass2;
+        vysegments[1] += fy[1]*DT_PARTICLE/mass2;
+        xsegments[1] += vxsegments[1]*DT_PARTICLE;
+        ysegments[1] += vysegments[1]*DT_PARTICLE;
     }
-    vxsegments += fx*DT_PARTICLE/(SEGMENTS_MASS);
-    vysegments += fy*DT_PARTICLE/(SEGMENTS_MASS);
-    xsegments += vxsegments*DT_PARTICLE;
-    ysegments += vysegments*DT_PARTICLE;
 
     /* to avoid numerical instabilities */
-    if (xsegments + 1.0 > BCXMAX) 
+    for (group=0; group<2; group++) if (xsegments[group] + 1.0 > BCXMAX) 
     {
-        xsegments = BCXMAX - 1.0;
-        vxsegments = 0.0;
+        xsegments[group] = BCXMAX - 1.0;
+        vxsegments[group] = 0.0;
     }
         
     translate_segments(segment, xsegments, ysegments);
@@ -783,7 +876,7 @@ void animation()
             thermostat_on = thermostat_schedule(i);
             printf("Termostat: %i\n", thermostat_on);
         }
-        if ((DEACTIVATE_SEGMENT)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME))
+        if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME))
             segment[nsegments-1].active = 0;
         
         blank();
@@ -793,6 +886,7 @@ void animation()
         pright = 0.0;
         if (RECORD_PRESSURES) for (j=0; j<N_PRESSURES; j++) pressure[j] = 0.0;
         
+//         printf("evolving particles\n");
         for(n=0; n<NVID; n++)
         {
             if (MOVE_OBSTACLE) 
@@ -853,10 +947,13 @@ void animation()
                 if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary);
                 else xwall = 0.0;
             }
-            if (MOVE_BOUNDARY) evolve_segments(segment);
+            if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment);
         } /* end of for (n=0; n<NVID; n++) */
+//         printf("evolved particles\n");
+
         
-        if (MOVE_BOUNDARY) printf("segments position (%.3lg, %.3lg), speed (%.3lg, %.3lg)\n", xsegments, ysegments, vxsegments, vysegments);
+        if (MOVE_BOUNDARY) 
+            printf("segments position (%.3lg, %.3lg), speed (%.3lg, %.3lg)\n", xsegments[0], ysegments[0], vxsegments[0], vysegments[0]);
         
 //         if ((PARTIAL_THERMO_COUPLING)) 
         if ((PARTIAL_THERMO_COUPLING)&&(i>N_T_AVERAGE)) 
@@ -918,11 +1015,11 @@ void animation()
 //         printf("Max number of neighbours: %i\n", max_nb); 
         
         if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-        draw_particles(particle, PLOT);
+        draw_particles(particle, PLOT, beta);
         draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
 
         /* add a particle */
-        if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
+        if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
             nadd_particle = add_particles(particle, px, py, nadd_particle);
         
         /* case of reaction-diffusion equation */
@@ -945,7 +1042,8 @@ void animation()
         
         if (PRINT_OMEGA) print_omega(angular_speed); 
         else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
-
+        else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments);
+        
 	glutSwapBuffers();
 
 
@@ -965,7 +1063,7 @@ void animation()
             if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
             {
                 if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-                draw_particles(particle, PLOT_B);
+                draw_particles(particle, PLOT_B, beta);
                 draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
                 print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
                          fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
@@ -973,6 +1071,7 @@ void animation()
                 else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
                 if (PRINT_OMEGA) print_omega(angular_speed); 
                 else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
+                else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments);
                 glutSwapBuffers();
                 save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
                 counter++;
@@ -995,7 +1094,7 @@ void animation()
         if (DOUBLE_MOVIE) 
         {
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT); 
+            draw_particles(particle, PLOT, beta); 
             draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
             print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
                          fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
@@ -1003,13 +1102,14 @@ void animation()
             else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
             if (PRINT_OMEGA) print_omega(angular_speed); 
             else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
+            else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments);
             glutSwapBuffers();
         }
         for (i=0; i<MID_FRAMES; i++) save_frame_lj();
         if (DOUBLE_MOVIE) 
         {
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT_B);
+            draw_particles(particle, PLOT_B, beta);
             draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
             print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
                          fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
@@ -1017,6 +1117,7 @@ void animation()
             else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
             if (PRINT_OMEGA) print_omega(angular_speed); 
             else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
+            else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments);
             glutSwapBuffers();
         }
         for (i=0; i<END_FRAMES; i++) save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
diff --git a/mangrove.c b/mangrove.c
index 100defa..e27c8f6 100644
--- a/mangrove.c
+++ b/mangrove.c
@@ -949,7 +949,7 @@ void animation()
     blank();
     glColor3f(0.0, 0.0, 0.0);
     draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
-    draw_billiard();
+    draw_billiard(0, 1.0);
 
     glutSwapBuffers();
 
@@ -1220,7 +1220,7 @@ void animation()
 
         
         printf("Drawing billiard\n");
-        draw_billiard();
+        draw_billiard(0, 1.0);
 
 	glutSwapBuffers();
 
diff --git a/particle_billiard.c b/particle_billiard.c
index 1f09a68..457adec 100644
--- a/particle_billiard.c
+++ b/particle_billiard.c
@@ -29,14 +29,15 @@
 #include <unistd.h>
 #include <sys/types.h>
 #include <tiffio.h>     /* Sam Leffler's libtiff library. */
+#include <time.h>
 
 #define MOVIE 0         /* set to 1 to generate movie */
 
-#define WINWIDTH 	1280  /* window width */
+#define WINWIDTH 	720  /* window width */
 #define WINHEIGHT 	720   /* window height */
 
-#define XMIN -2.0
-#define XMAX 2.0	/* x interval */
+#define XMIN -1.125
+#define XMAX 1.125	/* x interval */
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
@@ -47,22 +48,20 @@
 #define B_DOMAIN 30     /* choice of domain shape */
 
 #define CIRCLE_PATTERN 1    /* pattern of circles */
-#define POLYLINE_PATTERN 6  /* pattern of polyline */
+#define POLYLINE_PATTERN 8  /* pattern of polyline */
 
 #define ABSORBING_CIRCLES 1 /* set to 1 for circular scatterers to be absorbing */
 
 #define NMAXCIRCLES 100000     /* total number of circles (must be at least NCX*NCY for square grid) */
 #define NMAXPOLY 100000        /* total number of sides of polygonal line */   
-// #define NCX 10            /* number of circles in x direction */
-// #define NCY 10            /* number of circles in y direction */
 #define NCX 30            /* number of circles in x direction */
 #define NCY 20            /* number of circles in y direction */
 #define NPOISSON 500        /* number of points for Poisson C_RAND_POISSON arrangement */
 #define NGOLDENSPIRAL 2000  /* max number of points for C_GOLDEN_SPIRAL arrandement */
 #define SDEPTH 1            /* Sierpinski gastket depth */
 
-#define LAMBDA 0.8	/* parameter controlling shape of domain */
-#define MU 0.015        /* second parameter controlling shape of billiard */
+#define LAMBDA 1.5	/* parameter controlling shape of domain */
+#define MU 0.01          /* second parameter controlling shape of billiard */
 #define FOCI 1          /* set to 1 to draw focal points of ellipse */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
@@ -80,16 +79,17 @@
 #define LMAX 0.01       /* minimal segment length triggering resampling */ 
 #define DMIN 0.02       /* minimal distance to boundary for triggering resampling */ 
 #define CYCLE 1         /* set to 1 for closed curve (start in all directions) */
-#define SHOWTRAILS 1    /* set to 1 to keep trails of the particles */
+#define SHOWTRAILS 0    /* set to 1 to keep trails of the particles */
 #define SHOWZOOM 1      /* set to 1 to show zoom on specific area */
 #define PRINT_PARTICLE_NUMBER 1 /* set to 1 to print number of particles */
 #define PRINT_COLLISION_NUMBER 0 /* set to 1 to print number of collisions */
 #define TEST_ACTIVE 1   /* set to 1 to test whether particle is in billiard */
 
-#define NSTEPS 5000      /* number of frames of movie */
-#define TIME 1500         /* time between movie frames, for fluidity of real-time simulation */ 
-#define DPHI 0.00001    /* integration step */
-#define NVID 150         /* number of iterations between images displayed on screen */
+#define NSTEPS 1000      /* number of frames of movie */
+#define TIME 2500         /* time between movie frames, for fluidity of real-time simulation */ 
+#define DPHI 0.00001     /* integration step */
+#define NVID 100         /* number of iterations between images displayed on screen */
+#define END_FRAMES 25    /* number of still frames at the end of the movie */
 
 /* Decreasing TIME accelerates the animation and the movie                               */
 /* For constant speed of movie, TIME*DPHI should be kept constant                        */
@@ -99,14 +99,14 @@
 
 /* Colors and other graphical parameters */
 
-#define COLOR_PALETTE 11     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 17     /* Color palette, see list in global_pdes.c  */
 
-#define NCOLORS 32       /* number of colors */
+#define NCOLORS 64       /* number of colors */
 #define COLORSHIFT 0     /* hue of initial color */ 
 #define RAINBOW_COLOR 1  /* set to 1 to use different colors for all particles */
 #define FLOWER_COLOR 0   /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */
 #define NSEG 100         /* number of segments of boundary */
-#define LENGTH 0.05       /* length of velocity vectors */
+#define LENGTH 0.03       /* length of velocity vectors */
 #define BILLIARD_WIDTH 2    /* width of billiard */
 #define PARTICLE_WIDTH 3    /* width of particles */
 #define FRONT_WIDTH 3       /* width of wave front */
@@ -120,7 +120,7 @@
 #define PAUSE 200       /* number of frames after which to pause */
 #define PSLEEP 2         /* sleep time during pause */
 #define SLEEP1  1        /* initial sleeping time */
-#define SLEEP2  1000      /* final sleeping time */
+#define SLEEP2  1       /* final sleeping time */
 
 
 #include "global_particles.c"
@@ -437,10 +437,18 @@ void draw_config_showtrails(int color[NPARTMAX], double *configs[NPARTMAX], int
     }
     if (DRAW_BILLIARD) draw_billiard();    
     
-    if (SHOWZOOM) 
-    {
-        draw_zoom(color, configs, active, x_target, y_target, 0.1, 1.65, 0.75, 0.3, 0);
-        draw_zoom(color, configs, active, x_shooter, y_shooter, 0.1, -1.65, 0.75, 0.3, 1);    
+    if (SHOWZOOM) switch (POLYLINE_PATTERN) {
+        case (P_TOKA_PRIME):
+        {
+            draw_zoom(color, configs, active, x_target, y_target, 0.1, 1.65, 0.75, 0.3, 0);
+            draw_zoom(color, configs, active, x_shooter, y_shooter, 0.1, -1.65, 0.75, 0.3, 1);
+            break;
+        }
+        case (P_TOKA_NONSELF):
+        {
+            draw_zoom(color, configs, active, 0.0, 0.0, 0.1, 1.65, 0.75, 0.3, 0);
+            break;
+        }
     }
 //     if (SHOWZOOM) draw_zoom(color, configs, active, 0.95, 0.0, 0.1);
 }
@@ -547,10 +555,18 @@ void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPAR
     }
     if (DRAW_BILLIARD) draw_billiard();    
     
-    if (SHOWZOOM) 
-    {
-        draw_zoom(color, configs, active, x_target, y_target, 0.1, 1.65, 0.75, 0.3, 0);
-        draw_zoom(color, configs, active, x_shooter, y_shooter, 0.1, -1.65, 0.75, 0.3, 1);     
+    if (SHOWZOOM) switch (POLYLINE_PATTERN) {
+        case (P_TOKA_PRIME):
+        {
+            draw_zoom(color, configs, active, x_target, y_target, 0.1, 1.65, 0.75, 0.3, 0);
+            draw_zoom(color, configs, active, x_shooter, y_shooter, 0.1, -1.65, 0.75, 0.3, 1);
+            break;
+        }
+        case (P_TOKA_NONSELF):
+        {
+            draw_zoom(color, configs, active, 0.0, 0.0, 0.1, 0.82, 0.82, 0.25, 0);
+            break;
+        }
     }
 //     if (SHOWZOOM) draw_zoom(color, configs, active, 0.95, 0.0, 0.1);
 }
@@ -668,7 +684,7 @@ void animation()
     
 //     init_drop_config(-1.0 + 0.3*sqrt(2.0), -1.0 + 0.5*sqrt(2.0), 0.0, DPI, configs);
 
-    init_line_config(0.0, 0.3, 0.0, 0.9, 0.0, configs);   
+//     init_line_config(0.0, 0.0, 0.0, 0.9, 0.0, configs);   
 //     init_drop_config(-0.95, 0.0, -0.103 + DPI/15.0, -0.1 + DPI/15.0, configs);
     
     /* find long trajectory */
@@ -697,7 +713,9 @@ void animation()
 //     init_drop_config(x_shooter, y_shooter, alphamax, alphamax + DPI, configs);
     
     
-//     init_drop_config(-0.95, 0.0, 1.0, 1.0 + DPI, configs);    
+//     init_drop_config(-0.5, 0.0, 0.2, 0.4, configs);    
+    init_drop_config(0.0, 0.0, 0.0, DPI, configs);    
+    
 //     init_drop_config(-1.3, -0.1, 0.0, DPI, configs);    
 //     init_drop_config(1.4, 0.1, 0.0, DPI, configs);    
 //     init_drop_config(0.5, 0.5, -1.0, 1.0, configs);    
@@ -759,8 +777,8 @@ void animation()
     
     /* initialize drops in different colors */
 //     init_partial_drop_config(0.0, 0.0, 0.0, DPI, 0, 2*NPART/5, 0, configs, color, newcolor);
-//     init_partial_drop_config(0.0, 0.8, 0.0, DPI, 2*NPART/5, 4*NPART/5, 10, configs, color, newcolor);
-//     init_partial_drop_config(1.2, 0.1, 0.0, DPI, 4*NPART/5, NPART, 36, configs, color, newcolor);
+//     init_partial_drop_config(0.0, 0.8, 0.0, DPI, 2*NPART/5, 4*NPART/5, 30, configs, color, newcolor);
+//     init_partial_drop_config(LAMBDA - 0.05, 0.1, 0.0, DPI, 4*NPART/5, NPART, 60, configs, color, newcolor);
   
     for (i=0; i<=NSTEPS; i++)
     {
@@ -798,7 +816,7 @@ void animation()
  
     if (MOVIE) 
     {
-        for (i=0; i<20; i++) save_frame();
+        for (i=0; i<END_FRAMES; i++) save_frame();
         s = system("mv part*.tif tif_part/");
     }
     
@@ -811,6 +829,12 @@ void animation()
 
 void display(void)
 {
+    time_t rawtime;
+    struct tm * timeinfo;
+
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    
     glPushMatrix();
 
     blank();
@@ -827,6 +851,13 @@ void display(void)
     sleep(SLEEP2); 
 
     glPopMatrix();
+    
+    glutDestroyWindow(glutGetWindow());
+
+    printf("Start local time and date: %s", asctime(timeinfo));
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    printf("Current local time and date: %s", asctime(timeinfo));
 }
 
 
@@ -848,4 +879,3 @@ int main(int argc, char** argv)
     return 0;
 }
 
-
diff --git a/rde.c b/rde.c
index 48b9590..cfeeba7 100644
--- a/rde.c
+++ b/rde.c
@@ -44,37 +44,34 @@
 
 /* General geometrical parameters */
 
-// #define WINWIDTH 	1920  /* window width */
-// #define WINHEIGHT 	1000  /* window height */
-// #define NX 480          /* number of grid points on x axis */
-// #define NY 250          /* number of grid points on y axis */
-// // #define NX 1920          /* number of grid points on x axis */
-// // #define NY 1000          /* number of grid points on y axis */
-// 
-// #define XMIN -2.0
-// #define XMAX 2.0	/* x interval */
-// #define YMIN -1.041666667
-// #define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
-
-#define WINWIDTH 	1280  /* window width */
-#define WINHEIGHT 	720   /* window height */
-
-// #define NX 160          /* number of grid points on x axis */
-// #define NY 90          /* number of grid points on y axis */
-#define NX 320          /* number of grid points on x axis */
-#define NY 180          /* number of grid points on y axis */
-
-// #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 NX 720          /* number of grid points on x axis */
-// #define NY 720          /* number of grid points on y axis */
+#define WINWIDTH 	1920  /* window width */
+#define WINHEIGHT 	1000  /* window height */
+#define NX 480          /* number of grid points on x axis */
+#define NY 240          /* 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 YMIN -1.041666667
+#define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
+
+// #define WINWIDTH 	1280  /* window width */
+// #define WINHEIGHT 	720   /* window height */
+// 
+// // #define NX 160          /* number of grid points on x axis */
+// // #define NY 90          /* number of grid points on y axis */
+// #define NX 320          /* number of grid points on x axis */
+// #define NY 180          /* number of grid points on y axis */
+// 
+// // #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 XMIN -2.0
+// #define XMAX 2.0	/* x interval */
+// #define YMIN -1.125
+// #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
 /* Choice of simulated equation */
 
@@ -82,14 +79,14 @@
 #define NFIELDS 2       /* number of fields in reaction-diffusion equation */
 #define NLAPLACIANS 2   /* number of fields for which to compute Laplacian */
 
-#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodiner equation) */
-#define POTENTIAL 1     /* type of potential, see list in global_3d.c  */
+#define ADD_POTENTIAL 1 /* set to 1 to add a potential (for Schrodiner equation) */
+#define POTENTIAL 2     /* type of potential, see list in global_3d.c  */
 
 #define JULIA_SCALE 0.5 /* scaling for Julia sets */
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 9      /* choice of domain shape, see list in global_pdes.c  */
+#define B_DOMAIN 999      /* choice of domain shape, see list in global_pdes.c  */
 
 #define CIRCLE_PATTERN 99    /* pattern of circles, see list in global_pdes.c */
 
@@ -97,11 +94,11 @@
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.4	    /* parameter controlling the dimensions of domain */
-#define MU 0.06	            /* parameter controlling the dimensions of domain */
+#define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
+#define MU 1.0	            /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
-#define APOLY 1.0           /* angle by which to turn polygon, in units of Pi/2 */
-#define MDEPTH 5            /* depth of computation of Menger gasket */
+#define APOLY 0.333333333   /* angle by which to turn polygon, in units of Pi/2 */
+#define MDEPTH 7            /* depth of computation of Menger gasket */
 #define MRATIO 5            /* ratio defining Menger gasket */
 #define MANDELLEVEL 1000      /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0     /* limit value for approximation of Mandelbrot set */
@@ -125,38 +122,44 @@
 
 /* Physical patameters of wave equation */
 
-#define DT 0.00000001
+#define DT 0.00000002
 
-#define VISCOSITY 0.1
+#define VISCOSITY 2.0
 
 #define RPSA 0.75         /* parameter in Rock-Paper-Scissors-type interaction */
 #define RPSLZB 0.75       /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */
 
 #define EPSILON 0.8     /* time scale separation */
 #define DELTA 0.1       /* time scale separation */
-#define FHNA 1.0    /* parameter in FHN equation */
-#define FHNC -0.01    /* parameter in FHN equation */
-#define K_HARMONIC 1.5  /* spring constant of harmonic potential */
-#define BZQ 0.0008       /* parameter in BZ equation */
-#define BZF 1.2          /* parameter in BZ equation */
+#define FHNA 1.0        /* parameter in FHN equation */
+#define FHNC -0.01      /* parameter in FHN equation */
+#define K_HARMONIC 0.2  /* spring constant of harmonic potential */
+#define K_COULOMB 0.5   /* constant in Coulomb potential */
+#define BZQ 0.0008      /* parameter in BZ equation */
+#define BZF 1.2         /* parameter in BZ equation */
 
 #define T_OUT 2.0       /* outside temperature */
 #define T_IN 0.0        /* inside temperature */
 #define SPEED 0.0       /* speed of drift to the right */
 
 #define ADD_NOISE 0     /* set to 1 to add noise, set to 2 to add noise in right half */
-#define NOISE_INTENSITY 0.1     /* noise intensity */
+#define NOISE_INTENSITY 0.005      /* noise intensity */
+#define CHANGE_NOISE 1      /* set to 1 to increase noise intensity */
+#define NOISE_FACTOR 40.0   /* factor by which to increase noise intensity */
+#define NOISE_INITIAL_TIME 100  /* initial time during which noise remains constant */
 
 #define CHANGE_VISCOSITY 0      /* set to 1 to change the viscosity in the course of the simulation */
-#define ADJUST_INTSTEP 1        /* set to 1 to decrease integration step when viscosity increases */
+#define ADJUST_INTSTEP 0       /* set to 1 to decrease integration step when viscosity increases */
 #define VISCOSITY_INITIAL_TIME 10  /* initial time during which viscosity remains constant */
 #define VISCOSITY_FACTOR 100.0   /* factor by which to change viscosity */
 #define VISCOSITY_MAX 2.0        /* max value of viscosity beyond which NVID is increased and integration step is decrase, 
                                     for numerical stability */
+                                        
 #define CHANGE_RPSLZB 0         /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */
 #define RPSLZB_CHANGE 0.75      /* factor by which to rpslzb parameter */
 #define RPSLZB_INITIAL_TIME 0   /* initial time during which rpslzb remains constant */
 #define RPSLZB_FINAL_TIME 500   /* final time during which rpslzb remains constant */
+                                      
 
 /* Boundary conditions, see list in global_pdes.c  */
 
@@ -164,13 +167,13 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1500      /* number of frames of movie */
-#define NVID 900           /* number of iterations between images displayed on screen */
+#define NSTEPS 1150      /* number of frames of movie */
+#define NVID 850          /* number of iterations between images displayed on screen */
 #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */
 #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation  */
 #define MAX_DT 0.024     /* maximal value of integration step */
 #define NSEG 100         /* number of segments of boundary */
-#define BOUNDARY_WIDTH 1    /* width of billiard boundary */
+#define BOUNDARY_WIDTH 4    /* width of billiard boundary */
 
 #define PAUSE 100       /* number of frames after which to pause */
 #define PSLEEP 2         /* sleep time during pause */
@@ -178,11 +181,16 @@
 #define SLEEP2  1        /* final sleeping time */
 #define INITIAL_TIME 0  /* initial still time */
 #define MID_FRAMES 50    /* number of still frames between parts of two-part movie */
-#define END_FRAMES 100   /* number of still frames at end of movie */
+#define END_FRAMES 50    /* number of still frames at end of movie */
 #define FADE 1           /* set to 1 to fade at end of movie */
 
+/* Visualisation */
+
 #define PLOT_3D 1    /* controls whether plot is 2D or 3D */
 
+#define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
+#define ROTATE_ANGLE 360.0  /* total angle of rotation during simulation */
+
 /* Plot type - color scheme */
 
 #define CPLOT 30
@@ -190,7 +198,9 @@
 
 /* Plot type - height of 3D plot */
 
-#define ZPLOT 30      /* z coordinate in 3D plot */
+// #define ZPLOT 30     /* z coordinate in 3D plot */
+// #define ZPLOT_B 32    /* z coordinate in second 3D plot */
+#define ZPLOT 30     /* z coordinate in 3D plot */
 #define ZPLOT_B 30    /* z coordinate in second 3D plot */
 
 #define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */
@@ -198,13 +208,17 @@
 #define NON_DIRICHLET_BC 0      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
 #define WRAP_ANGLE 1            /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */
 #define FADE_IN_OBSTACLE 0      /* set to 1 to fade color inside obstacles */
-#define DRAW_OUTSIDE_GRAY 1     /* experimental - draw outside of billiard in gray */
+#define DRAW_OUTSIDE_GRAY 0     /* experimental - draw outside of billiard in gray */
+#define ADD_POTENTIAL_TO_Z 1    /* set to 1 to add the external potential to z-coordinate of plot */
+#define ADD_POT_CONSTANT 0.5    /* constant in front of added potential */
 
 #define PLOT_SCALE_ENERGY 0.05      /* vertical scaling in energy plot */
 
-#define PRINT_TIME 0    /* set to 1 to print running time */
-#define PRINT_VISCOSITY 0    /* set to 1 to print viscosity */
+#define PRINT_TIME 0        /* set to 1 to print running time */
+#define PRINT_VISCOSITY 0   /* set to 1 to print viscosity */
 #define PRINT_RPSLZB 0      /* set to 1 to print rpslzb parameter */
+#define PRINT_PROBABILITIES 0   /* set to 1 to print probabilities (for Ehrenfest urn configuration) */
+#define PRINT_NOISE 0       /* set to 1 to print noise intensity */
 
 #define DRAW_FIELD_LINES 0  /* set to 1 to draw field lines */
 #define FIELD_LINE_WIDTH 1  /* width of field lines */
@@ -223,22 +237,23 @@
 
 /* Color schemes, see list in global_pdes.c  */
 
-#define COLOR_PALETTE 11     /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 17     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 14     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 10     /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* black background */
 
 #define COLOR_SCHEME 3   /* choice of color scheme */
 
-#define COLOR_PHASE_SHIFT 1.0   /* phase shift of color scheme, in units of Pi */
+#define COLOR_PHASE_SHIFT 0.0   /* phase shift of color scheme, in units of Pi */
 
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
-#define SLOPE 1.0       /* sensitivity of color on wave amplitude */
-#define VSCALE_AMPLITUDE 20.0      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define SLOPE 1.0        /* sensitivity of color on wave amplitude */
+#define VSCALE_AMPLITUDE 7.5      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
 #define CURL_SCALE 0.000015   /* scaling factor for curl representation */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
-#define SLOPE_SCHROD_LUM 50.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
+#define SLOPE_SCHROD_LUM 15.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
+#define MIN_SCHROD_LUM 0.075      /* minimal luminosity in color scheme Z_ARGUMENT*/
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
 #define COLORDRIFT 0.0   /* how much the color hue drifts during the whole simulation */
@@ -250,9 +265,9 @@
 #define LOG_SCALE 1.0     /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.0 
 
-#define DRAW_COLOR_SCHEME 0     /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 0.55    /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 12.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 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 3.0   /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 
 /* only for compatibility with wave_common.c */
@@ -273,20 +288,23 @@ double u_3d[2] = {0.75, -0.45};     /* projections of basis vectors for REP_AXO_
 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] = {8.0, 8.0, 8.0};    /* location of observer for REP_PROJ_3D representation */ 
+double observer[3] = {8.0, 8.0, 6.0};    /* location of observer for REP_PROJ_3D representation */ 
+int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
-#define Z_SCALING_FACTOR 0.75    /* 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 Z_SCALING_FACTOR 0.75   /* 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.05          /* overall y shift for REP_PROJ_3D representation */
+#define XSHIFT_3D -0.1           /* overall x shift for REP_PROJ_3D representation */
+#define YSHIFT_3D 0.2          /* overall y shift for REP_PROJ_3D representation */
+
 
 /* For debugging purposes only */
 #define FLOOR 1         /* set to 1 to limit wave amplitude to VMAX */
-#define VMAX 10.0       /* max value of wave amplitude */
+#define VMAX 2.0        /* max value of wave amplitude */
 
 #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT)    /* to save computing time, to be improved */
 #define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)))
+#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE))
 
 #include "global_pdes.c"
 #include "sub_wave.c"
@@ -301,7 +319,7 @@ double observer[3] = {8.0, 8.0, 8.0};    /* location of observer for REP_PROJ_3D
 double potential(int i, int j)
 /* compute potential (e.g. for Schrödinger equation) */
 {
-    double x, y, xy[2];
+    double x, y, xy[2], r, small = 2.0e-1, kx, ky;
     
     ij_to_xy(i, j, xy);
     x = xy[0];
@@ -312,6 +330,19 @@ double potential(int i, int j)
         {
             return (K_HARMONIC*(x*x + y*y));
         }
+        case (POT_COULOMB):
+        {
+//             r = module2(x, y);
+            r = sqrt(x*x + y*y + small*small);
+//             if (r < small) r = small;
+            return (-K_COULOMB/r);
+        }
+        case (POT_PERIODIC):
+        {
+            kx = 4.0*DPI/(XMAX - XMIN);
+            ky = 2.0*DPI/(YMAX - YMIN);
+            return(-K_HARMONIC*cos(kx*x)*cos(ky*y));
+        }
         default:
         {
             return(0.0);
@@ -321,7 +352,7 @@ double potential(int i, int j)
 
 
 void initialize_potential(double potential_field[NX*NY])
-/* initialize the potential field, e.f. for the Schrödinger equation */
+/* initialize the potential field, e.g. for the Schrödinger equation */
 {
     int i, j;
     
@@ -455,10 +486,10 @@ void print_level(int level)
 
 
 
-void print_parameters(double time, short int left, double viscosity)
+void print_parameters(t_rde rde[NX*NY], short int xy_in[NX*NY], double time, short int left, double viscosity, double noise)
 {
     char message[100];
-    double density, hue, rgb[3], logratio, y, pos[2];
+    double density, hue, rgb[3], logratio, x, y, pos[2], probas[2];
     static double xbox, xtext, boxwidth, boxheight;
     static int first = 1;
     
@@ -498,17 +529,40 @@ void print_parameters(double time, short int left, double viscosity)
         first = 0;
     }
     
-    y = YMAX - 0.1;
-    erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0);
-    glColor3f(1.0, 1.0, 1.0);
-    if (PRINT_TIME) sprintf(message, "Time %.3f", time);
-    else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity);
-    else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb);
-    if (PLOT_3D) write_text(xtext, y, message);
+    if (PRINT_PROBABILITIES)
+    {
+        compute_probabilities(rde, xy_in, probas);
+        printf("pleft = %.3lg, pright = %.3lg\n", probas[0], probas[1]);
+        
+        x = XMIN + 0.15*(XMAX - XMIN);
+        y = YMIN + 0.3*(YMAX - YMIN);
+        erase_area_hsl(x, y, boxwidth, boxheight, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "Proba %.3f", probas[0]);
+        write_text(x, y, message);
+        
+        x = XMIN + 0.72*(XMAX - XMIN);
+        y = YMIN + 0.68*(YMAX - YMIN);
+        erase_area_hsl(x, y, boxwidth, boxheight, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "Proba %.3f", probas[1]);
+        write_text(x, y, message);
+    }
     else
     {
-        xy_to_pos(xtext, y, pos);
-        write_text(pos[0], pos[1], message);
+        y = YMAX - 0.1;
+        erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        if (PRINT_TIME) sprintf(message, "Time %.3f", time);
+        else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity);
+        else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb);
+        else if (PRINT_NOISE) sprintf(message, "noise %.3f", noise);
+        if (PLOT_3D) write_text(xtext, y, message);
+        else
+        {
+            xy_to_pos(xtext, y, pos);
+            write_text(pos[0], pos[1], message);
+        }
     }
 }
 
@@ -523,6 +577,18 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
         draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
 }
 
+double noise_schedule(int i)
+{
+    double ratio;
+    
+    if (i < NOISE_INITIAL_TIME) return (NOISE_INTENSITY);
+    else 
+    {
+        ratio = (double)(i - NOISE_INITIAL_TIME)/(double)(NSTEPS - NOISE_INITIAL_TIME);
+        return (NOISE_INTENSITY*(1.0 + ratio*(NOISE_FACTOR - 1.0)));
+    }
+}
+
 
 double viscosity_schedule(int i)
 {
@@ -550,10 +616,33 @@ double rpslzb_schedule(int i)
 }
 
 
+void viewpoint_schedule(int i)
+/* change position of observer */
+{
+    int j;
+    double angle, ca, sa;
+    static double observer_initial[3];
+    static int first = 1;
+    
+    if (first)
+    {
+        for (j=0; j<3; j++) observer_initial[j] = observer[j];
+        first = 0;
+    }
+    
+    angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
+    ca = cos(angle);
+    sa = sin(angle);
+    observer[0] = ca*observer_initial[0] - sa*observer_initial[1];
+    observer[1] = sa*observer_initial[0] + ca*observer_initial[1];
+    printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
+}
+
 
 void animation()
 {
-    double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, intstep0, viscosity_printed, fade_value;
+    double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, 
+        intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY;
     double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field;
     short int *xy_in;
     int i, j, k, s, nvid, field;
@@ -570,7 +659,11 @@ void animation()
     xy_in = (short int *)malloc(NX*NY*sizeof(short int));
     rde = (t_rde *)malloc(NX*NY*sizeof(t_rde));
     
-    if (ADD_POTENTIAL) potential_field = (double *)malloc(NX*NY*sizeof(double));
+    if (ADD_POTENTIAL) 
+    {
+        potential_field = (double *)malloc(NX*NY*sizeof(double));
+        initialize_potential(potential_field);
+    }
 
     npolyline = init_polyline(MDEPTH, polyline);
     for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
@@ -591,7 +684,7 @@ void animation()
 //     init_random(0.5, 0.4, phi, xy_in);
 //     init_random(0.0, 0.4, phi, xy_in);
 //     init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
-    init_coherent_state(-0.75, 0.0, 5.0, -7.0, 0.2, phi, xy_in);
+    init_coherent_state(-0.7, 0.0, 3.5, 0.0, 0.15, phi, xy_in);
     
     init_cfield_rde(phi, xy_in, CPLOT, rde, 0);
     if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0);
@@ -609,9 +702,9 @@ void animation()
     glutSwapBuffers();
     
     printf("Drawing wave\n");
-    draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+    draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
 //     draw_billiard();
-    if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, VISCOSITY);
+    if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, VISCOSITY, noise);
     if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
 
     glutSwapBuffers();
@@ -638,8 +731,14 @@ void animation()
         }
         if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i);
         
+        if (ROTATE_VIEW) 
+        {
+            viewpoint_schedule(i - INITIAL_TIME);
+            reset_view = 1;
+        }
+        
         printf("Drawing wave %i\n", i);
-        draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+        draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
         
 //         nvid = (int)((double)NVID*(1.0 + (ACCELERATION_FACTOR - 1.0)*(double)i/(double)NSTEPS));
         /* increase integration step */
@@ -660,22 +759,36 @@ void animation()
         
         if (ADD_NOISE == 1) 
         {
-            for (field=0; field<=NFIELDS; field++)
+//             #pragma omp parallel for private(field,j,k)
+            for (field=0; field<NFIELDS; field++)
                 for (j=0; j<NX; j++) 
                     for (k=0; k<NY; k++)
                         phi[field][j*NY+k] += sqrintstep*NOISE_INTENSITY*gaussian();
         }
         else if (ADD_NOISE == 2) 
         {
-            for (field=0; field<=NFIELDS; field++)
-                for (j=NX/2; j<NX; j++) 
-                    for (k=0; k<NY; k++)
-                        phi[field][j*NY+k] += sqrintstep*NOISE_INTENSITY*gaussian();
+            if (CHANGE_NOISE)
+            {
+                noise = noise_schedule(i);
+//                 #pragma omp parallel for private(field,j,k)
+                for (field=0; field<NFIELDS; field++)
+                    for (j=NX/2; j<NX; j++) 
+                        for (k=0; k<NY; k++)
+                            phi[field][j*NY+k] += sqrintstep*noise*gaussian();
+            }
+            else
+            {
+//                 #pragma omp parallel for private(field,j,k)
+                for (field=0; field<NFIELDS; field++)
+                    for (j=NX/2; j<NX; j++) 
+                        for (k=0; k<NY; k++)
+                            phi[field][j*NY+k] += sqrintstep*NOISE_INTENSITY*gaussian();
+            }
         }
         time += nvid*intstep;
         
 //         draw_billiard();
-        if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
+        if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
         if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); 
         
 //         print_level(MDEPTH);
@@ -693,9 +806,9 @@ void animation()
             
             if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
             {
-                draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
+                draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
 //                 draw_billiard();
-                if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
+                if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);  
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
@@ -719,9 +832,9 @@ void animation()
     {
         if (DOUBLE_MOVIE) 
         {
-            draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+            draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
 //             draw_billiard();
-            if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
+            if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);   
             glutSwapBuffers();
             
@@ -729,14 +842,14 @@ void animation()
             else for (i=0; i<MID_FRAMES; i++) 
             {
                 fade_value = 1.0 - (double)i/(double)MID_FRAMES;
-                draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
+                draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
 //                 draw_billiard();
-                if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
+                if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);   
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + i + 1);
             }
-            draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
+            draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); 
             glutSwapBuffers();
             
@@ -744,7 +857,7 @@ void animation()
             else for (i=0; i<END_FRAMES; i++) 
             {
                 fade_value = 1.0 - (double)i/(double)END_FRAMES;
-                draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
+                draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);   
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@@ -756,7 +869,7 @@ void animation()
             else for (i=0; i<END_FRAMES; i++) 
             {
                 fade_value = 1.0 - (double)i/(double)END_FRAMES;
-                draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
+                draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); 
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + 1 + counter + i);
diff --git a/schrodinger.c b/schrodinger.c
index ef78b25..e7921da 100644
--- a/schrodinger.c
+++ b/schrodinger.c
@@ -638,7 +638,7 @@ void animation()
         
         for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
         
-        draw_billiard();
+        draw_billiard(0, 1.0);
         
         if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT, COLORBAR_RANGE); 
         
diff --git a/sub_lj.c b/sub_lj.c
index 1cb2eb5..341f498 100644
--- a/sub_lj.c
+++ b/sub_lj.c
@@ -8,6 +8,7 @@
 // #define HUE_TYPE0 300.0     /* hue of particles of type 0 */
 // #define HUE_TYPE1 90.0      /* hue of particles of type 1 */
 
+
 int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression)
 {
   TIFF *file;
@@ -782,6 +783,18 @@ void init_particle_config(t_particle particles[NMAXCIRCLES])
                 }
             break;
         }
+        case (C_POOL_TABLE):
+        {
+            for (i=1; i<6; i++) for (j=0; j<i; j++)
+            {
+                particles[ncircles].xc = INITXMIN + (double)i*0.25*(INITXMAX - INITXMIN);
+                particles[ncircles].yc = 0.5*(INITYMIN + INITYMAX) + ((double)j - 0.5*(double)(i-1))*0.25*(INITYMAX - INITYMIN);
+                particles[ncircles].radius = MU;
+                particles[ncircles].active = 1;
+                ncircles += 1;
+            }
+            break;
+        }
         case (C_ONE):
         {
             particles[ncircles].xc = 0.0;
@@ -836,11 +849,28 @@ void init_people_config(t_person people[NMAXCIRCLES])
     
 }
 
+
+void add_obstacle(double x, double y, double radius, t_obstacle obstacle[NMAXOBSTACLES])
+/* add a circular obstacle to obstacle configuration */
+{
+    if (nobstacles + 1 < NMAXOBSTACLES)
+    {
+        obstacle[nobstacles].xc = x;
+        obstacle[nobstacles].yc = y;
+        obstacle[nobstacles].radius = radius;
+        obstacle[nobstacles].active = 1;
+        
+        nobstacles++;
+    }
+    else printf("Warning: NMAXOBSTACLES should be increased\n");
+}
+
+
 void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
 /* initialise circular obstacle configuration */
 {
     int i, j, n; 
-    double x, y, dx, dy;
+    double x, y, dx, dy, width, lpocket, xmid = 0.5*(BCXMIN + BCXMAX), radius;
     
     switch (OBSTACLE_PATTERN) {
         case (O_CORNERS):
@@ -890,19 +920,130 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
             nobstacles = n;
             break;
         }
+        case (O_POOL_TABLE):
+        {
+            lpocket = 0.1;
+            width = 0.5*MU;
+            radius = 2.0*width;
+            
+            add_obstacle(BCXMIN + lpocket, BCYMIN - width, radius, obstacle);
+            add_obstacle(xmid - lpocket, BCYMIN - width, radius, obstacle);
+
+            add_obstacle(xmid + lpocket, BCYMIN - width, radius, obstacle);
+            add_obstacle(BCXMAX - lpocket, BCYMIN - width, radius, obstacle);
+            
+            add_obstacle(BCXMAX + width, BCYMIN + lpocket, radius, obstacle);
+            add_obstacle(BCXMAX + width, BCYMAX - lpocket, radius, obstacle);
+            
+            add_obstacle(BCXMIN + lpocket, BCYMAX + width, radius, obstacle);
+            add_obstacle(xmid - lpocket, BCYMAX + width, radius, obstacle);
+
+            add_obstacle(xmid + lpocket, BCYMAX + width, radius, obstacle);
+            add_obstacle(BCXMAX - lpocket, BCYMAX + width, radius, obstacle);
+            
+            add_obstacle(BCXMIN - width, BCYMIN + lpocket, radius, obstacle);
+            add_obstacle(BCXMIN - width, BCYMAX - lpocket, radius, obstacle);
+            
+            break;
+        }
         default: 
         {
             printf("Function init_obstacle_config not defined for this pattern \n");
         }
     }
 }
+
+void add_rectangle_to_segments(double x1, double y1, double x2, double y2, t_segment segment[NMAXSEGMENTS])
+/* add four segements forming a rectangle to linear obstacle configuration */
+{
+    int i, n = nsegments; 
+    
+    if (nsegments + 4 < NMAXSEGMENTS)
+    {
+        segment[n].x1 = x1;
+        segment[n].y1 = y1;
+        segment[n].x2 = x2;
+        segment[n].y2 = y1;
+        
+        segment[n+1].x1 = x2;
+        segment[n+1].y1 = y1;
+        segment[n+1].x2 = x2;
+        segment[n+1].y2 = y2;
+        
+        segment[n+2].x1 = x2;
+        segment[n+2].y1 = y2;
+        segment[n+2].x2 = x1;
+        segment[n+2].y2 = y2;
+        
+        segment[n+3].x1 = x1;
+        segment[n+3].y1 = y2;
+        segment[n+3].x2 = x1;
+        segment[n+3].y2 = y1;
+        
+        for (i=0; i<4; i++) segment[n+i].concave = 1;
+        nsegments += 4;
+    }
+    else printf("Warning: NMAXSEGMENTS too small\n");
+}
+
+void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS])
+/* add four segments forming a rectangle, specified by two adjacent corners and width */
+{
+    double tx, ty, ux, uy, norm, x3, y3, x4, y4;
+    int i, n = nsegments; 
+    
+    tx = x2 - x1;
+    ty = y2 - y1;
+    norm = module2(tx, ty);
+    tx = tx/norm;
+    ty = ty/norm;
+    x3 = x2 + width*ty;
+    y3 = y2 - width*tx;
+    x4 = x1 + width*ty;
+    y4 = y1 - width*tx;
+    
+    if (nsegments + 4 < NMAXSEGMENTS)
+    {
+        segment[n].x1 = x1;
+        segment[n].y1 = y1;
+        segment[n].x2 = x2;
+        segment[n].y2 = y2;
+        
+        segment[n+1].x1 = x2;
+        segment[n+1].y1 = y2;
+        segment[n+1].x2 = x3;
+        segment[n+1].y2 = y3;
+        
+        segment[n+2].x1 = x3;
+        segment[n+2].y1 = y3;
+        segment[n+2].x2 = x4;
+        segment[n+2].y2 = y4;
+        
+        segment[n+3].x1 = x4;
+        segment[n+3].y1 = y4;
+        segment[n+3].x2 = x1;
+        segment[n+3].y2 = y1;
+        
+        for (i=0; i<4; i++) segment[n+i].concave = 1;
+        nsegments += 4;
+    }
+    else printf("Warning: NMAXSEGMENTS too small\n");
+}
  
+double nozzle_width(double x, double width)
+/* width of bell-shaped nozzle */
+{
+    double lam  = 0.5*LAMBDA;
+    
+    if (x >= 0.0) return(width);
+    else return(sqrt(width*width - 0.5*x));
+}
  
 void init_segment_config(t_segment segment[NMAXSEGMENTS])
 /* initialise linear obstacle configuration */
 {
-    int i, cycle = 0, iminus, iplus; 
-    double angle, angle2, dx, width, height, a, b, length;
+    int i, j, cycle = 0, iminus, iplus, nsides, n; 
+    double angle, angle2, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy;
     
     switch (SEGMENT_PATTERN) {
         case (S_RECTANGLE):
@@ -922,7 +1063,11 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             cycle = 1;
             nsegments = 4;
             
-            for (i=0; i<nsegments; i++) segment[i].concave = 0;
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 0;
+                segment[i].group = 0;
+            }
             break;
         }
         case (S_CUP):
@@ -945,7 +1090,12 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             cycle = 1;
             nsegments = 4;
             
-            for (i=0; i<nsegments; i++) segment[i].concave = 0;
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 0;
+                segment[i].group = 0;
+            }
+            
             break;
         }
         case (S_HOURGLASS):
@@ -988,6 +1138,8 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             
             cycle = 1;
             nsegments = 8;
+            for (i=0; i<nsegments; i++) segment[i].group = 0;
+            
             break;
         }
         case (S_PENTA):
@@ -1013,7 +1165,12 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             cycle = 1;
             nsegments = 5;
             
-            for (i=0; i<nsegments; i++) segment[i].concave = 0;
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 0;
+                segment[i].group = 0;
+            }
+            
             break;
         }
         case (S_CENTRIFUGE):
@@ -1041,6 +1198,9 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             
             cycle = 1;
             nsegments = 4*NPOLY;
+            
+            for (i=0; i<nsegments; i++) segment[i].group = 0;
+            
             break;
         }
         case (S_POLY_ELLIPSE):
@@ -1066,6 +1226,214 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             nsegments = 2*NPOLY;
             break;
         }
+        case (S_POOL_TABLE):
+        {
+            width = MU;
+            lpocket = 0.1;
+            
+            add_rectangle_to_segments(BCXMIN + lpocket, BCYMIN, xmid - lpocket, BCYMIN - width, segment);
+            add_rectangle_to_segments(xmid + lpocket, BCYMIN, BCXMAX - lpocket, BCYMIN - width, segment);
+            add_rectangle_to_segments(BCXMAX + width, BCYMIN + lpocket, BCXMAX, BCYMAX - lpocket, segment); 
+            
+            add_rectangle_to_segments(BCXMAX - lpocket, BCYMAX, xmid + lpocket, BCYMAX + width, segment);
+            add_rectangle_to_segments(xmid - lpocket, BCYMAX, BCXMIN + lpocket, BCYMAX + width, segment);
+            add_rectangle_to_segments(BCXMIN - width, BCYMAX - lpocket, BCXMIN, BCYMIN + lpocket, segment); 
+            
+            cycle = 0;
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 0;
+                segment[i].group = 0;
+            }
+            
+            break;
+        }   
+        case (S_CENTRIFUGE_RND):
+        {
+            angle = DPI/(double)NPOLY;
+            nsides = 24;
+            if ((nsides+2)*NPOLY > NMAXSEGMENTS)
+            {
+                printf("Error: NMAXSEGMENTS is too small\n");
+                exit(1);
+            }
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                segment[i*(nsides+2)].x1 = LAMBDA*cos(((double)i + 0.02)*angle);
+                segment[i*(nsides+2)].y1 = LAMBDA*sin(((double)i + 0.02)*angle);
+                segment[i*(nsides+2)].concave = 1;
+                
+                for (j=1; j<=nsides; j++)
+                {
+                    x = (double)j/(double)(nsides+1);
+                    r = 0.5 + sqrt(x*(1.0-x));
+                    angle2 = (double)i*angle + angle*(double)j/(double)(nsides+1);
+                    segment[i*(nsides+2) + j].x1 = r*cos(angle2);
+                    segment[i*(nsides+2) + j].y1 = r*sin(angle2);
+                    segment[i*(nsides+2) + j].concave = 0;
+                }
+                
+                segment[i*(nsides+2) + nsides + 1].x1 = LAMBDA*cos(((double)i + 0.98)*angle);
+                segment[i*(nsides+2) + nsides + 1].y1 = LAMBDA*sin(((double)i + 0.98)*angle);
+                segment[i*(nsides+2) + nsides + 1].concave = 1;
+            }
+            
+            cycle = 1;
+            nsegments = (nsides+2)*NPOLY;
+            for (i=0; i<nsegments; i++) segment[i].group = 0;
+            
+            break;
+        }
+        case (S_CENTRIFUGE_LEAKY):
+        {
+            angle = DPI/(double)NPOLY;
+            nsides = 20;
+            if (2*(nsides+2)*NPOLY > NMAXSEGMENTS)
+            {
+                printf("Error: NMAXSEGMENTS is too small\n");
+                exit(1);
+            }
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                angle2 = (double)i*angle;
+                x1 = LAMBDA*cos(angle2);
+                y1 = LAMBDA*sin(angle2);
+                x2 = 0.7*cos(angle2);
+                y2 = 0.7*sin(angle2);
+                add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment);
+                
+                for (j=0; j<nsides; j++) if (j!=nsides/2)
+                {
+                    x = (double)j/(double)(nsides);
+                    r = 0.5 + sqrt(x*(1.0-x) + 0.04);
+                    if (j < nsides/2) angle2 = (double)i*angle + angle*(double)j/((double)(nsides) - 0.15);
+                    else angle2 = (double)i*angle + angle*(double)j/((double)(nsides) + 0.15);
+                    x1 = r*cos(angle2);
+                    y1 = r*sin(angle2);
+                    
+                    x = (double)(j+1)/(double)(nsides);
+                    r = 0.5 + sqrt(x*(1.0-x) + 0.04);
+                    if (j < nsides/2) angle2 = (double)i*angle + angle*(double)(j+1)/((double)(nsides) - 0.15);
+                    else angle2 = (double)i*angle + angle*(double)(j+1)/((double)(nsides) + 0.15);
+                    x2 = r*cos(angle2);
+                    y2 = r*sin(angle2);
+                    
+                    add_rotated_angle_to_segments(x1, y1, x2, y2, 0.5*MU, segment);
+                }
+            }
+            
+            cycle = 0;
+            
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 0;
+                segment[i].group = 0;
+            }
+            
+            break;
+        }
+        case (S_CIRCLE_EXT):
+        {
+            angle = DPI/(double)NPOLY;
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                segment[i].x1 = SEGMENTS_X0 + LAMBDA*cos(((double)i)*angle);
+                segment[i].y1 = SEGMENTS_Y0 - LAMBDA*sin(((double)i)*angle);
+                segment[i].concave = 1;
+            }
+            
+            cycle = 1;
+            nsegments = NPOLY;
+            
+            for (i=0; i<nsegments; i++) segment[i].group = 0;
+            
+            break;
+        }
+        case (S_ROCKET_NOZZLE):
+        {
+            /* ellipse */
+            for (i=1; i<NPOLY-1; i++)
+            {
+                angle = -PI + (double)i*DPI/(double)NPOLY;
+                x1 = 0.7*LAMBDA*(1.0 + cos(angle));
+                y1 = 0.5*LAMBDA*sin(angle);
+                angle = -PI + (double)(i+1)*DPI/(double)NPOLY;
+                x2 = 0.7*LAMBDA*(1.0 + cos(angle));
+                y2 = 0.5*LAMBDA*sin(angle);                
+                add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
+            }
+             
+            /* compute intersection point of nozzle and ellipse */
+            angle = PI - DPI/(double)NPOLY;
+            nozx = 0.7*LAMBDA*(1.0 + cos(angle));
+            nozy = 0.5*LAMBDA*sin(angle);
+            
+            nsides = 10;
+            dx = LAMBDA/(double)(nsides);
+            
+            /* nozzle */
+            for (i=0; i<nsides; i++)
+            {
+                x1 = -LAMBDA + dx*(double)(i-1);
+                y1 = nozzle_width(x1, nozy);
+                x2 = x1 + dx;
+                y2 = nozzle_width(x2, nozy);
+                add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
+            }
+            add_rotated_angle_to_segments(x2, y2, nozx, nozy, 0.02, segment);
+            for (i=0; i<nsides; i++)
+            {
+                x1 = -LAMBDA + dx*(double)(i-1);
+                y1 = -nozzle_width(x1, nozy);
+                x2 = x1 + dx;
+                y2 = -nozzle_width(x2, nozy);
+                add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
+            }
+            add_rotated_angle_to_segments(x2, y2, nozx, -nozy, 0.02, segment);
+            
+            /* closing segment */
+            segment[nsegments].x1 = nozx;
+            segment[nsegments].y1 = nozy;
+            segment[nsegments].x2 = nozx;
+            segment[nsegments].y2 = -nozy;
+            nsegments++;
+        
+            cycle = 0;
+            for (i=0; i<nsegments; i++) segment[i].group = 0;
+            
+            break;
+        }
+        case (S_TWO_CIRCLES_EXT):
+        {
+            angle = DPI/(double)NPOLY;
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                segment[i].x1 = SEGMENTS_X0 + LAMBDA*cos(((double)i)*angle);
+                segment[i].y1 = SEGMENTS_Y0 - LAMBDA*sin(((double)i)*angle);
+                segment[i].x2 = SEGMENTS_X0 + LAMBDA*cos(((double)(i+1))*angle);
+                segment[i].y2 = SEGMENTS_Y0 - LAMBDA*sin(((double)(i+1))*angle);
+                segment[i].concave = 1;
+                segment[i].group = 0;
+            }
+            for (i=NPOLY; i<2*NPOLY; i++)
+            {
+                segment[i].x1 = -SEGMENTS_X0 + TWO_CIRCLES_RADIUS_RATIO*LAMBDA*cos(((double)i)*angle);
+                segment[i].y1 = SEGMENTS_Y0 - TWO_CIRCLES_RADIUS_RATIO*LAMBDA*sin(((double)i)*angle);
+                segment[i].x2 = -SEGMENTS_X0 + TWO_CIRCLES_RADIUS_RATIO*LAMBDA*cos(((double)(i+1))*angle);
+                segment[i].y2 = SEGMENTS_Y0 - TWO_CIRCLES_RADIUS_RATIO*LAMBDA*sin(((double)(i+1))*angle);
+                segment[i].concave = 1;
+                segment[i].group = 1;
+            }
+            
+            cycle = 0;
+            nsegments = 2*NPOLY;
+            
+            break;
+        }
         default: 
         {
             printf("Function init_segment_config not defined for this pattern \n");
@@ -1077,6 +1445,11 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
         segment[i].x2 = segment[(i+1)%(nsegments)].x1;
         segment[i].y2 = segment[(i+1)%(nsegments)].y1;
     }
+    else if (SEGMENT_PATTERN != S_TWO_CIRCLES_EXT) for (i=0; i<nsegments; i++) if (segment[i].cycle)
+    {
+        segment[i].x2 = segment[(i+1)%(nsegments)].x1;
+        segment[i].y2 = segment[(i+1)%(nsegments)].y1;
+    }
     
     /* add one segment for S_POLY_ELLIPSE configuration */
     if (SEGMENT_PATTERN == S_POLY_ELLIPSE)
@@ -1091,6 +1464,10 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
     /* activate all segments */
     for (i=0; i<nsegments; i++) segment[i].active = 1;
     
+    /* inactivate some segments of leaky centrifuge */
+//     if (SEGMENT_PATTERN == S_CENTRIFUGE_LEAKY)
+//         for (i=0; i<NPOLY; i++) segment[(nsides+2)*i + nsides/2].active = 0;
+    
     /* compute parameters for slope and normal of segments */
     for (i=0; i<nsegments; i++)
     {
@@ -1108,8 +1485,20 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
     /* deal with concave corners */
     for (i=0; i<nsegments; i++) if (segment[i].concave)
         {
-            iminus = i-1;  if (iminus == 0) iminus = nsegments - 1;
-            iplus = i+1;  if (iplus == nsegments - 1) iplus = 0;
+            iminus = i-1;  
+            iplus = i+1;  
+            if (SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)
+            {
+                if (iminus == -1) iminus += nsegments/2;
+                else if (iminus == nsegments/2 - 1) iminus += nsegments/2;
+                if (iplus == nsegments/2) iplus = 0;
+                else if (iplus == nsegments) iplus = nsegments/2;
+            }
+            else
+            {
+                if (iminus == 0) iminus = nsegments - 1;
+                if (iplus == nsegments - 1) iplus = 0;
+            }
             angle = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID;
             angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID;
             if (angle2 < angle) angle2 += DPI;
@@ -1141,7 +1530,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
 int in_segment_region(double x, double y)
 /* returns 1 if (x,y) is inside region delimited by obstacle segments */
 {
-    double angle, dx, height, width, theta;
+    double angle, dx, height, width, theta, lx, ly;
     
     if (x >= BCXMAX) return(0);
     if (x <= BCXMIN) return(0);
@@ -1194,6 +1583,51 @@ int in_segment_region(double x, double y)
             if (x*x + y*y/(LAMBDA*LAMBDA) < 0.95) return(1);
             else return(0);
         }
+        case (S_CENTRIFUGE_RND):
+        {
+            if (module2(x,y) > 0.75) return(0);
+            angle = argument(x,y);
+            theta = DPI/(double)NPOLY;
+            while (angle > theta) angle -= theta;
+            while (angle < 0.0) angle += theta;
+            if (angle < 0.1) return(0);
+            if (angle > 0.9) return(0);
+            return(1);
+        }
+        case (S_CENTRIFUGE_LEAKY):
+        {
+            if (module2(x,y) > 0.75) return(0);
+            angle = argument(x,y);
+            theta = DPI/(double)NPOLY;
+            while (angle > theta) angle -= theta;
+            while (angle < 0.0) angle += theta;
+            if (angle < 0.1) return(0);
+            if (angle > 0.9) return(0);
+            return(1);
+        }
+        case (S_CIRCLE_EXT):
+        {
+            if (module2(x - SEGMENTS_X0, y - SEGMENTS_Y0) > LAMBDA + 2.0*MU) return(1);
+            else return(0);
+        }
+        case (S_TWO_CIRCLES_EXT):
+        {
+            if ((module2(x - SEGMENTS_X0, y - SEGMENTS_Y0) > LAMBDA + 2.0*MU)&&(module2(x + SEGMENTS_X0, y - SEGMENTS_Y0) > TWO_CIRCLES_RADIUS_RATIO*LAMBDA + 2.0*MU)) return(1);
+            else return(0);
+        }
+        case (S_ROCKET_NOZZLE):
+        {
+            if (x < 0.0) return(0);
+            else if (x > 1.4*LAMBDA) return(0);
+            else 
+            {
+                lx = 0.7*LAMBDA;
+                ly = 0.5*LAMBDA;
+                x -= lx;
+                if (x*x/(lx*lx) + y*y/(ly*ly) < 0.95) return(1);
+            }
+            return(0);
+        }
         default: return(1);
     }
 }
@@ -1231,17 +1665,27 @@ void rotate_segments(t_segment segment[NMAXSEGMENTS], double angle)
     
 }
  
-void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax, double deltay)
+void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax[2], double deltay[2])
 /* rotates the repelling segments by given angle */
 {
-    int i;
+    int i, group;
     
     for (i=0; i<nsegments; i++) 
     {
-        segment[i].x1 = segment[i].x01 + deltax;
-        segment[i].y1 = segment[i].y01 + deltay;
-        segment[i].x2 = segment[i].x02 + deltax;
-        segment[i].y2 = segment[i].y02 + deltay;
+        group = segment[i].group;
+        if (group == 0)
+        {
+            segment[i].x1 = segment[i].x01 + deltax[group] - SEGMENTS_X0;
+            segment[i].x2 = segment[i].x02 + deltax[group] - SEGMENTS_X0;
+        }
+        else
+        {
+            segment[i].x1 = segment[i].x01 + deltax[group] + SEGMENTS_X0;
+            segment[i].x2 = segment[i].x02 + deltax[group] + SEGMENTS_X0;
+        }
+            
+        segment[i].y1 = segment[i].y01 + deltay[group] - SEGMENTS_Y0;
+        segment[i].y2 = segment[i].y02 + deltay[group] - SEGMENTS_Y0;
         segment[i].c = segment[i].nx*segment[i].x1 + segment[i].ny*segment[i].y1;
     }
 }
@@ -2074,6 +2518,63 @@ void compute_entropy(t_particle particle[NMAXCIRCLES], double entropy[2])
     else entropy[1] = -(p2*log(p2) + (1.0-p2)*log(1.0-p2)/log2);
 }
 
+
+void draw_triangles(t_particle particle[NMAXCIRCLES])
+/* fill triangles between neighboring particles */
+{
+    int i, j, k, p, q, t0, tj, tmax, nsame = 0, same_table[9*HASHMAX];
+    double rgb[3], hue, dx, dy, x, y;
+    
+    printf("Number of nbs: ");
+    for (i=0; i<ncircles; i++) if (particle[i].active)
+    {
+        nsame = 0;
+        t0 = particle[i].type;
+        
+        /* find neighbours of same type */
+        for (j=0; j<particle[i].hash_nneighb; j++)
+        {
+            k = particle[i].hashneighbour[j];
+            if ((k!=i)&&(particle[k].type == t0)) 
+            {
+                same_table[nsame] = j;
+                nsame++;
+            }
+        }
+                
+        /* draw the triangles */
+        if (nsame > 1)
+        {
+            if (t0 <= 1) hue = HUE_TYPE0;
+            else if (t0 == 2) hue = HUE_TYPE1;
+            else if (t0 == 3) hue = HUE_TYPE2;
+            else hue = HUE_TYPE3;
+                
+            hsl_to_rgb(hue, 0.9, 0.3, rgb);
+            glColor3f(rgb[0], rgb[1], rgb[2]);
+            for (p=-1; p<2; p++) for (q=-1; q<2; q++)
+            {
+                x = particle[i].xc + (double)p*(BCXMAX - BCXMIN);
+                y = particle[i].yc + (double)q*(BCYMAX - BCYMIN);
+            
+                glBegin(GL_TRIANGLE_FAN);
+                glVertex2d(x, y);
+//                 printf("%i | ", particle[i].hash_nneighb);
+                for (k=0; k<nsame; k++) 
+                {
+                    dx = particle[i].deltax[same_table[k]];
+                    dy = particle[i].deltay[same_table[k]];
+                    if (module2(dx, dy) < particle[i].radius*NBH_DIST_FACTOR)
+                        glVertex2d(x + dx, y + dy);
+                }
+                glEnd();
+            }
+        }
+    }
+    printf("\n");
+}
+
+
 void draw_one_particle_links(t_particle particle)
 /* draw links of one particle */
 {
@@ -2224,22 +2725,42 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES],
 }
     
 
-void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
+void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta)
 {
     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;
+    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;
     char message[100];
     
     if (!TRACER_PARTICLE) blank();
     glColor3f(1.0, 1.0, 1.0);
     
+    /* show region of partial thermostat */
+    if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_INBOX))
+    {
+        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);
+        
+    }
+    
     /* draw the bonds first */
-    if (plot == P_BONDS)
+    if ((DRAW_BONDS)||(plot == P_BONDS))
     {
         glLineWidth(LINK_WIDTH);
         for (j=0; j<ncircles; j++) if (particle[j].active) draw_one_particle_links(particle[j]);
     }
                 
+    /* fill triangles between particles */
+    if (FILL_TRIANGLES) draw_triangles(particle);
+    
+    
     /* determine particle color and size */
     for (j=0; j<ncircles; j++) if (particle[j].active)
     {
@@ -2741,6 +3262,11 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
         {
             hsl_to_rgb(300.0, 0.1, 0.5, rgb);
             draw_colored_rectangle(0.0, 0.0, BCXMAX, BCYMAX, rgb);
+            break;
+        }
+        default: 
+        {
+            /* do nothing */
         }
     }
     
@@ -2823,6 +3349,7 @@ void print_parameters(double beta, double temperature, double krepel, double len
     {
         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;
         erase_area_hsl_turbo(xbox, y + 0.025, 0.37, 0.05, hue, 0.9, 0.5);
@@ -3105,8 +3632,6 @@ void print_omega(double angular_speed)
     
     if (first)
     {
-//         xleftbox = XMIN + 0.4;
-//         xlefttext = xleftbox - 0.55;
         xrightbox = XMAX - 0.39;
         xrighttext = xrightbox - 0.55;
         first = 0;
@@ -3115,11 +3640,55 @@ void print_omega(double angular_speed)
 
     erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0);
     glColor3f(1.0, 1.0, 1.0);
-    sprintf(message, "Angular speed = %.4f", DPI*angular_speed*25.0/(double)(PERIOD_ROTATE_BOUNDARY));
+    sprintf(message, "Angular speed = %.4f", angular_speed);
+//     sprintf(message, "Angular speed = %.4f", DPI*angular_speed*25.0/(double)(PERIOD_ROTATE_BOUNDARY));
     write_text(xrighttext + 0.1, y, message);
 
 }
 
+void print_segments_speeds(double vx[2], double vy[2])
+{
+    char message[100];
+    double rgb[3], y = YMAX - 0.1;
+    static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05;
+    static int first = 1;
+    
+    if (first)
+    {
+        xleftbox = XMIN + 0.4;
+        xlefttext = xleftbox - 0.3;
+        xrightbox = XMAX - 0.39;
+        xrighttext = xrightbox - 0.3;
+        first = 0;
+    }
+    
+    erase_area_hsl(xrightbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0);
+    glColor3f(1.0, 1.0, 1.0);
+    sprintf(message, "Vx = %.2f", vx[0]);
+    write_text(xrighttext + 0.1, y, message);
+
+    y -= 0.1;
+    erase_area_hsl(xrightbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0);
+    glColor3f(1.0, 1.0, 1.0);
+    sprintf(message, "Vy = %.2f", vy[0]);
+    write_text(xrighttext + 0.1, y, message);
+    
+    if (TWO_OBSTACLES)
+    {
+        y = YMAX - 0.1;
+        erase_area_hsl(xleftbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "Vx = %.2f", vx[1]);
+        write_text(xlefttext + 0.1, y, message);
+
+        y -= 0.1;
+        erase_area_hsl(xleftbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "Vy = %.2f", vy[1]);
+        write_text(xlefttext + 0.1, y, message);
+    }
+}
+
 void print_particles_speeds(t_particle particle[NMAXCIRCLES])
 {
     char message[100];
@@ -3606,6 +4175,34 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
             }
             return(fperp);
         }
+        case (BC_ABSORBING):
+        {
+            /* add harmonic force outside screen */
+            padding = 0.1;
+            x = particle[j].xc;
+            y = particle[j].yc;
+            
+            if ((x > XMAX + padding)||(x < XMIN - padding)||(y > YMAX + padding)||(y < YMIN - padding)) 
+            {
+                particle[j].active = 0;
+                particle[j].vx = 0.0;
+                particle[j].vy = 0.0;
+                particle[j].xc = XMAX + 2.0*padding;
+                particle[j].yc = YMAX + 2.0*padding;
+            }
+            
+//             if ((x > XMAX)&&(x < XMAX + 0.5*padding)) particle[j].fx += KSPRING_BOUNDARY*(XMAX + 0.5*padding - x);
+//             else if (x < XMIN) particle[j].fx -= KSPRING_BOUNDARY;
+//             if (y > YMAX) particle[j].fy += KSPRING_BOUNDARY;
+//             else if (y < YMIN) particle[j].fy -= KSPRING_BOUNDARY;
+            
+//             if (x > XMAX + padding) particle[j].fx -= KSPRING_BOUNDARY*(x - XMAX - padding);
+//             else if (x < XMIN - padding) particle[j].fx += KSPRING_BOUNDARY*(XMIN - padding - x);
+//             if (y > YMAX + padding) particle[j].fy -= KSPRING_BOUNDARY*(y - YMAX - padding);
+//             else if (y < YMIN - padding) particle[j].fy += KSPRING_BOUNDARY*(YMIN - padding - y);
+            
+            return(fperp);
+        }
     }
 }
 
@@ -3907,13 +4504,15 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub
 //             px[ncircles - 1] = particle[ncircles - 1].vx;
 //             py[ncircles - 1] = particle[ncircles - 1].vy;
 //     add_particle(MU*(2.0*rand()/RAND_MAX - 1.0), YMAX + 2.0*MU, 0.0, 0.0, PARTICLE_MASS, 0, particle);
-            
+     
+//     add_particle(XMIN - 0.5*MU, 0.0, 50.0 + 5.0*(double)i, 0.0, 2.0*PARTICLE_MASS, 0, particle);
+    
     printf("Adding a particle\n\n");
     
-    add_particle(XMIN - 0.5*MU, 0.0, 50.0 + 5.0*(double)i, 0.0, 2.0*PARTICLE_MASS, 0, particle);
-    i++;
+    add_particle(BCXMIN + 0.1, 0.5*(BCYMIN + BCYMAX), 200.0, 0.0, PARTICLE_MASS, 0, particle);
+//     i++;
     
-    particle[ncircles - 1].radius = 0.5*MU;
+    particle[ncircles - 1].radius = MU;
     particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST;
     particle[ncircles - 1].thermostat = 0;
     px[ncircles - 1] = particle[ncircles - 1].vx;
@@ -3959,6 +4558,7 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin)
 /* only couple particles with x > xmin to thermostat */
 {
     int condition, i, nthermo = 0;
+    double x, y;
     
     for (i=0; i<ncircles; i++) 
     {
@@ -3970,7 +4570,17 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin)
             }
             case (TH_INSEGMENT):
             {
-                condition = (in_segment_region(particle[i].xc - xsegments, particle[i].yc - ysegments));
+                condition = (in_segment_region(particle[i].xc - xsegments[0], particle[i].yc - ysegments[0]));
+//                 condition = (in_segment_region(particle[i].xc - xsegments[0], particle[i].yc - ysegments[0]));
+//                 if (TWO_OBSTACLES) 
+//                     condition = condition||(in_segment_region(particle[i].xc - xsegments[1], particle[i].yc - ysegments[1]));
+                break;
+            }
+            case (TH_INBOX):
+            {
+                x = particle[i].xc;
+                y = particle[i].yc;
+                condition = ((y < YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN))&&(vabs(x) < PARTIAL_THERMO_WIDTH*XMAX));
                 break;
             }
             default: condition = 1;
diff --git a/sub_part_billiard.c b/sub_part_billiard.c
index 725e4de..15b66ce 100644
--- a/sub_part_billiard.c
+++ b/sub_part_billiard.c
@@ -276,6 +276,14 @@ void erase_rectangle_outside(double h, double s, double l)
     glEnd();
 }
 
+void draw_line(double x1, double y1, double x2, double y2)
+{
+    glBegin(GL_LINE_STRIP);
+    glVertex2d(x1, y1);
+    glVertex2d(x2, y2);
+    glEnd();    
+}
+
 void draw_circle(double x, double y, double r, int nseg)
 {
     int i;
@@ -1323,21 +1331,32 @@ void draw_billiard()      /* draws the billiard boundary */
             }
             if (FOCI)
             {
-                if (POLYLINE_PATTERN == P_TOKARSKY)
-                {
-                    glLineWidth(2);
-                    rgb[0] = 1.0;   rgb[1] = 0.0;   rgb[2] = 0.0;
-                    draw_colored_circle(-0.95, 0.0, r, NSEG, rgb);
-                    rgb[0] = 0.0;   rgb[1] = 0.8;   rgb[2] = 0.2;
-                    draw_colored_circle(0.95, 0.0, r, NSEG, rgb);
-                }
-                else if (POLYLINE_PATTERN == P_TOKA_PRIME)
-                {
-                    glLineWidth(2);
-                    rgb[0] = 1.0;   rgb[1] = 0.0;   rgb[2] = 0.0;
-                    draw_colored_circle(-polyline[84].x1, polyline[84].y1, r, NSEG, rgb);
-                    rgb[0] = 0.0;   rgb[1] = 0.8;   rgb[2] = 0.2;
-                    draw_colored_circle(polyline[84].x1, polyline[84].y1, r, NSEG, rgb);
+                switch (POLYLINE_PATTERN) {
+                    case (P_TOKARSKY):
+                    {
+                        glLineWidth(2);
+                        rgb[0] = 1.0;   rgb[1] = 0.0;   rgb[2] = 0.0;
+                        draw_colored_circle(-0.95, 0.0, r, NSEG, rgb);
+                        rgb[0] = 0.0;   rgb[1] = 0.8;   rgb[2] = 0.2;
+                        draw_colored_circle(0.95, 0.0, r, NSEG, rgb);
+                        break;
+                    }
+                    case (P_TOKA_PRIME):
+                    {
+                        glLineWidth(2);
+                        rgb[0] = 1.0;   rgb[1] = 0.0;   rgb[2] = 0.0;
+                        draw_colored_circle(-polyline[84].x1, polyline[84].y1, r, NSEG, rgb);
+                        rgb[0] = 0.0;   rgb[1] = 0.8;   rgb[2] = 0.2;
+                        draw_colored_circle(polyline[84].x1, polyline[84].y1, r, NSEG, rgb);
+                        break;
+                    }
+                    case (P_TOKA_NONSELF):
+                    {
+                        glLineWidth(2);
+                        rgb[0] = 0.0;   rgb[1] = 0.8;   rgb[2] = 0.2;
+                        draw_colored_circle(0.0, 0.0, r, NSEG, rgb);
+                        break;
+                    }
                 }
             }
             if (ABSORBING_CIRCLES)
@@ -6083,8 +6102,53 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES])
                 circles[i].radius = MU;
                 circles[i].active = 1;
             }
-    
+            break;
+        }
+        case (P_TOKA_NONSELF):
+        {
+            nsides = 12;
+            ncircles = 12;
             
+            polyline[0].x1 = 0.0;    polyline[0].y1 = 2.0;       polyline[0].angle = 3.0*PID;
+            polyline[1].x1 = 0.0;    polyline[1].y1 = 1.0;       polyline[1].angle = 0.0;
+            polyline[2].x1 = 1.0;    polyline[2].y1 = 1.0;       polyline[2].angle = 3.5*PID;
+            polyline[3].x1 = 2.0;    polyline[3].y1 = 0.0;       polyline[3].angle = PI;
+            polyline[4].x1 = 1.0;    polyline[4].y1 = 0.0;       polyline[4].angle = 3.0*PID;
+            polyline[5].x1 = 1.0;    polyline[5].y1 = -1.0;      polyline[5].angle = 2.5*PID;
+            polyline[6].x1 = 0.0;    polyline[6].y1 = -2.0;      polyline[6].angle = PID;
+            polyline[7].x1 = 0.0;    polyline[7].y1 = -1.0;      polyline[7].angle = PI;
+            polyline[8].x1 = -1.0;   polyline[8].y1 = -1.0;      polyline[8].angle = 1.5*PID;
+            polyline[9].x1 = -2.0;   polyline[9].y1 = 0.0;       polyline[9].angle = 0.0;
+            polyline[10].x1 = -1.0;  polyline[10].y1 = 0.0;      polyline[10].angle = PID;
+            polyline[11].x1 = -1.0;  polyline[11].y1 = 1.0;     polyline[11].angle = 0.5*PID;
+            
+            ratio = (YMAX - YMIN)/4.5;
+            for (i=0; i<nsides; i++)
+            {
+                polyline[i].x1 = ratio*(polyline[i].x1);
+                polyline[i].y1 = ratio*(polyline[i].y1);
+            }
+                
+            for (i=0; i<nsides; i++) if (i < nsides-1) 
+            {   
+                polyline[i].x2 = polyline[i+1].x1;
+                polyline[i].y2 = polyline[i+1].y1;
+            }
+            
+            polyline[nsides-1].x2 = polyline[0].x1;
+            polyline[nsides-1].y2 = polyline[0].y1;
+            
+            for (i=0; i<nsides; i++) 
+                polyline[i].length = module2(polyline[i].x2 - polyline[i].x1, polyline[i].y2 - polyline[i].y1);
+            
+            for (i=0; i<ncircles; i++)
+            {
+                circles[i].xc = polyline[i].x1;
+                circles[i].yc = polyline[i].y1;
+                circles[i].radius = MU;
+                circles[i].active = 1;
+            }
+            break;
         }
     }
 }
diff --git a/sub_rde.c b/sub_rde.c
index 98c029a..cd8f1e2 100644
--- a/sub_rde.c
+++ b/sub_rde.c
@@ -399,13 +399,29 @@ void compute_field_argument(double *phi[NFIELDS], t_rde rde[NX*NY])
         for (j=0; j<NY; j++)
         {
             arg = argument(phi[0][i*NY+j], phi[1][i*NY+j]) + COLOR_PHASE_SHIFT*PI;
-            if (arg < 0.0) arg += DPI;
-            if (arg >= DPI) arg -= DPI;
+            while (arg < 0.0) arg += DPI;
+            while (arg >= DPI) arg -= DPI;
             rde[i*NY+j].field_arg = arg;
         }
 }
 
-
+void compute_probabilities(t_rde rde[NX*NY], short int xy_in[NX*NY], double probas[2])
+/* compute probabilities for Ehrenfest urns */
+{
+    int i, j;
+    double pleft = 0.0, pright = 0.0, sum;
+    
+    #pragma omp parallel for private(j)
+    for (j=0; j<NY; j++)
+    {
+        for (i=0; i<NX/2; i++) if (xy_in[i*NY+j]) pleft += rde[i*NY+j].field_norm;
+        for (i=NX/2; i<NX; i++) if (xy_in[i*NY+j]) pright += rde[i*NY+j].field_norm;
+    }
+        
+    sum = pleft + pright;
+    probas[0] = pleft/sum;
+    probas[1] = pright/sum;
+}
 
 // void compute_field_norm(double phi_x[NX*NY], double phi_y[NX*NY], double phi_norm[NX*NY], double factor)
 // /* compute the norm of (phi_x, phi_y) */
@@ -487,7 +503,7 @@ void compute_laplacian(double phi_in[NX*NY], double phi_out[NX*NY], short int xy
     }
 }
 
-void compute_light_angle_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], int movie)
+void compute_light_angle_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], int movie)
 /* computes cosine of angle between normal vector and vector light */
 {
     int i, j;
@@ -510,6 +526,14 @@ void compute_light_angle_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], int movie
             {
                 gradx = (*rde[(i+1)*NY+j].p_zfield[movie] - *rde[(i-1)*NY+j].p_zfield[movie])/dx;
                 grady = (*rde[i*NY+j+1].p_zfield[movie] - *rde[i*NY+j-1].p_zfield[movie])/dy;
+                
+                /* case where the potential is added to the z coordinate */
+                if ((ADD_POTENTIAL)&&(ADD_POTENTIAL_TO_Z))
+                {
+                    gradx += ADD_POT_CONSTANT*(potential[(i+1)*NY+j] - potential[(i-1)*NY+j])/dx;
+                    grady += ADD_POT_CONSTANT*(potential[i*NY+j+1] - potential[i*NY+j-1])/dx;
+                }
+                
                 norm = sqrt(1.0 + gradx*gradx + grady*grady);
                 pscal = -gradx*light[0] - grady*light[1] + 1.0;
                 
@@ -603,7 +627,8 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
         }
         case (Z_NORM_GRADIENT): 
         {
-            color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
+//             color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
+            color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
             break;
         }
         case (Z_ANGLE_GRADIENT): 
@@ -631,6 +656,11 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
             hsl_to_rgb_palette(180.0*(1.0 - color_amplitude(value, 1.0, 0)), 0.9, 0.5, rgb, palette);
             break;
         }
+        case (Z_VORTICITY_ABS): 
+        {
+            hsl_to_rgb_palette(360.0*(1.0 - vabs(color_amplitude(value, 1.0, 0))), 0.9, 0.5, rgb, palette);
+            break;
+        }
         case (Z_MAXTYPE_RPSLZ): 
         {
             hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette);
@@ -653,25 +683,47 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
         }
         case (Z_ARGUMENT): 
         {
-            color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb);
+//             color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb);
+            hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette);
+            break;
+        }
+        case (Z_REALPART): 
+        {
+            color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*value, 1.0, 0, rgb);
             break;
         }
     }
 }
 
+double adjust_field(double z, double pot)
+/* add potential in case of option ADD_POTENTIAL_TO_Z */
+{
+    if ((ADD_POTENTIAL)&&(ADD_POTENTIAL_TO_Z)) return (z + ADD_POT_CONSTANT*pot);
+    else return(z);
+}  
 
-double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double palette, int cplot, 
+
+double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double potential[NX*NY], double palette, int cplot, 
                                  double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3],
+                                 double *z_sw, double *z_se, double *z_nw, double *z_ne, 
                                  int fade, double fade_value, int movie)
 {
     int k;
     double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid;
-    double *z_sw, *z_se, *z_nw, *z_ne, lum;
+    double lum;
     
-    z_sw = rde[i*NY+j].p_zfield[movie];
-    z_se = rde[(i+1)*NY+j].p_zfield[movie];
-    z_nw = rde[i*NY+j+1].p_zfield[movie];
-    z_ne = rde[(i+1)*NY+j+1].p_zfield[movie];
+    *z_sw = *rde[i*NY+j].p_zfield[movie];
+    *z_se = *rde[(i+1)*NY+j].p_zfield[movie];
+    *z_nw = *rde[i*NY+j+1].p_zfield[movie];
+    *z_ne = *rde[(i+1)*NY+j+1].p_zfield[movie];
+    
+    if ((ADD_POTENTIAL)&&(ADD_POTENTIAL_TO_Z))
+    {
+        *z_sw += ADD_POT_CONSTANT*potential[i*NY+j];
+        *z_se += ADD_POT_CONSTANT*potential[(i+1)*NY+j];
+        *z_nw += ADD_POT_CONSTANT*potential[i*NY+j+1];
+        *z_ne += ADD_POT_CONSTANT*potential[(i+1)*NY+j+1];
+    }
                             
     z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne);
     
@@ -707,9 +759,10 @@ double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double pa
     compute_field_color_rde(cn, cplot, palette, rgb_n);
     compute_field_color_rde(cs, cplot, palette, rgb_s);
     
+//     if ((cplot == Z_ARGUMENT)||(cplot == Z_REALPART))
     if (cplot == Z_ARGUMENT)
     {
-        lum = tanh(SLOPE_SCHROD_LUM*rde[i*NY+j].field_norm);
+        lum = tanh(SLOPE_SCHROD_LUM*rde[i*NY+j].field_norm) + MIN_SCHROD_LUM;
         for (k=0; k<3; k++) 
         {
             rgb_e[k] *= lum;
@@ -767,7 +820,7 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot,
                 compute_gradient_polar(rde, 0.03);
             }
     
-            if ((zplot == Z_VORTICITY)||(cplot == Z_VORTICITY))
+            if ((zplot == Z_VORTICITY)||(cplot == Z_VORTICITY)||(zplot == Z_VORTICITY_ABS)||(cplot == Z_VORTICITY_ABS))
             {
                 compute_gradient_theta(rde);
                 compute_curl(rde);
@@ -851,6 +904,12 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_
             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &rde[i*NY+j].curl;
             break;
         }
+        case (Z_VORTICITY_ABS):
+        {
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &rde[i*NY+j].curl;
+            break;
+        }
         case (Z_MAXTYPE_RPSLZ):
         {
             #pragma omp parallel for private(i,j)
@@ -881,6 +940,12 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_
             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &rde[i*NY+j].field_arg;
             break;
         }
+        case (Z_REALPART):
+        {
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &phi[0][i*NY+j];
+            break;
+        }
     }
 }
 
@@ -938,6 +1003,12 @@ void init_cfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int cplot, t_
             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].curl;
             break;
         }
+        case (Z_VORTICITY_ABS):
+        {
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].curl;
+            break;
+        }
         case (Z_MAXTYPE_RPSLZ):
         {
             #pragma omp parallel for private(i,j)
@@ -968,6 +1039,12 @@ void init_cfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int cplot, t_
             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].field_arg;
             break;
         }
+        case (Z_REALPART):
+        {
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &phi[0][i*NY+j];
+            break;
+        }
     }
 }
 
@@ -981,6 +1058,7 @@ void compute_cfield_rde(short int xy_in[NX*NY], int cplot, int palette, t_rde rd
     for (i=0; i<NX; i++) for (j=0; j<NY; j++)
     {
         compute_field_color_rde(*rde[i*NY+j].p_cfield[movie], cplot, palette, rde[i*NY+j].rgb);
+//         if ((cplot == Z_ARGUMENT)||(cplot == Z_REALPART))
         if (cplot == Z_ARGUMENT)
         {
             lum = tanh(SLOPE_SCHROD_LUM*rde[i*NY+j].field_norm);
@@ -1019,86 +1097,130 @@ void draw_wave_2d_rde(short int xy_in[NX*NY], t_rde rde[NX*NY])
             }
         }
     glEnd ();
-    if (DRAW_BILLIARD) draw_billiard();
+    if (DRAW_BILLIARD) draw_billiard(0, 1.0);
+}
+
+
+void draw_wave_3d_ij_rde(int i, int j, int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], 
+                         double potential[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value)
+{
+    int k, l, draw = 1;
+    double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; 
+    double z, z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0;
+    double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2];
+    double energy;
+    
+
+    if (NON_DIRICHLET_BC) 
+        draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
+    else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
+            
+    if (draw)
+    {
+        if (AMPLITUDE_HIGH_RES > 0)
+        {
+            z_mid = compute_interpolated_colors_rde(i, j, rde, potential, palette, cplot, 
+                                                        rgb_e, rgb_w, rgb_n, rgb_s, &z_sw, &z_se, &z_nw, &z_ne, 
+                                                        fade, fade_value, movie);
+            ij_to_xy(i, j, xy_sw);
+            ij_to_xy(i+1, j, xy_se);
+            ij_to_xy(i, j+1, xy_nw);
+            ij_to_xy(i+1, j+1, xy_ne);
+                    
+            for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]);
+                       
+            if (AMPLITUDE_HIGH_RES == 1)
+            {                        
+                glBegin(GL_TRIANGLE_FAN);
+                glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_nw, z_nw);
+                draw_vertex_xyz(xy_sw, z_sw);
+                    
+                glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
+                draw_vertex_xyz(xy_se, z_se);
+                    
+                glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
+                draw_vertex_xyz(xy_ne, z_ne);
+                    
+                glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
+                draw_vertex_xyz(xy_nw, z_nw);
+                glEnd ();
+            }
+        }
+        else
+        {
+            glColor3f(rde[i*NY+j].rgb[0], rde[i*NY+j].rgb[1], rde[i*NY+j].rgb[2]);
+            
+            glBegin(GL_TRIANGLE_FAN);
+            ij_to_xy(i, j, xy);
+            draw_vertex_xyz(xy, adjust_field(*rde[i*NY+j].p_zfield[movie], potential[i*NY+j]));
+            ij_to_xy(i+1, j, xy);
+            draw_vertex_xyz(xy, adjust_field(*rde[(i+1)*NY+j].p_zfield[movie], potential[(i+1)*NY+j]));
+            ij_to_xy(i+1, j+1, xy);
+            draw_vertex_xyz(xy, adjust_field(*rde[(i+1)*NY+j+1].p_zfield[movie], potential[(i+1)*NY+j+1]));
+            ij_to_xy(i, j+1, xy);
+            draw_vertex_xyz(xy, adjust_field(*rde[i*NY+j+1].p_zfield[movie], potential[i*NY+j+1]));
+            glEnd ();
+        }
+    }
 }
 
 
 
-void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], 
+void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], 
                   int zplot, int cplot, int palette, int fade, double fade_value)
 {
-    int i, j, k, l, draw = 1;
-    double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; 
-    double z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0;
-    double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2];
-    double energy;
+    int i, j;
+    double observer_angle;
     
     blank();
     if (DRAW_BILLIARD) draw_billiard_3d(fade, fade_value);
-            
-    for (i=0; i<NX-2; i++)
-        for (j=0; j<NY-2; j++)
-        {
-            if (NON_DIRICHLET_BC) 
-                draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
-            else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
-            
-            if (draw)
-            {
-                if (AMPLITUDE_HIGH_RES > 0)
-                {
-                    z_mid = compute_interpolated_colors_rde(i, j, rde, palette, cplot, 
-                                                             rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie);
-                    
-                    ij_to_xy(i, j, xy_sw);
-                    ij_to_xy(i+1, j, xy_se);
-                    ij_to_xy(i, j+1, xy_nw);
-                    ij_to_xy(i+1, j+1, xy_ne);
-                    
-                    for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]);
-                       
-                    if (AMPLITUDE_HIGH_RES == 1)
-                    {
-                        glBegin(GL_TRIANGLE_FAN);
-                        glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_nw, *rde[i*NY+j+1].p_zfield[movie]);
-                        draw_vertex_xyz(xy_sw, *rde[i*NY+j].p_zfield[movie]);
-                    
-                        glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
-                        draw_vertex_xyz(xy_se, *rde[(i+1)*NY+j].p_zfield[movie]);
-                    
-                        glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
-                        draw_vertex_xyz(xy_ne, *rde[(i+1)*NY+j+1].p_zfield[movie]);
-                    
-                        glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
-                        draw_vertex_xyz(xy_nw, *rde[i*NY+j+1].p_zfield[movie]);
-                        glEnd ();
-                    }
-                }
-                else
-                {
-                    glColor3f(rde[i*NY+j].rgb[0], rde[i*NY+j].rgb[1], rde[i*NY+j].rgb[2]);
     
-                    glBegin(GL_TRIANGLE_FAN);
-                    ij_to_xy(i, j, xy);
-                    draw_vertex_xyz(xy, *rde[i*NY+j].p_zfield[movie]);
-                    ij_to_xy(i+1, j, xy);
-                    draw_vertex_xyz(xy, *rde[(i+1)*NY+j].p_zfield[movie]);
-                    ij_to_xy(i+1, j+1, xy);
-                    draw_vertex_xyz(xy, *rde[(i+1)*NY+j+1].p_zfield[movie]);
-                    ij_to_xy(i, j+1, xy);
-                    draw_vertex_xyz(xy, *rde[i*NY+j+1].p_zfield[movie]);
-                    glEnd ();
-                }
-            }
-        }
+    if (!ROTATE_VIEW)
+    {
+        for (i=0; i<NX-2; i++)
+            for (j=0; j<NY-2; j++)
+                draw_wave_3d_ij_rde(i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
+    }
+    else    /* draw facets in an order depending on the position of the observer */
+    {
+        observer_angle = argument(observer[0], observer[1]);
+        observer_angle -= 0.5*PID;
+        if (observer_angle < 0.0) observer_angle += DPI;
+        printf("Observer_angle = %.3lg\n", observer_angle*360.0/DPI); 
         
+        if ((observer_angle > 0.0)&&(observer_angle < PID))
+        {
+            for (j=0; j<NY-2; j++)
+                for (i=0; i<NX-2; i++)
+                    draw_wave_3d_ij_rde(i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
+        }
+        else if (observer_angle < PI)
+        {
+            for (i=NX-3; i>0; i--)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_3d_ij_rde(i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
+        }
+        else if (observer_angle < 1.5*PI)
+        {
+             for (j=NY-3; j>0; j--)
+                for (i=0; i<NX-2; i++)
+                    draw_wave_3d_ij_rde(i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);   
+        }
+        else
+        {
+            for (i=0; i<NX-2; i++)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_3d_ij_rde(i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
+        }
+    }
+            
     if (DRAW_BILLIARD_FRONT) draw_billiard_3d_front(fade, fade_value);
 }
 
 
-void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], 
+void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], 
                    int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
 {
     int i, j, k, l, draw = 1;
@@ -1112,11 +1234,11 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd
 //         printf("Computing fields\n");
         compute_rde_fields(phi, xy_in, zplot, cplot, rde);      
 //         printf("Computed fields\n");
-        if ((PLOT_3D)&&(SHADE_3D)) compute_light_angle_rde(xy_in, rde, movie);       
+        if ((PLOT_3D)&&(SHADE_3D)) compute_light_angle_rde(xy_in, rde, potential, movie);       
     }
     compute_cfield_rde(xy_in, cplot, palette, rde, fade, fade_value, movie);
     
-    if (PLOT_3D) draw_wave_3d_rde(movie, phi, xy_in, rde, zplot, cplot, palette, fade, fade_value);
+    if (PLOT_3D) draw_wave_3d_rde(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
     else draw_wave_2d_rde(xy_in, rde);
 }
 
diff --git a/sub_wave.c b/sub_wave.c
index 8e428ac..6aaa59e 100644
--- a/sub_wave.c
+++ b/sub_wave.c
@@ -1461,6 +1461,56 @@ int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY])
     return(2*n);
 }
 
+int compute_noisepanel_rect_coordinates(t_vertex polyline[NMAXPOLY])
+/* compute positions of vertices of noise panel */
+{
+    int i, n, even;
+    double ymax, dy, x, y, x1, pos[2];
+    
+    /* find the leftmost point */
+    x = -NPWIDTH;
+    n = 0;
+    while (x > XMIN) 
+    {
+        x -= LAMBDA;
+        n++;
+    }
+    if (n%2 == 0) even = 1;
+    else even = 0;
+    
+    i = 0;
+    while (x <= 0.0)
+    {
+        if (even) y = YMIN + 0.1;
+        else y = YMIN + 0.1 + MU;
+        
+        x1 = x;
+        if (x1 > XMAX) x1 = XMAX;
+        else if (x1 < XMIN) x1 = XMIN;
+        
+        polyline[i].x = x1;
+        polyline[i].y = y;
+        
+        xy_to_pos(x1, y, pos);        
+        polyline[i].posi = pos[0];
+        polyline[i].posj = pos[1];
+        
+        x += LAMBDA;
+        even = 1 - even;
+        i++;
+    }
+    n = i;
+    for (i=0; i<n; i++)
+    {
+        polyline[n+i].x = polyline[n-i-1].x;
+        polyline[n+i].y = -polyline[n-i-1].y;
+        polyline[n+i].posi = polyline[n-i-1].posi;
+        polyline[n+i].posj = NY - polyline[n-i-1].posj;
+    }
+        
+    return(2*n);
+}
+
 int compute_qrd_coordinates(t_vertex polyline[NMAXPOLY])
 /* compute positions of quadratic noise diffuser */
 {
@@ -1549,6 +1599,10 @@ int init_polyline(int depth, t_vertex polyline[NMAXPOLY])
         {
             return(compute_noisepanel_coordinates(polyline));
         }
+        case (D_NOISEPANEL_RECT):
+        {
+            return(compute_noisepanel_rect_coordinates(polyline));
+        }
         case (D_QRD):
         {
             return(compute_qrd_coordinates(polyline));
@@ -1611,7 +1665,7 @@ int xy_in_billiard(double x, double y)
 /* returns 1 if (x,y) represents a point in the billiard */
 // double x, y;
 {
-    double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha;
+    double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha, s, a;
     int i, j, k, k1, k2, condition = 0, m;
     static int first = 1, nsides;
 
@@ -2047,6 +2101,22 @@ int xy_in_billiard(double x, double y)
             else y1 = 0.1 + 2.0*MU - MU*x1/LAMBDA;
             return((y > YMIN + y1)&&(y < YMAX - y1));
         }
+        case (D_NOISEPANEL_RECT):
+        {
+            x1 = -x;
+            if (x1 > NPWIDTH)
+            {
+                while (x1 > 2.0*LAMBDA) x1 -= 2.0*LAMBDA;
+                if (x1 <= LAMBDA) y1 = 0.1 + MU*x1/LAMBDA;
+                else y1 = 0.1 + 2.0*MU - MU*x1/LAMBDA;
+                return((y > YMIN + y1)&&(y < YMAX - y1)&&(x > XMIN + 0.1));
+            }
+            else if (x > NPWIDTH)
+            {
+                return((vabs(y) < YMAX - 0.1)&&(x < XMAX - 0.1));
+            }
+            else return(0);
+        }
         case (D_QRD):
         {
             x1 = vabs(x)/LAMBDA;
@@ -2055,6 +2125,25 @@ int xy_in_billiard(double x, double y)
             y1 = (MU/13.0)*(14.0 - (double)k1);
             return ((y > YMIN + y1)&&(y < YMAX - y1));
         }
+        case (D_QRD_ASYM):
+        {
+            if (y > 0.0)
+            {   
+                x1 = vabs(x)/LAMBDA;
+                k = (int)(x1 + 0.5);
+                k1 = (k*k) % 13;
+                y1 = (MU/13.0)*(14.0 - (double)k1);
+                return (y < YMAX - y1);
+            }
+            else
+            {   
+                x1 = vabs(x + 1.0)/LAMBDA;
+                k = (int)(x1 + 0.5);
+                k1 = (k*k) % 17;
+                y1 = (MU/17.0)*(18.0 - (double)k1);
+                return (y > YMIN + y1);
+            }
+        }
         case (D_CIRCLE_SEGMENT):
         {
             if (vabs(y) > 0.9*vabs(LAMBDA)) return(1);
@@ -2076,6 +2165,14 @@ int xy_in_billiard(double x, double y)
                 else return(1);                
             }
         }
+        case (D_GROOVE):
+        {
+            s = 0.85*LAMBDA;
+            a = 0.5*LAMBDA;
+            x1 = x - XMIN - (double)((int)((x - XMIN)/LAMBDA))*LAMBDA;
+            if (x1 < a) return (y > YMIN + LAMBDA);
+            else return (y > YMIN + LAMBDA + s);
+        }
         case (D_MENGER):       
         {
             x1 = 0.5*(x+1.0);
@@ -2255,14 +2352,22 @@ void tvertex_lineto(t_vertex z)
 }
 
 
-void draw_billiard()      /* draws the billiard boundary */
+void draw_billiard(int fade, double fade_value)      /* draws the billiard boundary */
 {
     double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax;
     int i, j, k, k1, k2, mr2;
     static int first = 1, nsides;
 
-    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(BOUNDARY_WIDTH);
 
     glEnable(GL_LINE_SMOOTH);
@@ -2298,7 +2403,8 @@ void draw_billiard()      /* draws the billiard boundary */
             /* draw foci */
             if (FOCI)
             {
-                glColor3f(0.3, 0.3, 0.3);
+                if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value);
+                else glColor3f(0.3, 0.3, 0.3);
                 x0 = sqrt(LAMBDA*LAMBDA-1.0);
 
                 glLineWidth(2);
@@ -3041,6 +3147,14 @@ void draw_billiard()      /* draws the billiard boundary */
             glEnd();
             break;
         }
+        case (D_NOISEPANEL_RECT):
+        {
+            glLineWidth(BOUNDARY_WIDTH);
+            glBegin(GL_LINE_STRIP);
+            for (i=0; i<npolyline; i++) tvertex_lineto(polyline[i]);
+            glEnd();
+            break;
+        }
         case (D_QRD):
         {
             glLineWidth(BOUNDARY_WIDTH);
@@ -3325,7 +3439,11 @@ void draw_billiard()      /* draws the billiard boundary */
             glEnd();
             break;
         }
-       default:
+        case (D_NOTHING):
+        {
+            break;
+        }   
+        default:
         {
             printf("Function draw_billiard not defined for this billiard \n");
         }
@@ -3544,4 +3662,111 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p
     draw_rectangle(x1, y1, x2, y2);
 }
 
+void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette, int fade, double fade_value)
+{
+    int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax;
+    double y, dy, dy_e, rgb[3], value, lum, amp;
+    
+    xy_to_ij(x1, y1, ij_botleft);
+    xy_to_ij(x2, y2, ij_topright);
+    
+    rgb[0] = 0.0;   rgb[1] = 0.0;   rgb[2] = 0.0;
+    erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb);
+
+    if (ROTATE_COLOR_SCHEME)
+    {
+        jmin = ij_botleft[0];
+        jmax = ij_topright[0];
+        imin = ij_botleft[1];
+        imax = ij_topright[1];    
+    }
+    else
+    {
+        imin = ij_botleft[0];
+        imax = ij_topright[0];
+        jmin = ij_botleft[1];
+        jmax = ij_topright[1];    
+    }
+        
+        
+    glBegin(GL_QUADS);
+    dy = (max - min)/((double)(jmax - jmin));
+    dy_e = max/((double)(jmax - jmin));
+    
+    for (j = jmin; j < jmax; j++)
+    {
+        switch (plot) {
+            case (P_AMPLITUDE):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_ENERGY):
+            {
+                value = dy_e*(double)(j - jmin)*100.0/E_SCALE;
+                if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_MEAN_ENERGY):
+            {
+                value = dy_e*(double)(j - jmin)*100.0/E_SCALE;
+                if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_LOG_ENERGY):
+            {
+                value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE);
+//                 if (value <= 0.0) value = 0.0;
+                color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_LOG_MEAN_ENERGY):
+            {
+                value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE);
+//                 if (value <= 0.0) value = 0.0;
+                color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_PHASE):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+//                 lum = (color_amplitude(value, 1.0, 1))*0.5;
+//                 if (lum < 0.0) lum = 0.0;
+//                 hsl_to_rgb(value*360.0, 0.9, 0.5, rgb);
+//                 color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb);
+//                 amp = color_amplitude_linear(value, 1.0, 1);
+                amp = 0.5*color_amplitude_linear(value, 1.0, 1);
+                while (amp > 1.0) amp -= 2.0;
+                while (amp < -1.0) amp += 2.0;
+                amp_to_rgb(0.5*(1.0 + amp), rgb);
+                break;
+            }
+        }
+        if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
+        glColor3f(rgb[0], rgb[1], rgb[2]);
+        if (ROTATE_COLOR_SCHEME)
+        {
+            glVertex2i(j, imin);
+            glVertex2i(j, imax);
+            glVertex2i(j+1, imax);
+            glVertex2i(j+1, imin);            
+        }
+        else
+        {
+            glVertex2i(imin, j);
+            glVertex2i(imax, j);
+            glVertex2i(imax, j+1);
+            glVertex2i(imin, j+1);
+        }
+    }
+    glEnd ();
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    glLineWidth(BOUNDARY_WIDTH);
+    draw_rectangle(x1, y1, x2, y2);
+}
 
diff --git a/sub_wave_3d.c b/sub_wave_3d.c
index 7ce5d68..00aac62 100644
--- a/sub_wave_3d.c
+++ b/sub_wave_3d.c
@@ -21,7 +21,7 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2])
     static double n2, m2, d, sm2, sn2, v[3], h[2], plane_ratio = 0.5;
     static int first = 1;
     
-    if ((first)&&(REPRESENTATION_3D == REP_PROJ_3D))
+    if (((first)&&(REPRESENTATION_3D == REP_PROJ_3D))||(reset_view))
     {
         m2 = observer[0]*observer[0] + observer[1]*observer[1];
         n2 = m2 + observer[2]*observer[2];
@@ -34,8 +34,9 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2])
         v[1] = -observer[1]*observer[2]/(sn2*sm2);
         v[2] = m2/(sn2*sm2);
         first = 0;
-        printf("h = (%.3lg, %.3lg)\n", h[0], h[1]);
-        printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]);
+        reset_view = 0;
+//         printf("h = (%.3lg, %.3lg)\n", h[0], h[1]);
+//         printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]);
     }
     
     switch (REPRESENTATION_3D) {
@@ -170,6 +171,16 @@ void draw_billiard_3d(int fade, double fade_value)      /* draws the billiard bo
     glEnable(GL_LINE_SMOOTH);
 
     switch (B_DOMAIN) {
+        case (D_RECTANGLE):
+        {
+            glBegin(GL_LINE_LOOP);
+            draw_vertex_x_y_z(LAMBDA, -1.0, 0.0);
+            draw_vertex_x_y_z(LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, -1.0, 0.0);
+            glEnd();
+            break;
+        }
         case (D_ELLIPSE):
         {
             glBegin(GL_LINE_LOOP);
@@ -714,6 +725,10 @@ void draw_billiard_3d(int fade, double fade_value)      /* draws the billiard bo
                 if (polygons[i].active) draw_tpolygon_3d(polygons[i]);
             break;
         }
+        case (D_NOTHING):
+        {
+            break;
+        }
         default:
         {
             break;
@@ -781,6 +796,15 @@ void draw_billiard_3d_front(int fade, double fade_value)
     glEnable(GL_LINE_SMOOTH);
 
     switch (B_DOMAIN) {
+        case (D_RECTANGLE):
+        {
+            glBegin(GL_LINE_STRIP);
+            draw_vertex_x_y_z(LAMBDA, -1.0, 0.0);
+            draw_vertex_x_y_z(LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, 1.0, 0.0);
+            glEnd();
+            break;
+        }
         case (D_YOUNG):
         {
             glBegin(GL_LINE_STRIP);
@@ -1150,6 +1174,43 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
                 }
                 break;
             }
+            case (IOR_MANDELBROT_LIN):
+            {
+                #pragma omp parallel for private(i,j)
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        u = 0.0;
+                        v = 0.0;
+                        k = 0;
+                        while ((k<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+                        {
+                            u1 = u*u - v*v + x;
+                            v = 2.0*u*v + y;
+                            u = u1;
+                            k++;
+                        }
+                        if (k >= MANDELLEVEL)
+                        {
+                            tc[i*NY+j] = COURANT;
+                            tcc[i*NY+j] = courant2;
+                            tgamma[i*NY+j] = GAMMA;
+                        }
+                        else 
+                        {
+                            speed = (double)k/(double)MANDELLEVEL;
+                            if (speed < 1.0e-10) speed = 1.0e-10;
+                            else if (speed > 10.0) speed = 10.0;
+                            tcc[i*NY+j] = courant2*speed;
+                            tc[i*NY+j] = COURANT*sqrt(speed);
+                            tgamma[i*NY+j] = GAMMA;
+                        }
+                    }
+                }
+                break;
+            }
             default:
             {
                 for (i=0; i<NX; i++){
@@ -1324,15 +1385,124 @@ void compute_cfield(short int xy_in[NX*NY], int cplot, int palette, t_wave wave[
 }
 
 
-void draw_wave_3d(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], 
-                  int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
+void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], 
+                  int zplot, int cplot, int palette, int fade, double fade_value)
+/* draw wave at simulation grid point (i,j) */
 {
-    int i, j, k, l, draw = 1;
+    int k, l, draw = 1;
     double xy[2], xy2[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; 
     double z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0;
     double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2];
     double energy, zfloor = -10.0;
     
+    if (NON_DIRICHLET_BC) 
+        draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
+    else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
+            
+    if (FLOOR_ZCOORD) 
+        draw = (draw)&&(*wave[i*NY+j].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j].p_zfield[movie] > zfloor)&&(*wave[i*NY+j+1].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j+1].p_zfield[movie] > zfloor);
+            
+    if (draw)
+    {
+        if (AMPLITUDE_HIGH_RES > 0)
+        {
+            z_mid = compute_interpolated_colors_wave(i, j, xy_in, wave, palette, cplot, 
+                                                        rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie);
+            ij_to_xy(i, j, xy_sw);
+            ij_to_xy(i+1, j, xy_se);
+            ij_to_xy(i, j+1, xy_nw);
+            ij_to_xy(i+1, j+1, xy_ne);
+                    
+            for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]);
+                       
+            if (AMPLITUDE_HIGH_RES == 1)
+            {
+                glBegin(GL_TRIANGLE_FAN);
+                glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
+                draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
+            
+                glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
+                draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
+                
+                glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
+                draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
+                    
+                glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
+                draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
+                glEnd ();
+            }
+            else /* experimental */
+            {
+                glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
+                glBegin(GL_TRIANGLE_STRIP);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
+                draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
+                glEnd ();
+                    
+                glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
+                glBegin(GL_TRIANGLE_STRIP);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
+                draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
+                glEnd ();
+                    
+                glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
+                glBegin(GL_TRIANGLE_STRIP);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
+                draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
+                glEnd ();
+                    
+                glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
+                glBegin(GL_TRIANGLE_STRIP);
+                draw_vertex_xyz(xy_mid, z_mid);
+                draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
+                draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
+                glEnd ();
+            }
+        }
+        else
+        {
+            glColor3f(wave[i*NY+j].rgb[0], wave[i*NY+j].rgb[1], wave[i*NY+j].rgb[2]);
+            glBegin(GL_TRIANGLE_FAN);
+            ij_to_xy(i, j, xy);
+            draw_vertex_xyz(xy, *wave[i*NY+j].p_zfield[movie]);
+            ij_to_xy(i+1, j, xy);
+            draw_vertex_xyz(xy, *wave[(i+1)*NY+j].p_zfield[movie]);
+            ij_to_xy(i+1, j+1, xy);
+            draw_vertex_xyz(xy, *wave[(i+1)*NY+j+1].p_zfield[movie]);
+            ij_to_xy(i, j+1, xy);
+            draw_vertex_xyz(xy, *wave[i*NY+j+1].p_zfield[movie]);
+            glEnd ();
+        }
+    }
+    
+    if ((DRAW_OUTSIDE_GRAY)&&((!xy_in[i*NY+j])))
+    {
+        glColor3f(0.5, 0.5, 0.5);
+        glBegin(GL_TRIANGLE_FAN);
+        ij_to_xy(i, j, xy);
+        draw_vertex_xyz(xy, 0.0);
+        ij_to_xy(i+1, j, xy);
+        draw_vertex_xyz(xy, 0.0);
+        ij_to_xy(i+1, j+1, xy);
+        draw_vertex_xyz(xy, 0.0);
+        ij_to_xy(i, j+1, xy);
+        draw_vertex_xyz(xy, 0.0);
+        glEnd ();
+    }
+}
+
+
+void draw_wave_3d(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], 
+                  int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
+{
+    int i, j;
+    double observer_angle;
+    
     blank();
     if (DRAW_BILLIARD) draw_billiard_3d(fade, fade_value);
             
@@ -1343,111 +1513,44 @@ void draw_wave_3d(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_
         compute_cfield(xy_in, cplot, palette, wave, fade, fade_value, movie);
     }
     
-    for (i=0; i<NX-2; i++)
-        for (j=0; j<NY-2; j++)
+    if (!ROTATE_VIEW)
+    {
+        for (i=0; i<NX-2; i++)
+            for (j=0; j<NY-2; j++)
+                draw_wave_3d_ij(i, j, movie, phi, psi, xy_in, wave, zplot, cplot, palette, fade, fade_value);
+    }
+    else    /* draw facets in an order depending on the position of the observer */
+    {
+        observer_angle = argument(observer[0], observer[1]);
+        observer_angle -= 0.5*PID;
+        if (observer_angle < 0.0) observer_angle += DPI;
+        printf("Observer_angle = %.3lg\n", observer_angle*360.0/DPI); 
+        
+        if ((observer_angle > 0.0)&&(observer_angle < PID))
         {
-            if (NON_DIRICHLET_BC) 
-                draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
-            else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
-            
-            if (FLOOR_ZCOORD) 
-                draw = (draw)&&(*wave[i*NY+j].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j].p_zfield[movie] > zfloor)&&(*wave[i*NY+j+1].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j+1].p_zfield[movie] > zfloor);
-            
-            if (draw)
-            {
-                if (AMPLITUDE_HIGH_RES > 0)
-                {
-                    z_mid = compute_interpolated_colors_wave(i, j, xy_in, wave, palette, cplot, 
-                                                             rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie);
-                    
-                    ij_to_xy(i, j, xy_sw);
-                    ij_to_xy(i+1, j, xy_se);
-                    ij_to_xy(i, j+1, xy_nw);
-                    ij_to_xy(i+1, j+1, xy_ne);
-                    
-                    for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]);
-                       
-                    if (AMPLITUDE_HIGH_RES == 1)
-                    {
-                        glBegin(GL_TRIANGLE_FAN);
-                        glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
-                        draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
-                    
-                        glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
-                        draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
-                    
-                        glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
-                        draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
-                    
-                        glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
-                        draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
-                        glEnd ();
-                    }
-                    else /* experimental */
-                    {
-                        glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
-                        glBegin(GL_TRIANGLE_STRIP);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
-                        draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
-                        glEnd ();
-                    
-                        glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
-                        glBegin(GL_TRIANGLE_STRIP);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]);
-                        draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
-                        glEnd ();
-                    
-                        glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
-                        glBegin(GL_TRIANGLE_STRIP);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]);
-                        draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
-                        glEnd ();
-                    
-                        glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
-                        glBegin(GL_TRIANGLE_STRIP);
-                        draw_vertex_xyz(xy_mid, z_mid);
-                        draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]);
-                        draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]);
-                        glEnd ();
-                    }
-                }
-                else
-                {
-                    glColor3f(wave[i*NY+j].rgb[0], wave[i*NY+j].rgb[1], wave[i*NY+j].rgb[2]);
-    
-                    glBegin(GL_TRIANGLE_FAN);
-                    ij_to_xy(i, j, xy);
-                    draw_vertex_xyz(xy, *wave[i*NY+j].p_zfield[movie]);
-                    ij_to_xy(i+1, j, xy);
-                    draw_vertex_xyz(xy, *wave[(i+1)*NY+j].p_zfield[movie]);
-                    ij_to_xy(i+1, j+1, xy);
-                    draw_vertex_xyz(xy, *wave[(i+1)*NY+j+1].p_zfield[movie]);
-                    ij_to_xy(i, j+1, xy);
-                    draw_vertex_xyz(xy, *wave[i*NY+j+1].p_zfield[movie]);
-                    glEnd ();
-                }
-            }
-            
-            if ((DRAW_OUTSIDE_GRAY)&&((!xy_in[i*NY+j])))
-            {
-                glColor3f(0.5, 0.5, 0.5);
-                glBegin(GL_TRIANGLE_FAN);
-                ij_to_xy(i, j, xy);
-                draw_vertex_xyz(xy, 0.0);
-                ij_to_xy(i+1, j, xy);
-                draw_vertex_xyz(xy, 0.0);
-                ij_to_xy(i+1, j+1, xy);
-                draw_vertex_xyz(xy, 0.0);
-                ij_to_xy(i, j+1, xy);
-                draw_vertex_xyz(xy, 0.0);
-                glEnd ();
-            }
+            for (j=0; j<NY-2; j++)
+                for (i=0; i<NX-2; i++)
+                    draw_wave_3d_ij(i, j, movie, phi, psi, xy_in, wave, zplot, cplot, palette, fade, fade_value);
         }
+        else if (observer_angle < PI)
+        {
+            for (i=NX-3; i>0; i--)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_3d_ij(i, j, movie, phi, psi, xy_in, wave, zplot, cplot, palette, fade, fade_value);
+        }
+        else if (observer_angle < 1.5*PI)
+        {
+             for (j=NY-3; j>0; j--)
+                for (i=0; i<NX-2; i++)
+                    draw_wave_3d_ij(i, j, movie, phi, psi, xy_in, wave, zplot, cplot, palette, fade, fade_value);   
+        }
+        else
+        {
+            for (i=0; i<NX-2; i++)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_3d_ij(i, j, movie, phi, psi, xy_in, wave, zplot, cplot, palette, fade, fade_value);
+        }
+    }
         
     if (DRAW_BILLIARD_FRONT) draw_billiard_3d_front(fade, fade_value);
 }
diff --git a/sub_wave_3d_rde.c b/sub_wave_3d_rde.c
index 2158c0d..74bbfc6 100644
--- a/sub_wave_3d_rde.c
+++ b/sub_wave_3d_rde.c
@@ -21,7 +21,7 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2])
     static double n2, m2, d, sm2, sn2, v[3], h[2], plane_ratio = 0.5;
     static int first = 1;
     
-    if ((first)&&(REPRESENTATION_3D == REP_PROJ_3D))
+    if (((first)&&(REPRESENTATION_3D == REP_PROJ_3D))||(reset_view))
     {
         m2 = observer[0]*observer[0] + observer[1]*observer[1];
         n2 = m2 + observer[2]*observer[2];
@@ -34,8 +34,9 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2])
         v[1] = -observer[1]*observer[2]/(sn2*sm2);
         v[2] = m2/(sn2*sm2);
         first = 0;
-        printf("h = (%.3lg, %.3lg)\n", h[0], h[1]);
-        printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]);
+        reset_view = 0;
+//         printf("h = (%.3lg, %.3lg)\n", h[0], h[1]);
+//         printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]);
     }
     
     switch (REPRESENTATION_3D) {
@@ -154,14 +155,24 @@ void draw_billiard_3d(int fade, double fade_value)      /* draws the billiard bo
     glEnable(GL_LINE_SMOOTH);
 
     switch (B_DOMAIN) {
+        case (D_RECTANGLE):
+        {
+            glBegin(GL_LINE_LOOP);
+            draw_vertex_x_y_z(LAMBDA, -1.0, 0.0);
+            draw_vertex_x_y_z(LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, -1.0, 0.0);
+            glEnd();
+            break;
+        }
         case (D_ELLIPSE):
         {
             glBegin(GL_LINE_LOOP);
             for (i=0; i<=NSEG; i++)
             {
                 phi = (double)i*DPI/(double)NSEG;
-                x = LAMBDA*cos(phi);
-                y = sin(phi);
+                x = 1.05*LAMBDA*cos(phi);
+                y = 1.05*sin(phi);
                 draw_vertex_x_y_z(x, y, 0.0);
             }
             glEnd ();
@@ -205,6 +216,19 @@ void draw_billiard_3d(int fade, double fade_value)      /* draws the billiard bo
             draw_circle_3d(0.0, 0.0, LAMBDA, NSEG);
             break;
         }
+        case (D_POLYGON):
+        {
+            omega = DPI/((double)NPOLY);
+            glBegin(GL_LINE_LOOP);
+            for (i=0; i<=NPOLY; i++)
+            {
+                x = 1.0075*cos(i*omega + APOLY*PID);
+                y = 1.0075*sin(i*omega + APOLY*PID);
+                draw_vertex_x_y_z(x, y, 0.0);
+            }
+            glEnd ();
+            break;
+        }
         case (D_YOUNG):
         {
             if (FILL_BILLIARD_COMPLEMENT)
@@ -281,6 +305,14 @@ void draw_billiard_3d(int fade, double fade_value)      /* draws the billiard bo
                 y = (LAMBDA + 0.01)*sin(phi);
                 draw_vertex_x_y_z(x, y, 0.0);
             }
+            phi = PI - alpha;
+            x = 1.0 + (LAMBDA + 0.01)*cos(phi);
+            y = 0.01 + (LAMBDA + 0.01)*sin(phi);
+            draw_vertex_x_y_z(x, y, 0.0);
+            phi = alpha;
+            x = -1.0 + (LAMBDA + 0.01)*cos(phi);
+            y = 0.01 + (LAMBDA + 0.01)*sin(phi);
+            draw_vertex_x_y_z(x, y, 0.0);
             for (i=0; i<=NSEG; i++)
             {
                 phi = alpha + (double)i*dphi;
@@ -771,6 +803,15 @@ void draw_billiard_3d_front(int fade, double fade_value)
     glEnable(GL_LINE_SMOOTH);
 
     switch (B_DOMAIN) {
+        case (D_RECTANGLE):
+        {
+            glBegin(GL_LINE_STRIP);
+            draw_vertex_x_y_z(LAMBDA, -1.0, 0.0);
+            draw_vertex_x_y_z(LAMBDA, 1.0, 0.0);
+            draw_vertex_x_y_z(-LAMBDA, 1.0, 0.0);
+            glEnd();
+            break;
+        }
         case (D_STADIUM):
         {
             glBegin(GL_LINE_STRIP);
@@ -791,6 +832,19 @@ void draw_billiard_3d_front(int fade, double fade_value)
             glEnd();
             break;
         }
+        case (D_POLYGON):
+        {
+            omega = DPI/((double)NPOLY);
+            glBegin(GL_LINE_STRIP);
+            for (i=0; i<=NPOLY; i++)
+            {
+                x = cos(i*omega + APOLY*PID);
+                y = sin(i*omega + APOLY*PID);
+                draw_vertex_visible(x, y, 0.0);
+            }
+            glEnd ();
+            break;
+        }
         case (D_YOUNG):
         {
             glBegin(GL_LINE_STRIP);
@@ -1399,9 +1453,16 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in
             case (Z_ARGUMENT):
             {
                 value = dy_phase*(double)(j - jmin);
+//                 hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette);
                 color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
                 break;
             }
+            case (Z_REALPART):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
+            }
        }
         if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
         glColor3f(rgb[0], rgb[1], rgb[2]);
diff --git a/wave_3d.c b/wave_3d.c
index f6a7057..e283aff 100644
--- a/wave_3d.c
+++ b/wave_3d.c
@@ -50,15 +50,11 @@
 
 // #define WINWIDTH 	1920  /* window width */
 // #define WINHEIGHT 	1000  /* window height */
-// // #define NX 900          /* number of grid points on x axis */
-// // #define NY 450         /* number of grid points on y axis */
 // #define NX 1800          /* number of grid points on x axis */
 // #define NY 900          /* number of grid points on y axis */
-// // #define NX 1920          /* number of grid points on x axis */
-// // #define NY 1000          /* number of grid points on y axis */
 // 
-// // #define XMIN -2.0
-// // #define XMAX 2.0	/* x interval  */
+// #define XMIN -2.2
+// #define XMAX 1.8	/* x interval  */
 // #define YMIN -1.041666667
 // #define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
 
@@ -66,13 +62,10 @@
 
 #define WINWIDTH 	1280  /* window width */
 #define WINHEIGHT 	720   /* window height */
-
+ 
 #define NX 1280          /* number of grid points on x axis */
 #define NY 720          /* number of grid points on y axis */
 
-// #define NX 640          /* number of grid points on x axis */
-// #define NY 360          /* number of grid points on y axis */
- 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval  */
 #define YMIN -1.125
@@ -82,20 +75,20 @@
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 44       /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 24       /* choice of domain shape, see list in global_pdes.c */
 
-#define CIRCLE_PATTERN 202   /* pattern of circles or polygons, see list in global_pdes.c */
+#define CIRCLE_PATTERN 202    /* pattern of circles or polygons, see list in global_pdes.c */
 
-#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
+#define VARIABLE_IOR 1      /* set to 1 for a variable index of refraction */
 #define IOR 1               /* choice of index of refraction, see list in global_pdes.c */
-#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependece of IoR on Mandelbrot escape speed */
+#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
-#define MU 0.5             /* parameter controlling the dimensions of domain */
+#define LAMBDA 2.0	    /* parameter controlling the dimensions of domain */
+#define MU 0.5              /* parameter controlling the dimensions of domain */
 #define NPOLY 3             /* number of sides of polygon */
 #define APOLY 2.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 3            /* depth of computation of Menger gasket */
@@ -128,11 +121,11 @@
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 
 #define OMEGA 0.017         /* frequency of periodic excitation */
-#define AMPLITUDE 0.9      /* amplitude of periodic excitation */ 
-#define COURANT 0.1        /* Courant number */
-#define COURANTB 0.05      /* Courant number in medium B */
-#define GAMMA 0.0          /* damping factor in wave equation */
-#define GAMMAB 2.0e-5           /* damping factor in wave equation */
+#define AMPLITUDE 0.9       /* amplitude of periodic excitation */ 
+#define COURANT 0.08        /* Courant number */
+#define COURANTB 0.025      /* 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 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 */
@@ -152,8 +145,8 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1500        /* number of frames of movie */
-#define NVID 8            /* number of iterations between images displayed on screen */
+#define NSTEPS 2100        /* number of frames of movie */
+#define NVID 6            /* 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 1    /* width of billiard boundary */
@@ -165,6 +158,7 @@
 #define MID_FRAMES 100    /* number of still frames between parts of two-part movie */
 #define END_FRAMES 100   /* number of still frames at end of movie */
 #define FADE 1           /* set to 1 to fade at end of movie */
+#define ROTATE_VIEW_WHILE_FADE 0    /* set to 1 to keep rotating viewpoint during fade */
 
 /* Parameters of initial condition */
 
@@ -174,18 +168,17 @@
 
 /* Plot type, see list in global_pdes.c  */
 
-#define ZPLOT 108     /* wave height */
-#define CPLOT 108     /* color scheme */
+#define ZPLOT 103     /* wave height */
+#define CPLOT 103     /* color scheme */
 
-#define ZPLOT_B 109        
+#define ZPLOT_B 109 
 #define CPLOT_B 109        /* plot type for second movie */
 
-
 #define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of plot */
 #define SHADE_3D 1              /* set to 1 to change luminosity according to normal vector */
 #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 */
-#define DRAW_BILLIARD 1         /* set to 1 to draw boundary */
+#define FLOOR_ZCOORD 0          /* set to 1 to draw only facets with z not too negative */
+#define DRAW_BILLIARD 0         /* set to 1 to draw boundary */
 #define DRAW_BILLIARD_FRONT 0   /* set to 1 to draw front of boundary after drawing wave */
 #define DRAW_CONSTRUCTION_LINES 0   /* set to 1 to draw construction lines of certain domains */
 #define FADE_IN_OBSTACLE 1      /* set to 1 to fade color inside obstacles */
@@ -201,25 +194,28 @@
 #define REP_AXO_3D 0        /* linear projection (axonometry) */
 #define REP_PROJ_3D 1       /* projection on plane orthogonal to observer line of sight */
 
+#define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
+#define ROTATE_ANGLE 540.0  /* total angle of rotation during simulation */
+
 /* Color schemes */
 
-#define COLOR_PALETTE 14      /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 11     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 11      /* 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 */
 
 #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 VSCALE_AMPLITUDE 1.0   /* additional scaling factor for color scheme P_3D_AMPLITUDE */
-#define VSCALE_ENERGY 10.0     /* additional scaling factor for color scheme P_3D_ENERGY */
+#define SLOPE 2.0        /* sensitivity of color on wave amplitude */
+#define VSCALE_AMPLITUDE 0.2   /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSCALE_ENERGY 7.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 ATTENUATION 0.0    /* exponential attenuation coefficient of contrast with time */
 #define E_SCALE 1000.0     /* scaling factor for energy representation */
-#define LOG_SCALE 0.25    /* scaling factor for energy log representation */
-#define LOG_SHIFT 0.0    /* shift of colors on log scale */
+#define LOG_SCALE 0.2      /* scaling factor for energy log representation */
+#define LOG_SHIFT -0.5     /* shift of colors on log scale */
 #define LOG_ENERGY_FLOOR -1000.0    /* floor value for log of (total) energy */
 #define LOG_MEAN_ENERGY_SHIFT 5.0   /* additional shift for log of mean energy */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@@ -229,12 +225,12 @@
 #define LUMMEAN 0.5      /* amplitude of luminosity variation for scheme C_LUM */
 #define LUMAMP 0.3       /* amplitude of luminosity variation for scheme C_LUM */
 #define HUEMEAN 240.0    /* mean value of hue for color scheme C_HUE */
-#define HUEAMP -200.0      /* amplitude of variation of hue for color scheme C_HUE */
+#define HUEAMP -200.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 20.0      /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 100.0    /* scale of color scheme bar for 2nd part */
-#define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
+#define DRAW_COLOR_SCHEME 1       /* set to 1 to plot the color scheme */
+#define COLORBAR_RANGE 2.5      /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 50.0    /* scale of color scheme bar for 2nd part */
+#define ROTATE_COLOR_SCHEME 0     /* set to 1 to draw color scheme horizontally */
 
 #define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
 
@@ -248,13 +244,14 @@ double u_3d[2] = {0.75, -0.45};     /* projections of basis vectors for REP_AXO_
 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] = {6.0, 8.0, 6.0};    /* location of observer for REP_PROJ_3D representation */ 
+double observer[3] = {8.0, 8.0, 8.0};    /* location of observer for REP_PROJ_3D representation */ 
+int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
-#define Z_SCALING_FACTOR 0.05     /* 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 ZMAX_FACTOR 1.0           /* max value of z coordinate for REP_PROJ_3D representation */
-#define XSHIFT_3D -0.1           /* overall x shift for REP_PROJ_3D representation */
-#define YSHIFT_3D 0.15           /* overall y shift for REP_PROJ_3D representation */
+#define Z_SCALING_FACTOR 0.05   /* 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.1          /* overall x shift for REP_PROJ_3D representation */
+#define YSHIFT_3D 0.0           /* overall y shift for REP_PROJ_3D representation */
 
 
 #include "global_pdes.c"        /* constants and global variables */
@@ -519,6 +516,28 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
         draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
 }
 
+void viewpoint_schedule(int i)
+/* change position of observer */
+{
+    int j;
+    double angle, ca, sa;
+    static double observer_initial[3];
+    static int first = 1;
+    
+    if (first)
+    {
+        for (j=0; j<3; j++) observer_initial[j] = observer[j];
+        first = 0;
+    }
+    
+    angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
+    ca = cos(angle);
+    sa = sin(angle);
+    observer[0] = ca*observer_initial[0] - sa*observer_initial[1];
+    observer[1] = sa*observer_initial[0] + ca*observer_initial[1];
+    printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
+}
+
 
 void animation()
 {
@@ -589,7 +608,8 @@ void animation()
 
 //     init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in);
     init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in);
-//     add_circular_wave(-1.0, -0.2, -0.4, phi, psi, xy_in);
+//     init_wave_flat_mod(phi, psi, xy_in);
+//     add_circular_wave_mod(1.0, 1.0, 0.0, phi, psi, xy_in);
 
 //     printf("Wave initialized\n");
 
@@ -634,6 +654,12 @@ void animation()
         }
         else scale = 1.0;
         
+        if (ROTATE_VIEW) 
+        {
+            viewpoint_schedule(i - INITIAL_TIME);
+            reset_view = 1;
+        }
+        
         draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
         for (j=0; j<NVID; j++) 
         {
@@ -702,6 +728,11 @@ void animation()
             if (!FADE) for (i=0; i<MID_FRAMES; i++) save_frame();
             else for (i=0; i<MID_FRAMES; i++) 
             {
+                if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
+                {
+                    viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
+                    reset_view = 1;
+                }
                 fade_value = 1.0 - (double)i/(double)MID_FRAMES;
                 draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);   
@@ -715,6 +746,11 @@ void animation()
             if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
             else for (i=0; i<END_FRAMES; i++) 
             {
+                if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE)) 
+                {
+                    viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
+                    reset_view = 1;
+                }
                 fade_value = 1.0 - (double)i/(double)END_FRAMES;
                 draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);   
diff --git a/wave_billiard.c b/wave_billiard.c
index 4f1936d..97508d5 100644
--- a/wave_billiard.c
+++ b/wave_billiard.c
@@ -41,6 +41,7 @@
 #include <sys/types.h>
 #include <tiffio.h>     /* Sam Leffler's libtiff library. */
 #include <omp.h>
+#include <time.h>
 
 #define MOVIE 0         /* set to 1 to generate movie */
 #define DOUBLE_MOVIE 0  /* set to 1 to produce movies for wave height and energy simultaneously */
@@ -49,17 +50,17 @@
 
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1000  /* window height */
-// #define NX 1920          /* number of grid points on x axis */
-// #define NY 1000          /* number of grid points on y axis */
-#define NX 3840          /* number of grid points on x axis */
-#define NY 2000          /* number of grid points on y axis */
+#define NX 1920          /* number of grid points on x axis */
+#define NY 1000          /* number of grid points on y axis */
+// #define NX 3840          /* number of grid points on x axis */
+// #define NY 2000          /* number of grid points on y axis */
 
-#define XMIN -1.25
-#define XMAX 2.75	/* x interval  */
+#define XMIN -2.0
+#define XMAX 2.0	/* x interval  */
 #define YMIN -1.041666667
 #define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
 
-#define HIGHRES 1        /* set to 1 if resolution of grid is double that of displayed image */
+#define HIGHRES 0        /* set to 1 if resolution of grid is double that of displayed image */
 
 // #define WINWIDTH 	1280  /* window width */
 // #define WINHEIGHT 	720   /* window height */
@@ -69,8 +70,8 @@
 // #define NX 2560          /* number of grid points on x axis */
 // #define NY 1440          /* number of grid points on y axis */
 // 
-// #define XMIN -1.25
-// #define XMAX 2.75	/* x interval  */
+// #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 */
 
@@ -159,6 +160,7 @@
 #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 FADE 1           /* set to 1 to fade at end of movie */
 
 /* Parameters of initial condition */
 
@@ -169,7 +171,7 @@
 /* Plot type, see list in global_pdes.c  */
 
 #define PLOT 0
-// #define PLOT 3
+// #define PLOT 0
 // #define PLOT 1
 
 #define PLOT_B 3        /* plot type for second movie */
@@ -491,18 +493,30 @@ void draw_color_bar(int plot, double range)
 //     else draw_color_scheme(1.7, YMIN + 0.25, 1.9, YMAX - 0.25, plot, -range, range);
 }
 
-void draw_color_bar_palette(int plot, double range, int palette)
+// void draw_color_bar_palette(int plot, double range, int palette)
+// {
+//     if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette);
+//     else draw_color_scheme_palette(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range, palette);
+// }
+
+void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
 {
-    if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette);
-    else draw_color_scheme_palette(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range, palette);
+    double width = 0.14;
+//     double width = 0.2;
+    
+    if (ROTATE_COLOR_SCHEME) 
+        draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
+    else 
+        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, r2, xy[2]; 
+    double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; 
     double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX], *color_scale[NX];
     short int *xy_in[NX];
-    int i, j, s, sample_left[2], sample_right[2], period = 0;
+    int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
     static int counter = 0;
     long int wave_value;
     
@@ -569,6 +583,9 @@ 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);
+
+//     init_circular_wave(sqrt(LAMBDA*LAMBDA - 1.0), 0.0, phi, psi, xy_in);
+//     init_circular_wave(0.0, 0.0, phi, psi, xy_in);
     
 //     init_wave_plus(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
 //     init_wave(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
@@ -607,12 +624,12 @@ void animation()
     glColor3f(0.0, 0.0, 0.0);
 //     draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
     if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, 1.0, 0, PLOT, COLOR_PALETTE);
-    else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE);
+    else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
 
-    draw_billiard();
+    draw_billiard(0, 1.0);
     
     
-    if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE);
+    if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
 
     glutSwapBuffers();
 
@@ -634,7 +651,7 @@ void animation()
 
 //         draw_wave(phi, psi, xy_in, scale, i, PLOT);
         if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT, COLOR_PALETTE);
-        else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE);
+        else draw_wave_epalette(phi, psi, total_energy, 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);
@@ -650,9 +667,9 @@ void animation()
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
-        draw_billiard();
+        draw_billiard(0, 1.0);
         
-        if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE); 
+        if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value); 
         
         /* add oscillating waves */
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
@@ -673,9 +690,9 @@ void animation()
             {
 //                 draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
                 if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B);
-                else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B);
-                draw_billiard();
-                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B);  
+                else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
+                draw_billiard(0, 1.0);
+                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);  
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
 //                 save_frame_counter(NSTEPS + 21 + counter);
@@ -700,22 +717,42 @@ void animation()
         {
 //             draw_wave(phi, psi, xy_in, scale, i, PLOT);
             if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE);
-            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE);
-            draw_billiard();
-            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE);   
+            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
+            draw_billiard(0, 1.0);
+            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);   
             glutSwapBuffers();
         }
-        for (i=0; i<MID_FRAMES; i++) save_frame();
+        if (!FADE) for (i=0; i<MID_FRAMES; i++) save_frame();
+        else for (i=0; i<MID_FRAMES; i++) 
+        {
+            fade_value = 1.0 - (double)i/(double)MID_FRAMES;
+            if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE);
+            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
+            draw_billiard(1, fade_value);
+            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);   
+            glutSwapBuffers();
+            save_frame_counter(NSTEPS + i + 1);
+        }
         if (DOUBLE_MOVIE) 
         {
 //             draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
             if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B);
-            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B);
-            draw_billiard();
-            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B); 
+            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
+            draw_billiard(0, 1.0);
+            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); 
             glutSwapBuffers();
         }
-        for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+        if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+        else for (i=0; i<END_FRAMES; i++) 
+        {
+            fade_value = 1.0 - (double)i/(double)END_FRAMES;
+            if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B);
+            else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
+            draw_billiard(1, fade_value);
+            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);   
+            glutSwapBuffers();
+            save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+        }
         
         s = system("mv wave*.tif tif_wave/");
     }
@@ -742,6 +779,12 @@ void animation()
 
 void display(void)
 {
+    time_t rawtime;
+    struct tm * timeinfo;
+
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    
     glPushMatrix();
 
     blank();
@@ -755,7 +798,11 @@ void display(void)
     glPopMatrix();
 
     glutDestroyWindow(glutGetWindow());
-
+    
+    printf("Start local time and date: %s", asctime(timeinfo));
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    printf("Current local time and date: %s", asctime(timeinfo));
 }
 
 
diff --git a/wave_common.c b/wave_common.c
index 831274c..5b198ff 100644
--- a/wave_common.c
+++ b/wave_common.c
@@ -683,10 +683,10 @@ void draw_wave_highres_diss(int size, double *phi[NX], double *psi[NX], double *
 
 
 void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], 
-                        double scale, int time, int plot, int palette)
+                        double scale, int time, int plot, int palette, int fade, double fade_value)
 /* same as draw_wave_e, but with color scheme specification */
 {
-    int i, j, iplus, iminus, jplus, jminus;
+    int i, j, k, iplus, iminus, jplus, jminus;
     double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2;
     static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
 
@@ -756,6 +756,7 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
                         break;
                     }
                 }
+                if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
                 glColor3f(rgb[0], rgb[1], rgb[2]);
 
                 glVertex2i(i, j);