From f56b8a0ab488ed724734bb31883e38accaa8f679 Mon Sep 17 00:00:00 2001
From: nilsberglund-orleans
 <83530463+nilsberglund-orleans@users.noreply.github.com>
Date: Sun, 18 Jul 2021 15:06:59 +0200
Subject: [PATCH] Add files via upload

---
 drop_billiard.c     |  68 +++++-----
 global_particles.c  |  33 +++++
 global_pdes.c       |  90 ++++++++++++++
 heat.c              | 209 ++++++++++++++++++-------------
 particle_billiard.c | 215 ++++++++++++++++++++++----------
 schrodinger.c       | 179 ++++++++++-----------------
 sub_wave.c          | 294 +++++++++++++++++++++++++++++++++++++++++++-
 wave_billiard.c     | 288 +++++++++++++++++++++++--------------------
 8 files changed, 951 insertions(+), 425 deletions(-)
 create mode 100644 global_particles.c
 create mode 100644 global_pdes.c

diff --git a/drop_billiard.c b/drop_billiard.c
index a45d841..9a13700 100644
--- a/drop_billiard.c
+++ b/drop_billiard.c
@@ -33,32 +33,28 @@
 #define WINWIDTH 	1280  /* window width */
 #define WINHEIGHT 	720   /* window height */
 
-#define XMIN -1.8
-#define XMAX 1.8	/* x interval */
-#define YMIN -0.91
-#define YMAX 1.115	/* y interval for 9/16 aspect ratio */
-// #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 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 XMIN -1.8
+// #define XMAX 1.8	/* x interval */
+// #define YMIN -0.91
+// #define YMAX 1.115	/* y interval for 9/16 aspect ratio */
 
 #define SCALING_FACTOR 1.0       /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */
 
-/* Choice of the billiard table */
+/* Choice of the billiard table, see global_particles.c  */
 
-#define B_DOMAIN 9      /* choice of domain shape */
+#define B_DOMAIN 13      /* choice of domain shape */
 
-#define D_RECTANGLE 0   /* rectangular domain */
-#define D_ELLIPSE 1     /* elliptical domain */
-#define D_STADIUM 2     /* stadium-shaped domain */
-#define D_SINAI 3       /* Sinai billiard */
-#define D_DIAMOND 4     /* diamond-shaped billiard */
-#define D_TRIANGLE 5    /* triangular billiard */
-#define D_ANNULUS 7     /* annulus */
-#define D_POLYGON 8     /* polygon */
-#define D_REULEAUX 9    /* Reuleaux and star shapes */
-#define D_FLOWER 10     /* Bunimovich flower */
-#define D_ALT_REU 11    /* alternating between star and Reuleaux */
+#define CIRCLE_PATTERN 0    /* pattern of circles */
+
+#define NMAXCIRCLES 1000        /* total number of circles (must be at least NCX*NCY for square grid) */
+// #define NCX 10            /* number of circles in x direction */
+// #define NCY 15            /* number of circles in y direction */
+#define NCX 15            /* number of circles in x direction */
+#define NCY 20            /* number of circles in y direction */
 
 // #define LAMBDA 1.0	/* parameter controlling shape of billiard */
 #define LAMBDA 1.124950941	/* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */
@@ -72,16 +68,17 @@
 #define APOLY -1.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define DRAW_BILLIARD 1     /* set to 1 to draw billiard */
 #define DRAW_CONSTRUCTION_LINES 1   /* set to 1 to draw additional construction lines for billiard */
+#define PERIODIC_BC 0       /* set to 1 to enforce periodic boundary conditions when drawing particles */
 
 #define RESAMPLE 0      /* set to 1 if particles should be added when dispersion too large */
 
-#define NPART 100000	/* number of particles */
+#define NPART 50000	/* number of particles */
 #define NPARTMAX 100000	/* maximal number of particles after resampling */
 
-#define NSTEPS 4000         /* number of frames of movie */
-#define TIME 15             /* time between movie frames, for fluidity of real-time simulation */ 
+#define NSTEPS 5500         /* number of frames of movie */
+#define TIME 40             /* time between movie frames, for fluidity of real-time simulation */ 
 #define DPHI 0.0001         /* integration step */
-#define NVID 10             /* number of iterations between images displayed on screen */
+#define NVID 20             /* number of iterations between images displayed on screen */
 
 /* Decreasing TIME accelerates the animation and the movie               */
 /* For constant speed of movie, TIME*DPHI should be kept constant        */
@@ -91,16 +88,17 @@
 /* simulation parameters */
 
 #define LMAX 0.01       /* minimal segment length triggering resampling */ 
-#define LPERIODIC 1.0   /* lines longer than this are not drawn (useful for Sinai billiard) */
-#define LCUT 1000.0        /* controls the max size of segments not considered as being cut */
+#define LPERIODIC 0.1   /* lines longer than this are not drawn (useful for Sinai billiard) */
+#define LCUT 2.0        /* controls the max size of segments not considered as being cut */
 #define DMIN 0.02       /* minimal distance to boundary for triggering resampling */ 
 #define CYCLE 0         /* set to 1 for closed curve (start in all directions) */
 #define ORDER_COLORS 1  /* set to 1 if colors should be drawn in order */ 
 
 /* color and other graphical parameters */
 
-#define NCOLORS 10          /* number of colors */
+#define NCOLORS 16          /* number of colors */
 #define COLORSHIFT 200      /* hue of initial color */ 
+#define RAINBOW_COLOR 1     /* set to 1 to use different colors for all particles */
 #define NSEG 100            /* number of segments of boundary */
 #define BILLIARD_WIDTH 4    /* width of billiard */
 #define FRONT_WIDTH 3       /* width of wave front */
@@ -120,6 +118,7 @@
 #define DPI 	6.283185307
 #define PID 	1.570796327
 
+#include "global_particles.c"
 #include "sub_part_billiard.c"
 
 
@@ -394,8 +393,11 @@ double *configs[NPARTMAX];
             if (configs[i][2]<0.0) 
             {    
                 vbilliard(configs[i]);
-                color[i]++;
-                if (color[i] >= NCOLORS) color[i] -= NCOLORS;
+                if (!RAINBOW_COLOR)
+                {
+                    color[i]++;
+                    if (color[i] >= NCOLORS) color[i] -= NCOLORS;
+                }
             }
 
             configs[i][2] += DPHI; 
@@ -449,7 +451,8 @@ void animation()
     for (i=0; i<NPARTMAX; i++)
         configs[i] = (double *)malloc(8*sizeof(double));
   
-    init_drop_config(0.0, 0.1, 0.0, DPI, configs);
+    init_drop_config(-1.0 + 0.3*sqrt(2.0), -1.0 + 0.5*sqrt(2.0), 0.0, DPI, configs);
+//     init_drop_config(-0.5, -0.5, 0.0, DPI, configs);
 //     init_boundary_config(1.5, 1.5, 0.0, PI, configs);
 
   
@@ -463,6 +466,9 @@ void animation()
   
   
     for (i=0; i<NPARTMAX; i++) color[i] = 0;
+    
+    if (RAINBOW_COLOR)      /* rainbow color scheme */
+        for (i=0; i<NPART; i++) color[i] = (i*NCOLORS)/NPART;
   
     sleep(SLEEP1);
   
diff --git a/global_particles.c b/global_particles.c
new file mode 100644
index 0000000..d830c8a
--- /dev/null
+++ b/global_particles.c
@@ -0,0 +1,33 @@
+double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES];      /* position and radius of circular scatters */
+short int circleactive[NMAXCIRCLES];                                      /* tells which circular scatters are active */
+int ncircles = NMAXCIRCLES;            /* actual number of circles, can be decreased e.g. for random patterns */
+
+
+/* some basic math */
+
+#define PI 	3.141592654
+#define DPI 	6.283185307
+#define PID 	1.570796327
+
+/* Choice of the billiard table */
+
+#define D_RECTANGLE 0   /* rectangular domain */
+#define D_ELLIPSE 1     /* elliptical domain */
+#define D_STADIUM 2     /* stadium-shaped domain */
+#define D_SINAI 3       /* Sinai billiard */
+#define D_DIAMOND 4     /* diamond-shaped billiard */
+#define D_TRIANGLE 5    /* triangular billiard */
+#define D_ANNULUS 7     /* annulus */
+#define D_POLYGON 8     /* polygon */
+#define D_REULEAUX 9    /* Reuleaux and star shapes */
+#define D_FLOWER 10     /* Bunimovich flower */
+#define D_ALT_REU 11    /* alternating between star and Reuleaux */
+#define D_ANGLE 12      /* angular sector */
+#define D_LSHAPE 13     /* L-shaped billiard for conical singularity */
+#define D_GENUSN 14     /* polygon with identifies opposite sides */
+
+#define D_CIRCLES 20     /* several circles */
+
+#define C_FOUR_CIRCLES 0  /* four circles almost touching each other */
+#define C_SQUARE 1        /* square grid of circles */
+#define C_HEX 2           /* hexagonal/triangular grid of circles */
diff --git a/global_pdes.c b/global_pdes.c
new file mode 100644
index 0000000..9af3673
--- /dev/null
+++ b/global_pdes.c
@@ -0,0 +1,90 @@
+/* Global variables and parameters for wave_billiard, heat and schrodinger */
+
+/* Basic math */
+
+#define PI 	3.141592654
+#define DPI 	6.283185307
+#define PID 	1.570796327
+
+/* shape of domain */
+
+#define D_RECTANGLE 0   /* rectangular domain */
+#define D_ELLIPSE 1     /* elliptical domain */
+#define D_STADIUM 2     /* stadium-shaped domain */
+#define D_SINAI 3       /* Sinai billiard */
+#define D_DIAMOND 4     /* diamond-shaped billiard */
+#define D_TRIANGLE 5    /* triangular billiard */
+#define D_FLAT 6        /* flat interface */
+#define D_ANNULUS 7     /* annulus */
+#define D_POLYGON 8     /* polygon */
+#define D_YOUNG 9       /* Young diffraction slits */
+#define D_GRATING 10    /* diffraction grating */
+#define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
+
+/* The following 3 types are superseded by D_CIRCLES and can be removed later */
+#define D_DISK_GRID 12  /* grid of disks */
+#define D_DISK_HEX 13   /* haxagonl grid of disks */
+#define D_DISK_PERCOL 14    /* random grid of percolation type */
+
+#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
+#define D_JULIA_INT 16  /* interior of Julia set */ 
+#define D_MENGER_ROTATED 17  /* rotated Menger-Sierpinski carpet */
+
+#define D_CIRCLES 20    /* several circles */
+
+#define NMAXCIRCLES 1000        /* total number of circles (must be at least NCX*NCY for square grid) */
+
+#define C_SQUARE 0          /* square grid of circles */
+#define C_HEX 1             /* hexagonal/triangular grid of circles */
+#define C_RAND_DISPLACED 2  /* randomly displaced square grid */
+#define C_RAND_PERCOL 3     /* random percolation arrangement */
+#define C_RAND_POISSON 4    /* random Poisson point process */
+
+/* Billiard tables for heat equation */
+
+#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
+#define D_MENGER_HEATED 22  /* Menger gasket with different temperatures */
+#define D_MENGER_H_OPEN 23  /* Menger gasket with different temperatures and larger domain */
+#define D_MANDELBROT 24     /* Mandelbrot set */
+#define D_JULIA 25          /* Julia set */
+#define D_MANDELBROT_CIRCLE 26     /* Mandelbrot set with circular conductor */
+
+/* Boundary conditions */
+
+#define BC_DIRICHLET 0   /* Dirichlet boundary conditions */
+#define BC_PERIODIC 1    /* periodic boundary conditions */
+#define BC_ABSORBING 2   /* absorbing boundary conditions (beta version) */
+#define BC_VPER_HABS 3   /* vertically periodic and horizontally absorbing boundary conditions */
+
+/* For debugging purposes only */
+#define FLOOR 0         /* set to 1 to limit wave amplitude to VMAX */
+#define VMAX 10.0       /* max value of wave amplitude */
+
+/* Plot types */
+
+/* For wave equation */
+#define P_AMPLITUDE 0    /* plot amplitude of wave */
+#define P_ENERGY 1       /* plot energy of wave */
+#define P_MIXED 2        /* plot amplitude in upper half, energy in lower half */
+
+/* For Schrodinger equation */
+#define P_MODULE 10        /* plot module of wave function squared */
+#define P_PHASE 11         /* plot phase of wave function */
+#define P_REAL 12          /* plot real part */
+#define P_IMAGINARY 13     /* plot imaginary part */
+
+
+/* Color schemes */
+
+#define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
+#define C_HUE 1          /* color scheme modifies hue */
+#define C_PHASE 2        /* color scheme shows phase */
+
+
+double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES];      /* position and radius of circular scatterers */
+short int circleactive[NMAXCIRCLES];                                      /* tells which circular scatters are active */
+int ncircles = NMAXCIRCLES;            /* actual number of circles, can be decreased e.g. for random patterns */
+
+double julia_x = -0.5, julia_y = 0.5;    /* parameters for Julia sets */
+// double julia_x = 0.33267, julia_y = 0.06395;    /* parameters for Julia sets */
+// double julia_x = 0.37468, julia_y = 0.21115;    /* parameters for Julia sets */
diff --git a/heat.c b/heat.c
index 375af8b..c50a33f 100644
--- a/heat.c
+++ b/heat.c
@@ -50,46 +50,25 @@
 /* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
 /* but will multiply run time by 4                                 */
 
-#define XMIN -2.0
-#define XMAX 2.0	/* x interval */
-// #define XMIN -1.5
-// #define XMAX 2.5	/* x interval */
+// #define XMIN -2.0
+// #define XMAX 2.0	/* x interval */
+#define XMIN -2.5
+#define XMAX 1.5	/* x interval */
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
-#define JULIA_SCALE 1.1 /* scaling for Julia sets */
-// #define JULIA_SCALE 0.8 /* scaling for Julia sets */
+#define JULIA_SCALE 1.0 /* scaling for Julia sets */
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 25      /* choice of domain shape */
+#define B_DOMAIN 26      /* choice of domain shape, see list in global_pdes.c  */
 
-#define D_RECTANGLE 0   /* rectangular domain */
-#define D_ELLIPSE 1     /* elliptical domain */
-#define D_STADIUM 2     /* stadium-shaped domain */
-#define D_SINAI 3       /* Sinai billiard */
-#define D_DIAMOND 4     /* diamond-shaped billiard */
-#define D_TRIANGLE 5    /* triangular billiard */
-#define D_FLAT 6        /* flat interface */
-#define D_ANNULUS 7     /* annulus */
-#define D_POLYGON 8     /* polygon */
-#define D_YOUNG 9       /* Young diffraction slits */
-#define D_GRATING 10    /* diffraction grating */
-#define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
+#define CIRCLE_PATTERN 0    /* pattern of circles, see list in global_pdes.c */
 
-#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
-#define D_JULIA_INT 16  /* interior of Julia set */ 
+#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 */
 
-/* Billiard tables for heat equation */
-
-#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
-#define D_MENGER_HEATED 22  /* Menger gasket with different temperatures */
-#define D_MENGER_H_OPEN 23  /* Menger gasket with different temperatures and larger domain */
-#define D_MANDELBROT 24     /* Mandelbrot set */
-#define D_JULIA 25          /* Julia set */
-#define D_MANDELBROT_CIRCLE 26     /* Mandelbrot set with circular conductor */
-
-#define LAMBDA 0.7	    /* parameter controlling the dimensions of domain */
+#define LAMBDA -1.0	    /* parameter controlling the dimensions of domain */
 #define MU 0.1	            /* 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 */
@@ -98,6 +77,8 @@
 #define MANDELLEVEL 1000      /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0     /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
+#define NGRIDX 15            /* number of grid point for grid of disks */
+#define NGRIDY 20           /* number of grid point for grid of disks */
 
 /* You can add more billiard tables by adapting the functions */
 /* xy_in_billiard and draw_billiard in sub_wave.c */
@@ -116,13 +97,9 @@
 // #define T_IN 2.0        /* inside temperature */
 #define SPEED 0.0       /* speed of drift to the right */
 
-/* Boundary conditions */
+/* Boundary conditions, see list in global_pdes.c  */
 
-#define B_COND 0
-
-#define BC_DIRICHLET 0   /* Dirichlet boundary conditions */
-#define BC_PERIODIC 1    /* periodic boundary conditions */
-#define BC_ABSORBING 2   /* absorbing boundary conditions (beta version) */
+#define B_COND 1
 
 /* Parameters for length and speed of simulation */
 
@@ -142,26 +119,22 @@
 
 /* Field representation */
 
-#define FIELD_REP 0
+#define FIELD_REP 1
 
 #define F_INTENSITY 0   /* color represents intensity */
 #define F_GRADIENT 1    /* color represents norm of gradient */ 
 
 #define DRAW_FIELD_LINES 1  /* set to 1 to draw field lines */
 #define FIELD_LINE_WIDTH 1  /* width of field lines */
-#define N_FIELD_LINES 200   /* number of field lines */
+#define N_FIELD_LINES 50   /* number of field lines */
 #define FIELD_LINE_FACTOR 100 /* factor controlling precision when computing origin of field lines */
 
-/* Color schemes */
+/* Color schemes, see list in global_pdes.c  */
 
 #define BLACK 1          /* black background */
 
 #define COLOR_SCHEME 1   /* choice of color scheme */
 
-#define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
-#define C_HUE 1          /* color scheme modifies hue */
-#define C_PHASE 2        /* color scheme shows phase */
-
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 // #define SLOPE 0.1        /* sensitivity of color on wave amplitude */
 #define SLOPE 0.3        /* sensitivity of color on wave amplitude */
@@ -171,20 +144,13 @@
 #define COLORDRIFT 0.0   /* how much the color hue drifts during the whole simulation */
 #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 280.0    /* mean value of hue for color scheme C_HUE */
-#define HUEAMP -110.0      /* amplitude of variation of hue for color scheme C_HUE */
+#define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
+#define HUEAMP -220.0      /* amplitude of variation of hue for color scheme C_HUE */
 // #define HUEMEAN 270.0    /* mean value of hue for color scheme C_HUE */
 // #define HUEAMP -130.0      /* amplitude of variation of hue for color scheme C_HUE */
 
-/* Basic math */
-
-#define PI 	3.141592654
-#define DPI 	6.283185307
-#define PID 	1.570796327
-
-double julia_x = 0.0, julia_y = 0.0;    /* parameters for Julia sets */
-
 
+#include "global_pdes.c"
 #include "sub_wave.c"
 
 double courant2;  /* Courant parameter squared */
@@ -320,12 +286,7 @@ short int *xy_in[NX];
             x1 = x1 + delta*nabx/norm;
             y1 = y1 + delta*naby/norm;
         }
-        else 
-        {
-            cont = 0;
-//             nablax[ij[0]][ij[1]] = 0.0;
-//             nablay[ij[0]][ij[1]] = 0.0;            
-        }
+        else cont = 0;
         
         if (!xy_in[ij[0]][ij[1]]) cont = 0;
         
@@ -361,7 +322,6 @@ int time;
     }
     
     /* compute the gradient */
-//     if (FIELD_REP > 0) 
     compute_gradient(phi, nablax, nablay);
     
     /* compute the position of origins of field lines */
@@ -371,25 +331,24 @@ int time;
         
         printf("computing linex\n");
         
-        x1 = sqrt(3.58);
-        y1 = 0.0;
+        x1 = LAMBDA + MU*1.01;
+        y1 = 1.0;
         linex[0] = x1;
         liney[0] = y1;
         dangle = DPI/((double)(N_FIELD_LINES*FIELD_LINE_FACTOR));
             
         for (i = 1; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
         {
-//             angle = PID + (double)i*dangle;
             angle = (double)i*dangle;
-            x2 = sqrt(3.58)*cos(angle);
-            y2 = sqrt(1.18)*sin(angle);
+            x2 = LAMBDA + MU*1.01*cos(angle);
+            y2 = 0.5 + MU*1.01*sin(angle);
             linex[i] = x2;
             liney[i] = y2;
             distance[i-1] = module2(x2-x1,y2-y1);
             x1 = x2;
             y1 = y2;
         }
-        distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-sqrt(3.58),y2);
+        distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-LAMBDA,y2-0.5);
     }
 
     dx = (XMAX-XMIN)/((double)NX);
@@ -404,7 +363,6 @@ int time;
                 value = module2(nablax[i][j], nablay[i][j]);
             }
             
-//             if ((phi[i][j] - T_IN)*(phi[i][j] - T_OUT) < 0.0) 
             if (xy_in[i][j] == 1) 
             {
                 color_scheme(COLOR_SCHEME, value, scale, time, rgb);
@@ -431,18 +389,14 @@ int time;
             else integral[i] = intens;
         }
         deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
-//         deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES + 1.0);
-//         deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
         
 //         printf("delta = %.5lg\n", deltaintens);
         
         i = 0;
-//         draw_field_line(linex[0], liney[0], nablax, nablay, 0.00002, 100000);
         draw_field_line(linex[0], liney[0], xy_in, nablax, nablay, 0.00002, 100000);
         for (j = 1; j < N_FIELD_LINES+1; j++)
         {
             while ((integral[i] <= j*deltaintens)&&(i < N_FIELD_LINES*FIELD_LINE_FACTOR)) i++; 
-//             draw_field_line(linex[i], liney[i], nablax, nablay, 0.00002, 100000);
             draw_field_line(linex[i], liney[i], xy_in, nablax, nablay, 0.00002, 100000);
             counter++;
         }
@@ -459,12 +413,101 @@ int time;
 
 
 
-void evolve_wave(phi, xy_in)
+void evolve_wave_half(phi_in, phi_out, xy_in)
+/* time step of field evolution */
+    double *phi_in[NX], *phi_out[NX]; short int *xy_in[NX];
+{
+    int i, j, iplus, iminus, jplus, jminus;
+    double delta1, delta2, x, y;
+    
+    #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
+    for (i=0; i<NX; i++){
+        for (j=0; j<NY; j++){
+            if (xy_in[i][j] == 1){
+                /* discretized Laplacian depending on boundary conditions */
+                if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
+                {
+                    iplus = (i+1);   if (iplus == NX) iplus = NX-1;
+                    iminus = (i-1);  if (iminus == -1) iminus = 0;
+                    jplus = (j+1);   if (jplus == NY) jplus = NY-1;
+                    jminus = (j-1);  if (jminus == -1) jminus = 0;
+                }
+                else if (B_COND == BC_PERIODIC)
+                {
+                    iplus = (i+1) % NX;
+                    iminus = (i-1) % NX;
+                    if (iminus < 0) iminus += NX;
+                    jplus = (j+1) % NY;
+                    jminus = (j-1) % NY;
+                    if (jminus < 0) jminus += NY;
+                }
+                
+                delta1 = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
+
+                x = phi_in[i][j];
+
+                /* evolve phi */
+                if (B_COND != BC_ABSORBING)
+                {
+                    phi_out[i][j] = x + intstep*(delta1 - SPEED*(phi_in[iplus][j] - phi_in[i][j]));
+                }
+                else        /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
+                {
+                    /* in the bulk */
+                    if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
+                    {
+                        phi_out[i][j] = x - intstep*delta2;
+                    }
+                     /* right border */
+                    else if (i==NX-1) 
+                    {
+                        phi_out[i][j] = x - intstep1*(x - phi_in[i-1][j]);
+                    }
+                    /* upper border */
+                    else if (j==NY-1) 
+                    {
+                        phi_out[i][j] = x - intstep1*(x - phi_in[i][j-1]);
+                    }
+                    /* left border */
+                    else if (i==0) 
+                    {
+                        phi_out[i][j] = x - intstep1*(x - phi_in[1][j]);
+                    }
+                   /* lower border */
+                    else if (j==0) 
+                    {
+                        phi_out[i][j] = x - intstep1*(x - phi_in[i][1]);
+                    }
+                }
+
+
+                if (FLOOR)
+                {
+                    if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
+                    if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
+                }
+            }
+        }
+    }
+    
+//     printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
+}
+
+void evolve_wave(phi, phi_tmp, xy_in)
+/* time step of field evolution */
+    double *phi[NX], *phi_tmp[NX]; short int *xy_in[NX];
+{
+    evolve_wave_half(phi, phi_tmp, xy_in);
+    evolve_wave_half(phi_tmp, phi, xy_in);
+}
+
+
+void old_evolve_wave(phi, xy_in)
 /* time step of field evolution */
     double *phi[NX]; short int *xy_in[NX];
 {
     int i, j, iplus, iminus, jplus, jminus;
-    double delta1, delta2, x, y, *newphi[NX];;
+    double delta1, delta2, x, y, *newphi[NX];
     
     for (i=0; i<NX; i++) newphi[i] = (double *)malloc(NY*sizeof(double));
 
@@ -552,9 +595,6 @@ void evolve_wave(phi, xy_in)
 //     printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
 }
 
-
-
-
 double compute_variance(phi, xy_in)
 /* compute the variance (total probability) of the field */
 double *phi[NX]; short int * xy_in[NX];
@@ -658,7 +698,6 @@ short int *xy_in[NX];
     sinj = sin(jangle);
     julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
     julia_y = 0.5*sinj*(1.0-cosj) + yshift;
-    /* need to decrease 0.05 for i > 2000 */
 //     julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
 //     julia_y = 0.5*sinj*(1.0-cosj);
     init_julia_set(phi, xy_in);
@@ -669,7 +708,7 @@ short int *xy_in[NX];
 void animation()
 {
     double time, scale, dx, var, jangle, cosj, sinj;
-    double *phi[NX];
+    double *phi[NX], *phi_tmp[NX];
     short int *xy_in[NX];
     int i, j, s;
 
@@ -677,6 +716,7 @@ void animation()
     for (i=0; i<NX; i++)
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
+        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
 
@@ -687,7 +727,7 @@ void animation()
 //     julia_x = 0.1; 
 //     julia_y = 0.6; 
     
-    set_Julia_parameters(0, phi, xy_in);
+//     set_Julia_parameters(0, phi, xy_in);
     
     printf("Integration step %.3lg\n", intstep);
 
@@ -710,7 +750,7 @@ void animation()
     
     draw_wave(phi, xy_in, 1.0, 0);
     draw_billiard();
-    print_Julia_parameters(i);
+//     print_Julia_parameters(i);
     
 //     print_level(MDEPTH);
 
@@ -734,17 +774,17 @@ void animation()
         
         draw_wave(phi, xy_in, scale, i);
         
-        for (j=0; j<NVID; j++) evolve_wave(phi, xy_in);
+        for (j=0; j<NVID; j++) evolve_wave(phi, phi_tmp, xy_in);
 
         draw_billiard();
         
 //         print_level(MDEPTH);
-        print_Julia_parameters(i);
+//         print_Julia_parameters(i);
 
 	glutSwapBuffers();
         
         /* modify Julia set */
-        set_Julia_parameters(i, phi, xy_in);
+//         set_Julia_parameters(i, phi, xy_in);
 
 	if (MOVIE)
         {
@@ -770,6 +810,7 @@ void animation()
     for (i=0; i<NX; i++)
     {
         free(phi[i]);
+        free(phi_tmp[i]);
     }
 
 }
diff --git a/particle_billiard.c b/particle_billiard.c
index c884e99..e0c9b17 100644
--- a/particle_billiard.c
+++ b/particle_billiard.c
@@ -47,23 +47,20 @@
 // #define YMIN -0.91
 // #define YMAX 1.115	/* y interval for 9/16 aspect ratio */
 
-/* Choice of the billiard table */
+/* Choice of the billiard table, see global_particles.c */
 
-#define B_DOMAIN 9      /* choice of domain shape */
+#define B_DOMAIN 14      /* choice of domain shape */
 
-#define D_RECTANGLE 0   /* rectangular domain */
-#define D_ELLIPSE 1     /* elliptical domain */
-#define D_STADIUM 2     /* stadium-shaped domain */
-#define D_SINAI 3       /* Sinai billiard */
-#define D_DIAMOND 4     /* diamond-shaped billiard */
-#define D_TRIANGLE 5    /* triangular billiard */
-#define D_ANNULUS 7     /* annulus */
-#define D_POLYGON 8     /* polygon */
-#define D_REULEAUX 9    /* Reuleaux and star shapes */
-#define D_FLOWER 10     /* Bunimovich flower */
-#define D_ALT_REU 11    /* alternating between star and Reuleaux */
+#define CIRCLE_PATTERN 0    /* pattern of circles */
 
-#define LAMBDA -3.346065215	/* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */
+#define NMAXCIRCLES 1000        /* total number of circles (must be at least NCX*NCY for square grid) */
+// #define NCX 10            /* number of circles in x direction */
+// #define NCY 15            /* number of circles in y direction */
+#define NCX 15            /* number of circles in x direction */
+#define NCY 20            /* number of circles in y direction */
+
+#define LAMBDA 0.75	/* parameter controlling shape of billiard */
+// #define LAMBDA -3.346065215	/* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */
 // #define LAMBDA 3.0	/* parameter controlling shape of billiard */
 // #define LAMBDA 0.6	/* parameter controlling shape of billiard */
 // #define LAMBDA 0.4175295	/* sin(20°)/sin(55°) for 9-star shape with 30° angles */
@@ -71,28 +68,30 @@
 // #define LAMBDA 3.75738973	/* sin(36°)/sin(9°) for 5-star shape with 90° angles */
 // #define LAMBDA -1.73205080756888	/* -sqrt(3) for Reuleaux triangle */
 // #define LAMBDA 1.73205080756888	/* sqrt(3) for triangle tiling plane */
-#define MU 0.1          /* second parameter controlling shape of billiard */
+#define MU 0.035          /* 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 2.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define NPOLY 8             /* number of sides of polygon */
+#define APOLY 0.25           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define DRAW_BILLIARD 1     /* set to 1 to draw billiard */
 #define DRAW_CONSTRUCTION_LINES 1   /* set to 1 to draw additional construction lines for billiard */
+#define PERIODIC_BC 0       /* set to 1 to enforce periodic boundary conditions when drawing particles */
 
 #define RESAMPLE 0      /* set to 1 if particles should be added when dispersion too large */
 #define DEBUG 0         /* draw trajectories, for debugging purposes */
 
 /* Simulation parameters */
 
-#define NPART 20000     /* number of particles */
+#define NPART 5000       /* number of particles */
 #define NPARTMAX 100000	/* maximal number of particles after resampling */
 #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 NSTEPS 6000     /* number of frames of movie */
+#define NSTEPS 3000     /* number of frames of movie */
 #define TIME 1000       /* time between movie frames, for fluidity of real-time simulation */ 
-// #define DPHI 0.000005    /* integration step */
-#define DPHI 0.00002    /* integration step */
+// #define DPHI 0.000002    /* integration step */
+// #define DPHI 0.00002    /* integration step */
+#define DPHI 0.000005    /* integration step */
 #define NVID 150         /* number of iterations between images displayed on screen */
 
 /* Decreasing TIME accelerates the animation and the movie                               */
@@ -103,12 +102,13 @@
 
 /* Colors and other graphical parameters */
 
-#define NCOLORS -10      /* number of colors */
-#define COLORSHIFT 200     /* hue of initial color */ 
+#define NCOLORS 16      /* 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.04      /* length of velocity vectors */
-#define BILLIARD_WIDTH 3    /* width of billiard */
+#define LENGTH 0.02      /* length of velocity vectors */
+#define BILLIARD_WIDTH 2    /* width of billiard */
 #define PARTICLE_WIDTH 2    /* width of particles */
 #define FRONT_WIDTH 3       /* width of wave front */
 
@@ -123,10 +123,8 @@
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1000      /* final sleeping time */
 
-#define PI 	3.141592654
-#define DPI 	6.283185307
-#define PID 	1.570796327
 
+#include "global_particles.c"
 #include "sub_part_billiard.c"
 
 
@@ -169,7 +167,8 @@ double *configs[NPARTMAX];
     double conf[2], pos[2];
   
     while (angle2 < angle1) angle2 += DPI;
-    dalpha = (angle2 - angle1)/((double)(NPART-1));
+    if (NPART > 1) dalpha = (angle2 - angle1)/((double)(NPART-1));
+    else dalpha = 0.0;
     for (i=0; i<NPART; i++) 
     {
         alpha = angle1 + dalpha*((double)i);  
@@ -253,8 +252,11 @@ double *configs[NPARTMAX];
         if (configs[i][2]<0.0) 
         {    
             vbilliard(configs[i]);
-            color[i]++;
-            if (color[i] >= NCOLORS) color[i] -= NCOLORS;
+            if (!RAINBOW_COLOR)
+            {
+                color[i]++;
+                if (color[i] >= NCOLORS) color[i] -= NCOLORS;
+            }
         }
 
         configs[i][2] += DPHI; 
@@ -280,36 +282,39 @@ double *configs[NPARTMAX];
             glEnd ();
         
             /* taking care of boundary conditions - only needed for periodic boundary conditions */
-            if (SCALING_FACTOR*x2 > XMAX)
+            if (PERIODIC_BC)
             {
-                glBegin(GL_LINE_STRIP);
-                glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
-                glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
-                glEnd ();
-            }
+                if (SCALING_FACTOR*x2 > XMAX)
+                {
+                    glBegin(GL_LINE_STRIP);
+                    glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
+                    glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
+                    glEnd ();
+                }
         
-            if (SCALING_FACTOR*x2 < XMIN)
-            {
-                glBegin(GL_LINE_STRIP);
-                glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
-                glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
-                glEnd ();
-            }
+                if (SCALING_FACTOR*x2 < XMIN)
+                {
+                    glBegin(GL_LINE_STRIP);
+                    glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
+                    glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
+                    glEnd ();
+                }
 
-            if (SCALING_FACTOR*y2 > YMAX)
-            {
-                glBegin(GL_LINE_STRIP);
-                glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
-                glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
-                glEnd ();
-            }
+                if (SCALING_FACTOR*y2 > YMAX)
+                {
+                    glBegin(GL_LINE_STRIP);
+                    glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
+                    glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
+                    glEnd ();
+                }
 
-            if (SCALING_FACTOR*y2 < YMIN)
-            {
-                glBegin(GL_LINE_STRIP);
-                glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
-                glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
-                glEnd ();
+                if (SCALING_FACTOR*y2 < YMIN)
+                {
+                    glBegin(GL_LINE_STRIP);
+                    glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
+                    glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
+                    glEnd ();
+                }
             }
         }
         
@@ -343,11 +348,14 @@ double *configs[NPARTMAX];
         {      
             if (configs[i][2]<0.0) 
             {    
+//                 printf("reflecting particle %i\n", i);
                 c = vbilliard(configs[i]);
 //                 if (c>=0) color[i]++;
-                color[i]++;
-                if (color[i] >= NCOLORS) color[i] -= NCOLORS;
-                
+                if (!RAINBOW_COLOR)
+                {
+                    color[i]++;
+                    if (color[i] >= NCOLORS) color[i] -= NCOLORS;
+                }                
             }
 
             configs[i][2] += DPHI; 
@@ -363,7 +371,74 @@ double *configs[NPARTMAX];
 }
 
 
+void init_circle_config()
+{
+    int i, j, n; 
+    double dx, dy;
+    
+    switch (CIRCLE_PATTERN) {
+        case (C_FOUR_CIRCLES):
+        {
+            ncircles = 4;
+            
+            circlex[0] = 1.0;
+            circley[0] = 0.0;
+            circlerad[0] = 0.8;
+            
+            circlex[1] = -1.0;
+            circley[1] = 0.0;
+            circlerad[1] = 0.8;
+            
+            circlex[2] = 0.0;
+            circley[2] = 0.8;
+            circlerad[2] = 0.4;
+            
+            circlex[3] = 0.0;
+            circley[3] = -0.8;
+            circlerad[3] = 0.4;
+            
+            for (i=0; i<4; i++) circleactive[i] = 1;
 
+            break;
+        }
+        case (C_SQUARE):
+        {
+            ncircles = NCX*NCY;
+            dy = (YMAX - YMIN)/((double)NCY);
+            for (i = 0; i < NCX; i++)
+                for (j = 0; j < NCY; j++)
+                {
+                    n = NCY*i + j;
+                    circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
+                    circley[n] = YMIN + ((double)j + 0.5)*dy;
+                    circlerad[n] = MU;
+                    circleactive[n] = 1;
+                }
+            break;
+        }
+        case (C_HEX):
+        {
+            ncircles = NCX*(NCY+1);
+            dy = (YMAX - YMIN)/((double)NCY);
+            dx = dy*0.5*sqrt(3.0);
+            for (i = 0; i < NCX; i++)
+                for (j = 0; j < NCY+1; j++)
+                {
+                    n = (NCY+1)*i + j;
+                    circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
+                    circley[n] = YMIN + ((double)j - 0.5)*dy;
+                    if ((i+NCX)%2 == 1) circley[n] += 0.5*dy;
+                    circlerad[n] = MU;
+                    circleactive[n] = 1;
+                }
+            break;
+        }
+        default: 
+        {
+            printf("Function init_circle_config not defined for this pattern \n");
+        }
+    }
+}
 
 
 void animation()
@@ -379,19 +454,24 @@ void animation()
     active = malloc(sizeof(int)*(NPARTMAX));
     for (i=0; i<NPARTMAX; i++)
         configs[i] = (double *)malloc(8*sizeof(double));
+    
+    /* init circle configuration if the domain is D_CIRCLES */
+    if (B_DOMAIN == D_CIRCLES) init_circle_config();
+    
       
     /* initialize system by putting particles in a given point with a range of velocities */
     r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY);
 
-//     init_drop_config(0.0, 0.0, -0.2, 0.2, configs);    
+//     init_drop_config(0.0, 0.0, 0.0, PI, configs);    
 //     init_drop_config(0.5, 0.5, -1.0, 1.0, configs);    
 //     init_sym_drop_config(-1.0, 0.5, -PID, PID, configs);
 //     init_drop_config(-0.999, 0.0, -alpha, alpha, configs);
 
 //  other possible initial conditions :
-//     init_line_config(-0.6, 0.2, -0.6, 0.7, 0.0, configs);
+//     init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs);
+    init_line_config(0.0, -1.0, -1.0, 1.0, 0.25*PID, configs);
 //     init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs);
-    init_line_config(0.0, -0.3, 0.0, 0.3, 0.0, configs);
+//     init_line_config(-1.5, 0.1, -0.1, 1.0, -0.5*PID, configs);
   
     blank();  
     glColor3f(0.0, 0.0, 0.0);
@@ -425,7 +505,14 @@ void animation()
         }
     }
   
-    sleep(SLEEP1);
+    if (RAINBOW_COLOR)      /* rainbow color scheme */
+        for (i=0; i<NPART; i++) 
+        {
+            color[i] = (i*NCOLORS)/NPART;
+            newcolor[i] = (i*NCOLORS)/NPART;  
+        }
+
+        sleep(SLEEP1);
   
     for (i=0; i<=NSTEPS; i++)
     {
diff --git a/schrodinger.c b/schrodinger.c
index edea4c0..77430c0 100644
--- a/schrodinger.c
+++ b/schrodinger.c
@@ -42,8 +42,10 @@
 #define WINWIDTH 	1280  /* window width */
 #define WINHEIGHT 	720   /* window height */
 
-#define NX 640          /* number of grid points on x axis */
-#define NY 360          /* number of grid points on y axis */
+#define NX 1280          /* number of grid points on x axis */
+#define NY 720          /* number of grid points on y axis */
+// #define NX 640          /* number of grid points on x axis */
+// #define NY 360          /* number of grid points on y axis */
 
 /* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
 /* but will multiply run time by 4                                 */
@@ -55,67 +57,45 @@
 
 #define JULIA_SCALE 1.0 /* scaling for Julia sets */
 
-/* Choice of the billiard table */
+/* Choice of the billiard table, see list in global_pdes.c  */
 
-#define B_DOMAIN 3      /* choice of domain shape */
+#define B_DOMAIN 9      /* choice of domain shape */
 
-#define D_RECTANGLE 0   /* rectangular domain */
-#define D_ELLIPSE 1     /* elliptical domain */
-#define D_STADIUM 2     /* stadium-shaped domain */
-#define D_SINAI 3       /* Sinai billiard */
-#define D_DIAMOND 4     /* diamond-shaped billiard */
-#define D_TRIANGLE 5    /* triangular billiard */
-#define D_FLAT 6        /* flat interface */
-#define D_ANNULUS 7     /* annulus */
-#define D_POLYGON 8     /* polygon */
-#define D_YOUNG 9       /* Young diffraction slits */
-#define D_GRATING 10    /* diffraction grating */
-#define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
+#define CIRCLE_PATTERN 0    /* pattern of circles, see list in global_pdes.c */
 
-#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
-#define D_JULIA_INT 16  /* interior of Julia set */ 
+#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 */
 
-/* Billiard tables for heat equation */
-
-#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
-#define D_MENGER_HEATED 22  /* Menger gasket with different temperatures */
-#define D_MENGER_H_OPEN 23  /* Menger gasket with different temperatures and larger domain */
-#define D_MANDELBROT 24     /* Mandelbrot set */
-#define D_JULIA 25          /* Julia set */
-#define D_MANDELBROT_CIRCLE 26     /* Mandelbrot set with circular conductor */
-
-#define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.3	    /* parameter controlling the dimensions of domain */
 #define MU 0.05	            /* 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 2            /* depth of computation of Menger gasket */
-#define MRATIO 5            /* ratio defining Menger gasket */
+#define MDEPTH 3            /* depth of computation of Menger gasket */
+#define MRATIO 3            /* ratio defining Menger gasket */
 #define MANDELLEVEL 1000      /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0     /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
+#define NGRIDX 15            /* number of grid point for grid of disks */
+#define NGRIDY 20           /* number of grid point for grid of disks */
 
 /* You can add more billiard tables by adapting the functions */
 /* xy_in_billiard and draw_billiard in sub_wave.c */
 
 /* Physical patameters of wave equation */
 
-#define DT 0.00000005
-// #define DT 0.00000002
+// #define DT 0.00000005
+#define DT 0.000000005
 // #define DT 0.000000005
 #define HBAR 1.0
 
-/* Boundary conditions */
+/* Boundary conditions, see list in global_pdes.c  */
 
 #define B_COND 1
 
-#define BC_DIRICHLET 0   /* Dirichlet boundary conditions */
-#define BC_PERIODIC 1    /* periodic boundary conditions */
-#define BC_ABSORBING 2   /* absorbing boundary conditions (beta version) */
-
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 4500      /* number of frames of movie */
-#define NVID 250         /* number of iterations between images displayed on screen */
+#define NSTEPS 200      /* number of frames of movie */
+#define NVID 1200         /* number of iterations between images displayed on screen */
 #define NSEG 100         /* number of segments of boundary */
 
 #define PAUSE 1000       /* number of frames after which to pause */
@@ -128,25 +108,17 @@
 #define VMAX 10.0       /* max value of wave amplitude */
 
 
-/* Plot type */
+/* Plot type, see list in global_pdes.c  */
 
-#define PLOT 0
+#define PLOT 10
 
-#define P_MODULE 0        /* plot module of wave function squared */
-#define P_PHASE 1         /* plot phase of wave function */
-#define P_REAL 2          /* plot real part */
-#define P_IMAGINARY 3     /* plot imaginary part */
 
-/* Color schemes */
+/* Color schemes, see list in global_pdes.c  */
 
 #define BLACK 1          /* black background */
 
 #define COLOR_SCHEME 1   /* choice of color scheme */
 
-#define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
-#define C_HUE 1          /* color scheme modifies hue */
-#define C_PHASE 2        /* color scheme shows phase */
-
 #define SCALE 1          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
@@ -158,14 +130,7 @@
 #define HUEMEAN 150.0    /* mean value of hue for color scheme C_HUE */
 #define HUEAMP -150.0      /* amplitude of variation of hue for color scheme C_HUE */
 
-/* Basic math */
-
-#define PI 	3.141592654
-#define DPI 	6.283185307
-#define PID 	1.570796327
-
-double julia_x = 0.0, julia_y = 0.0;    /* parameters for Julia sets */
-
+#include "global_pdes.c"
 #include "sub_wave.c"
 
 double courant2;  /* Courant parameter squared */
@@ -269,22 +234,14 @@ int time;
     glEnd ();
 }
 
-void evolve_wave(phi, psi, xy_in)
+void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
 /* time step of field evolution */
 /* phi is real part, psi is imaginary part */
-    double *phi[NX], *psi[NX]; 
-    short int *xy_in[NX];
+    double *phi_in[NX], *psi_in[NX], *phi_out[NX], *psi_out[NX]; short int *xy_in[NX];
 {
     int i, j, iplus, iminus, jplus, jminus;
-    double delta1, delta2, x, y, *newphi[NX], *newpsi[NX];
+    double delta1, delta2, x, y;
     
-    for (i=0; i<NX; i++) 
-    {
-        newphi[i] = (double *)malloc(NY*sizeof(double));
-        newpsi[i] = (double *)malloc(NY*sizeof(double));        
-    }
-    
-
     #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
     for (i=0; i<NX; i++){
         for (j=0; j<NY; j++){
@@ -307,81 +264,76 @@ void evolve_wave(phi, psi, xy_in)
                     if (jminus < 0) jminus += NY;
                 }
                 
-                delta1 = phi[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
-                delta2 = psi[iplus][j] + psi[iminus][j] + psi[i][jplus] + psi[i][jminus] - 4.0*psi[i][j];
+                delta1 = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
+                delta2 = psi_in[iplus][j] + psi_in[iminus][j] + psi_in[i][jplus] + psi_in[i][jminus] - 4.0*psi_in[i][j];
 
-                x = phi[i][j];
-		y = psi[i][j];
+                x = phi_in[i][j];
+		y = psi_in[i][j];
 
                 /* evolve phi and psi */
                 if (B_COND != BC_ABSORBING)
                 {
-                    newphi[i][j] = x - intstep*delta2;
-                    newpsi[i][j] = y + intstep*delta1;
+                    phi_out[i][j] = x - intstep*delta2;
+                    psi_out[i][j] = y + intstep*delta1;
                 }
                 else        /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
                 {
                     /* in the bulk */
                     if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
                     {
-                        newphi[i][j] = x - intstep*delta2;
-                        newpsi[i][j] = y + intstep*delta1;
+                        phi_out[i][j] = x - intstep*delta2;
+                        psi_out[i][j] = y + intstep*delta1;
                     }
                      /* right border */
                     else if (i==NX-1) 
                     {
-                        newphi[i][j] = x - intstep1*(y - psi[i-1][j]);
-                        newpsi[i][j] = y + intstep1*(x - phi[i-1][j]);
+                        phi_out[i][j] = x - intstep1*(y - psi_in[i-1][j]);
+                        psi_out[i][j] = y + intstep1*(x - phi_in[i-1][j]);
                     }
                     /* upper border */
                     else if (j==NY-1) 
                     {
-                        newphi[i][j] = x - intstep1*(y - psi[i][j-1]);
-                        newpsi[i][j] = y + intstep1*(x - phi[i][j-1]);
+                        phi_out[i][j] = x - intstep1*(y - psi_in[i][j-1]);
+                        psi_out[i][j] = y + intstep1*(x - phi_in[i][j-1]);
                     }
                     /* left border */
                     else if (i==0) 
                     {
-                        newphi[i][j] = x - intstep1*(y - psi[1][j]);
-                        newpsi[i][j] = y + intstep1*(x - phi[1][j]);
+                        phi_out[i][j] = x - intstep1*(y - psi_in[1][j]);
+                        psi_out[i][j] = y + intstep1*(x - phi_in[1][j]);
                     }
                    /* lower border */
                     else if (j==0) 
                     {
-                        newphi[i][j] = x - intstep1*(y - psi[i][1]);
-                        newpsi[i][j] = y + intstep1*(x - phi[i][1]);
+                        phi_out[i][j] = x - intstep1*(y - psi_in[i][1]);
+                        psi_out[i][j] = y + intstep1*(x - phi_in[i][1]);
                     }
                 }
 
 
                 if (FLOOR)
                 {
-                    if (newphi[i][j] > VMAX) newphi[i][j] = VMAX;
-                    if (newphi[i][j] < -VMAX) newphi[i][j] = -VMAX;
-                    if (newpsi[i][j] > VMAX) newpsi[i][j] = VMAX;
-                    if (newpsi[i][j] < -VMAX) newpsi[i][j] = -VMAX;
+                    if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
+                    if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
+                    if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
+                    if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
                 }
             }
         }
     }
     
-    
-    for (i=0; i<NX; i++){
-        for (j=0; j<NY; j++){
-            if (xy_in[i][j] == 1) phi[i][j] = newphi[i][j];
-            if (xy_in[i][j] == 1) psi[i][j] = newpsi[i][j];
-        }
-    }
-    
-    for (i=0; i<NX; i++)
-    {
-        free(newphi[i]);
-        free(newpsi[i]);
-    }
-
 //     printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
 }
 
+void evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in)
+/* time step of field evolution */
+/* phi is value of field at time t, psi at time t-1 */
+    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX]; short int *xy_in[NX];
+{
+    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
+    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
+}
+
 
 double compute_variance(phi, psi, xy_in)
 /* compute the variance (total probability) of the field */
@@ -428,7 +380,7 @@ short int * xy_in[NX];
 void animation()
 {
     double time, scale, dx, var;
-    double *phi[NX], *psi[NX];
+    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
     short int *xy_in[NX];
     int i, j, s;
 
@@ -437,6 +389,8 @@ void animation()
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
         psi[i] = (double *)malloc(NY*sizeof(double));
+        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
+        psi_tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
 
@@ -447,7 +401,7 @@ void animation()
     printf("Integration step %.3lg\n", intstep);
 
     /* initialize wave wave function */
-    init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.2, phi, psi, xy_in);
+    init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.25, phi, psi, xy_in);
 //     init_coherent_state(0.0, 0.0, 0.0, 5.0, 0.03, phi, psi, xy_in);
 //     init_coherent_state(-0.5, 0.0, 1.0, 1.0, 0.05, phi, psi, xy_in);
     
@@ -488,7 +442,7 @@ void animation()
         
 //         printf("Wave drawn\n");
         
-        for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
+        for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
         
         draw_billiard();
         
@@ -510,15 +464,18 @@ void animation()
 
     }
 
-    if (MOVIE)
-    {
-        for (i=0; i<20; i++) save_frame();
-        s = system("mv wave*.tif tif_schrod/");
-    }
+//     if (MOVIE)
+//     {
+//         for (i=0; i<20; i++) save_frame();
+//         s = system("mv wave*.tif tif_schrod/");
+//     }
     for (i=0; i<NX; i++)
     {
         free(phi[i]);
         free(psi[i]);
+        free(phi_tmp[i]);
+        free(psi_tmp[i]);
+        free(xy_in[i]);
     }
 
 }
diff --git a/sub_wave.c b/sub_wave.c
index 187b31b..11f8e37 100644
--- a/sub_wave.c
+++ b/sub_wave.c
@@ -315,13 +315,129 @@ double x1, y1, x2, y2;
 }
 
 
+void draw_rotated_rectangle(x1, y1, x2, y2)
+double x1, y1, x2, y2;
+{
+    double pos[2];
+    double xa, ya, xb, yb, xc, yc;
+    
+    xa = 0.5*(x1 - y2);
+    xb = 0.5*(x2 - y1);
+    xc = 0.5*(x1 - y1);
+    ya = 0.5*(x1 + y1);
+    yb = 0.5*(x2 + y2);
+    yc = 0.5*(x2 + y1);
+    
+    glBegin(GL_LINE_LOOP);
+    xy_to_pos(xc, ya, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(xb, yc, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(xc, yb, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(xa, yc, pos);
+    glVertex2d(pos[0], pos[1]);
+    glEnd();    
+}
+
+void init_circle_config()
+/* initialise the arrays circlex, circley, circlerad and circleactive */
+/* for billiard shape D_CIRCLES */
+{
+    int i, j, n; 
+    double dx, dy, p;
+    
+    switch (CIRCLE_PATTERN) {
+        case (C_SQUARE):
+        {
+            ncircles = NGRIDX*NGRIDY;
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    n = NGRIDY*i + j;
+                    circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
+                    circley[n] = YMIN + ((double)j + 0.5)*dy;
+                    circlerad[n] = MU;
+                    circleactive[n] = 1;
+                }
+            break;
+        }
+        case (C_HEX):
+        {
+            ncircles = NGRIDX*(NGRIDY+1);
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            dx = dy*0.5*sqrt(3.0);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY+1; j++)
+                {
+                    n = (NGRIDY+1)*i + j;
+                    circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
+                    circley[n] = YMIN + ((double)j - 0.5)*dy;
+                    if ((i+NGRIDX)%2 == 1) circley[n] += 0.5*dy;
+                    circlerad[n] = MU;
+                    circleactive[n] = 1;
+                }
+            break;
+        }
+        case (C_RAND_DISPLACED):
+        {
+            ncircles = NGRIDX*NGRIDY;
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    n = NGRIDY*i + j;
+                    circlex[n] = ((double)(i-NGRIDX/2) + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy;
+                    circley[n] = YMIN + ((double)j + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy;
+                    circlerad[n] = MU*(1.0 + 0.35*(double)rand()/RAND_MAX);
+                    circleactive[n] = 1;
+                }
+            break;
+        }
+        case (C_RAND_PERCOL):
+        {
+            ncircles = NGRIDX*NGRIDY;
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    n = NGRIDY*i + j;
+                    circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
+                    circley[n] = YMIN + ((double)j + 0.5)*dy;
+                    circlerad[n] = MU;
+                    p = (double)rand()/RAND_MAX;
+                    if (p < P_PERCOL) circleactive[n] = 1;
+                    else circleactive[n] = 0;
+                }
+            break;
+        }
+        case (C_RAND_POISSON):
+        {
+            ncircles = NPOISSON;
+            for (n = 0; n < NPOISSON; n++)
+            {
+                circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
+                circley[n] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN;
+                circlerad[n] = MU;
+                circleactive[n] = 1;
+            }
+            break;
+        }
+        default: 
+        {
+            printf("Function init_circle_config not defined for this pattern \n");
+        }
+    }
+}
+
 
 int xy_in_billiard(x, y)
 /* returns 1 if (x,y) represents a point in the billiard */
 double x, y;
 {
-    double l2, r2, r2mu, omega, c, angle, z, x1, y1, u, v, u1, v1;
-    int i, k, k1, k2, condition;
+    double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy;
+    int i, j, k, k1, k2, condition;
 
     switch (B_DOMAIN) {
         case D_RECTANGLE:
@@ -416,6 +532,44 @@ double x, y;
         {
             return(((x-1.0)*(x-1.0) + y*y < LAMBDA*LAMBDA)||((x+1.0)*(x+1.0) + y*y < LAMBDA*LAMBDA)||((vabs(x) < 1.0)&&(vabs(y) < MU)));
         }
+        case D_DISK_GRID:
+        {
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            for (i = -NGRIDX/2; i < NGRIDX/2; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    x1 = ((double)i + 0.5)*dy;
+                    y1 = YMIN + ((double)j + 0.5)*dy;
+                    if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0); 
+                }
+            return(1);
+        }
+        case D_DISK_HEX:
+        {
+            dy = (YMAX - YMIN)/((double)NGRIDY);
+            dx = dy*0.5*sqrt(3.0);
+            for (i = -NGRIDX/2; i < NGRIDX/2; i++)
+                for (j = -1; j < NGRIDY; j++)
+                {
+                    x1 = ((double)i + 0.5)*dy;
+                    y1 = YMIN + ((double)j + 0.5)*dy;
+                    if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy;
+                    if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0); 
+                }
+            return(1);
+        }
+        case D_CIRCLES:
+        {
+            for (i = 0; i < ncircles; i++)
+                if (circleactive[i]) 
+                {
+                    x1 = circlex[i];
+                    y1 = circley[i];
+                    r2 = circlerad[i]*circlerad[i];
+                    if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0); 
+                }
+            return(1);
+        }
         case D_MENGER:       
         {
             x1 = 0.5*(x+1.0);
@@ -443,6 +597,23 @@ double x, y;
             if (u*u + v*v < MANDELLIMIT) return(1);
             else return(0);
         }
+        case D_MENGER_ROTATED:       
+        {
+            x2 = 1.0*(x + y);
+            y2 = 1.0*(x - y);
+            if ((vabs(x2) < 1.0)&&(vabs(y2) < 1.0))
+            {
+                x1 = 0.5*(x2 + 1.0);
+                y1 = 0.5*(y2 + 1.0);
+                for (k=0; k<MDEPTH; k++)
+                {
+                    x1 = x1*(double)MRATIO;
+                    y1 = y1*(double)MRATIO;
+                    if ((vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(0);
+                }                
+            }
+            return(1);
+        }
         case D_ANNULUS_HEATED:      /* returns 2 if in inner circle */ 
         {
             l2 = LAMBDA*LAMBDA;
@@ -558,8 +729,8 @@ int i, j;
 
 void draw_billiard()      /* draws the billiard boundary */
 {
-    double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
-    int i, j, k1, k2, mr2;
+    double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
+    int i, j, k, k1, k2, mr2;
 
     if (BLACK) glColor3f(1.0, 1.0, 1.0);
     else glColor3f(0.0, 0.0, 0.0);
@@ -845,6 +1016,71 @@ void draw_billiard()      /* draws the billiard boundary */
             glEnd ();
             break;
         }
+        case D_DISK_GRID:
+        {
+            glLineWidth(2);
+            for (i = -NGRIDX/2; i < NGRIDX/2; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    dy = (YMAX - YMIN)/((double)NGRIDY);
+                    dx = dy*0.5*sqrt(3.0);
+                    x1 = ((double)i + 0.5)*dy;
+                    y1 = YMIN + ((double)j + 0.5)*dy;
+                    glBegin(GL_LINE_LOOP);
+                    for (k=0; k<=NSEG; k++)
+                    {
+                        phi = (double)k*DPI/(double)NSEG;
+                        x = x1 + MU*cos(phi);
+                        y = y1 + MU*sin(phi);
+                        xy_to_pos(x, y, pos);
+                        glVertex2d(pos[0], pos[1]);
+                    }
+                    glEnd ();
+                }
+            break;
+        }
+        case D_DISK_HEX:
+        {
+            glLineWidth(2);
+            for (i = -NGRIDX/2; i < NGRIDX/2; i++)
+                for (j = -1; j < NGRIDY; j++)
+                {
+                    dy = (YMAX - YMIN)/((double)NGRIDY);
+                    x1 = ((double)i + 0.5)*dy;
+                    y1 = YMIN + ((double)j + 0.5)*dy;
+                    if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy;
+                    glBegin(GL_LINE_LOOP);
+                    for (k=0; k<=NSEG; k++)
+                    {
+                        phi = (double)k*DPI/(double)NSEG;
+                        x = x1 + MU*cos(phi);
+                        y = y1 + MU*sin(phi);
+                        xy_to_pos(x, y, pos);
+                        glVertex2d(pos[0], pos[1]);
+                    }
+                    glEnd ();
+                }
+            break;
+        }
+        case (D_CIRCLES):
+        {
+            glLineWidth(2);
+            for (i = 0; i < ncircles; i++) 
+                if (circleactive[i]) 
+                {
+                    glBegin(GL_LINE_LOOP);
+                    for (k=0; k<=NSEG; k++)
+                    {
+                        phi = (double)k*DPI/(double)NSEG;
+                        x = circlex[i] + circlerad[i]*cos(phi);
+                        y = circley[i] + circlerad[i]*sin(phi);
+                        xy_to_pos(x, y, pos);
+                        glVertex2d(pos[0], pos[1]);
+                    }
+                    glEnd ();
+                }
+            break;
+        }
         case (D_MENGER):
         {
             glLineWidth(3);
@@ -893,11 +1129,59 @@ void draw_billiard()      /* draws the billiard boundary */
             
             break;
         }        
-         case (D_JULIA_INT):
+        case (D_JULIA_INT):
         {
             /* Do nothing */
             break;
         }
+        case (D_MENGER_ROTATED):
+        {
+            glLineWidth(3);
+//             draw_rectangle(XMIN, -1.0, XMAX, 1.0);
+
+            /* level 1 */
+            if (MDEPTH > 0)
+            {
+                glLineWidth(2);
+                x = 1.0/((double)MRATIO);
+                draw_rotated_rectangle(x, x, -x, -x);
+            }
+            
+            /* level 2 */
+            if (MDEPTH > 1)
+            {
+                glLineWidth(1);
+                mr2 = MRATIO*MRATIO;
+                l = 2.0/((double)mr2);
+                
+                for (i=0; i<MRATIO; i++)
+                    for (j=0; j<MRATIO; j++)
+                        if ((i!=MRATIO/2)||(j!=MRATIO/2))
+                        {
+                            x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)MRATIO);
+                            y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)MRATIO);
+                            draw_rotated_rectangle(x, y, x+l, y+l);
+                        }
+            }
+            
+            /* level 3 */
+            if (MDEPTH > 2)
+            {
+                glLineWidth(1);
+                l = 2.0/((double)(mr2*MRATIO));
+                
+                for (i=0; i<mr2; i++)
+                    for (j=0; j<mr2; j++)
+                        if ( (((i%MRATIO!=MRATIO/2))||(j%MRATIO!=MRATIO/2)) && (((i/MRATIO!=MRATIO/2))||(j/MRATIO!=MRATIO/2)) )
+                        {
+                            x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)mr2);
+                            y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)mr2);
+                            draw_rotated_rectangle(x, y, x+l, y+l);
+                        }
+            }
+            
+            break;
+        }        
         case (D_ANNULUS_HEATED):
         {
             glBegin(GL_LINE_LOOP);
diff --git a/wave_billiard.c b/wave_billiard.c
index 6ddfea5..4741448 100644
--- a/wave_billiard.c
+++ b/wave_billiard.c
@@ -58,130 +58,105 @@
 #define XMAX 2.0	/* x interval */
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
+// #define XMIN -1.8
+// #define XMAX 1.8	/* x interval */
+// #define YMIN -1.0125
+// #define YMAX 1.0125	/* y interval for 9/16 aspect ratio */
 
-#define JULIA_SCALE 1.1 /* scaling for Julia sets */
+#define JULIA_SCALE 1.0 /* scaling for Julia sets */
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 16      /* choice of domain shape */
+#define B_DOMAIN 3      /* choice of domain shape, see list in global_pdes.c */
 
-#define D_RECTANGLE 0   /* rectangular domain */
-#define D_ELLIPSE 1     /* elliptical domain */
-#define D_STADIUM 2     /* stadium-shaped domain */
-#define D_SINAI 3       /* Sinai billiard */
-#define D_DIAMOND 4     /* diamond-shaped billiard */
-#define D_TRIANGLE 5    /* triangular billiard */
-#define D_FLAT 6        /* flat interface */
-#define D_ANNULUS 7     /* annulus */
-#define D_POLYGON 8     /* polygon */
-#define D_YOUNG 9       /* Young diffraction slits */
-#define D_GRATING 10    /* diffraction grating */
-#define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
+#define CIRCLE_PATTERN 4    /* pattern of circles, see list in global_pdes.c */
 
-#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
-#define D_JULIA_INT 16  /* interior of Julia set */ 
+#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 */
 
-/* Billiard tables for heat equation */
-
-#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
-#define D_MENGER_HEATED 22  /* Menger gasket with different temperatures */
-#define D_MENGER_H_OPEN 23  /* Menger gasket with different temperatures and larger domain */
-#define D_MANDELBROT 24     /* Mandelbrot set */
-#define D_JULIA 25          /* Julia set */
-#define D_MANDELBROT_CIRCLE 26     /* Mandelbrot set with circular conductor */
-
-#define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
-#define MU 0.05	            /* parameter controlling the dimensions of domain */
-#define NPOLY 8             /* number of sides of polygon */
+#define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
+#define MU 0.025	    /* parameter controlling the dimensions of domain */
+#define NPOLY 3             /* number of sides of polygon */
 #define APOLY 1.0           /* 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 */
 #define MANDELLEVEL 1000      /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0     /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
+#define NGRIDX 15            /* number of grid point for grid of disks */
+#define NGRIDY 20           /* number of grid point for grid of disks */
 
 /* You can add more billiard tables by adapting the functions */
 /* xy_in_billiard and draw_billiard below */
 
 /* Physical parameters of wave equation */
 
+#define TWOSPEEDS 1         /* set to 1 to replace hardcore boundary by medium with different speed */
+
 #define OMEGA 0.9           /* frequency of periodic excitation */
 #define COURANT 0.01       /* Courant number */
+#define COURANTB 0.003      /* Courant number in medium B */
 #define GAMMA 0.0      /* damping factor in wave equation */
+#define GAMMA_BC 1.0e-4      /* damping factor on boundary */
 // #define GAMMA 5.0e-10      /* damping factor in wave equation */
 #define KAPPA 0.0       /* "elasticity" term enforcing oscillations */
-// #define KAPPA 5.0e-6       /* "elasticity" term enforcing oscillations */
-// #define KAPPA 5.0e-9       /* "elasticity" term enforcing oscillations */
-// #define KAPPA 5.0e-8       /* "elasticity" term enforcing oscillations */
+#define KAPPA_BC 5.0e-4       /* "elasticity" term on absorbing boundary */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
-/* Boundary conditions */
+/* Boundary conditions, see list in global_pdes.c  */
 
-#define B_COND 2
-
-#define BC_DIRICHLET 0   /* Dirichlet boundary conditions */
-#define BC_PERIODIC 1    /* periodic boundary conditions */
-#define BC_ABSORBING 2   /* absorbing boundary conditions (beta version) */
-
-
-/* For debugging purposes only */
-#define FLOOR 0         /* set to 1 to limit wave amplitude to VMAX */
-#define VMAX 10.0       /* max value of wave amplitude */
+#define B_COND 3
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 4000      /* number of frames of movie */
-#define NVID 25          /* number of iterations between images displayed on screen */
+#define NSTEPS 7000      /* number of frames of movie */
+#define NVID 50          /* number of iterations between images displayed on screen */
 #define NSEG 100         /* number of segments of boundary */
+#define INITIAL_TIME 200    /* time after which to start saving frames */
 
 #define PAUSE 1000         /* number of frames after which to pause */
 #define PSLEEP 1         /* sleep time during pause */
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1   /* final sleeping time */
 
+/* Plot type, see list in global_pdes.c  */
+
+#define PLOT 1
+
 /* Color schemes */
 
 #define BLACK 1          /* background */
 
-#define COLOR_SCHEME 1   /* choice of color scheme */
-
-#define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
-#define C_HUE 1          /* color scheme modifies hue */
+#define COLOR_SCHEME 1   /* choice of color scheme, see list in global_pdes.c  */
 
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
+// #define SLOPE 0.5       /* sensitivity of color on wave amplitude */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
+#define E_SCALE 750.0     /* scaling factor for energy representation */
 
 #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 */
 #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 230.0    /* mean value of hue for color scheme C_HUE */
-#define HUEAMP 50.0      /* amplitude of variation of hue for color scheme C_HUE */
+#define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
+#define HUEAMP -220.0      /* amplitude of variation of hue for color scheme C_HUE */
 // #define HUEMEAN 320.0    /* mean value of hue for color scheme C_HUE */
 // #define HUEAMP 100.0      /* amplitude of variation of hue for color scheme C_HUE */
 
-/* Basic math */
-
-#define PI 	3.141592654
-#define DPI 	6.283185307
-#define PID 	1.570796327
-
-double julia_x = -0.5, julia_y = 0.5;    /* parameters for Julia sets */
-// double julia_x = 0.33267, julia_y = 0.06395;    /* parameters for Julia sets */
-// double julia_x = 0.37468, julia_y = 0.21115;    /* parameters for Julia sets */
 
+#include "global_pdes.c"
 #include "sub_wave.c"
 
-double courant2;  /* Courant parameter squared */
+double courant2, courantb2;  /* Courant parameters squared */
 
 void init_wave(x, y, phi, psi, xy_in)
 /* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */
     double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX];
-
+ 
 {
     int i, j;
     double xy[2], dist2;
@@ -192,11 +167,36 @@ void init_wave(x, y, phi, psi, xy_in)
             ij_to_xy(i, j, xy);
             dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
 	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
-	    phi[i][j] = 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01);
+            
+	    if ((xy_in[i][j])||(TWOSPEEDS)) phi[i][j] = 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01);
+            else phi[i][j] = 0.0;
             psi[i][j] = 0.0;
         }
 }
 
+void init_planar_wave(x, y, phi, psi, xy_in)
+/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */
+/* beta version, works for vertical planar wave only so far */
+    double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX];
+ 
+{
+    int i, j;
+    double xy[2], dist2;
+
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+            dist2 = (xy[0]-x)*(xy[0]-x);
+	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
+            
+	    if ((xy_in[i][j])||(TWOSPEEDS)) phi[i][j] = 0.01*exp(-dist2/0.0005)*cos(-sqrt(dist2)/0.01);
+            else phi[i][j] = 0.0;
+            psi[i][j] = 0.0;
+        }
+}
+
+
 void init_wave_flat(phi, psi, xy_in)
 /* initialise flat field - phi is wave height, psi is phi at time t-1 */
     double *phi[NX], *psi[NX]; short int * xy_in[NX];
@@ -232,7 +232,7 @@ double factor, x, y, *phi[NX], *psi[NX];
 }
 
 void oscillate_linear_wave(amplitude, t, x1, y1, x2, y2, phi, psi)
-/* oscillating boundary condition at (x,y) */
+/* oscillating boundary condition at (x,y), beta version */
 double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
 {
     int i, j, ij1[2], ij2[2], imin, imax, jmin, jmax, d = 5;
@@ -265,6 +265,30 @@ double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
 /* animation part    */
 /*********************/
 
+double compute_energy(phi, psi, xy_in, i, j)
+double *phi[NX], *psi[NX];
+short int *xy_in[NX];
+int i, j;
+{
+    double velocity, energy, gradientx2, gradienty2;
+    int iplus, iminus, jplus, jminus;
+    
+    velocity = (phi[i][j] - psi[i][j]);
+                    
+    iplus = (i+1);   if (iplus == NX) iplus = NX-1;
+    iminus = (i-1);  if (iminus == -1) iminus = 0;
+    jplus = (j+1);   if (jplus == NY) jplus = NY-1;
+    jminus = (j-1);  if (jminus == -1) jminus = 0;
+                        
+    gradientx2 = (phi[iplus][j]-phi[i][j])*(phi[iplus][j]-phi[i][j]) 
+        + (phi[i][j] - phi[iminus][j])*(phi[i][j] - phi[iminus][j]);
+    gradienty2 = (phi[i][jplus]-phi[i][j])*(phi[i][jplus]-phi[i][j]) 
+        + (phi[i][j] - phi[i][jminus])*(phi[i][j] - phi[i][jminus]);
+    if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANT*COURANT*(gradientx2+gradienty2)));
+    else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANTB*COURANTB*(gradientx2+gradienty2)));
+    else return(0);
+}
+
 
 void draw_wave(phi, psi, xy_in, scale, time)
 /* draw the field */
@@ -272,17 +296,28 @@ double *phi[NX], *psi[NX], scale;
 short int *xy_in[NX];
 int time;
 {
-    int i, j;
-    double rgb[3], xy[2], x1, y1, x2, y2;
+    int i, j, iplus, iminus, jplus, jminus;
+    double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2;
+    static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
 
     glBegin(GL_QUADS);
+    
+//     printf("dtinverse = %.5lg\n", dtinverse);
 
     for (i=0; i<NX; i++)
         for (j=0; j<NY; j++)
         {
-            if (xy_in[i][j])
+            if ((TWOSPEEDS)||(xy_in[i][j]))
             {
-                color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
+                if (PLOT == P_AMPLITUDE)
+                    color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
+                else if (PLOT == P_ENERGY)
+                    color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
+                else if (PLOT == P_MIXED)
+                {
+                    if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
+                    else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
+                }
                 glColor3f(rgb[0], rgb[1], rgb[2]);
 
                 glVertex2i(i, j);
@@ -303,14 +338,25 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
     int i, j, iplus, iminus, jplus, jminus;
     double delta, x, y, c, cc;
     
-    c = COURANT;
-    cc = courant2;
+//     c = COURANT;
+//     cc = courant2;
 
     #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
     for (i=0; i<NX; i++){
         for (j=0; j<NY; j++){
-            if (xy_in[i][j]){
-                /* discretized Laplacian */
+            if (xy_in[i][j])
+            {
+                c = COURANT;
+                cc = courant2;
+            }
+            else if (TWOSPEEDS)
+            {
+                c = COURANTB;
+                cc = courantb2;                    
+            }
+
+            if ((TWOSPEEDS)||(xy_in[i][j])){
+                /* discretized Laplacian for various boundary conditions */
                 if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
                 {
                     iplus = (i+1);   if (iplus == NX) iplus = NX-1;
@@ -327,33 +373,55 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
                     jminus = (j-1) % NY;
                     if (jminus < 0) jminus += NY;
                 }
+                else if (B_COND == BC_VPER_HABS)
+                {
+                    iplus = (i+1);   if (iplus == NX) iplus = NX-1;
+                    iminus = (i-1);  if (iminus == -1) iminus = 0;
+                    jplus = (j+1) % NY;
+                    jminus = (j-1) % NY;
+                    if (jminus < 0) jminus += NY;
+                }
                 delta = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
 
                 x = phi_in[i][j];
 		y = psi_in[i][j];
 
                 /* evolve phi */
-                if (B_COND != BC_ABSORBING) phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
-                else
+                if ((B_COND == BC_PERIODIC)||(B_COND == BC_DIRICHLET)) 
+                    phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
+                else if (B_COND == BC_ABSORBING)
                 {
                     if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
-                        phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
+                        phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA_BC*(x-y);
                 
                     /* upper border */
                     else if (j==NY-1) 
-                        phi_out[i][j] = x - c*(x - phi_in[i][NY-2]) - KAPPA*x - GAMMA*(x-y);
+                        phi_out[i][j] = x - c*(x - phi_in[i][NY-2]) - KAPPA_BC*x - GAMMA_BC*(x-y);
                     
                     /* lower border */
                     else if (j==0) 
-                        phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA*x - GAMMA*(x-y);
+                        phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA_BC*x - GAMMA_BC*(x-y);
                 
                     /* right border */
                     if (i==NX-1) 
-                        phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA*x - GAMMA*(x-y);
+                        phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
                     
                     /* left border */
                     else if (i==0) 
-                        phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA*x - GAMMA*(x-y);
+                        phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
+                }
+                else if (B_COND == BC_VPER_HABS)
+                {
+                    if ((i>0)&&(i<NX-1))
+                        phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
+                
+                    /* right border */
+                    else if (i==NX-1) 
+                        phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
+                    
+                    /* left border */
+                    else if (i==0) 
+                        phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
                 }
                 psi_out[i][j] = x;
 
@@ -380,51 +448,7 @@ void evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in)
     evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
 }
 
-void old_evolve_wave(phi, psi, xy_in)
-/* time step of field evolution */
-/* phi is value of field at time t, psi at time t-1 */
-    double *phi[NX], *psi[NX]; short int *xy_in[NX];
-{
-    int i, j, iplus, iminus, jplus, jminus;
-    double delta, x, y;
 
-    #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
-    for (i=0; i<NX; i++){
-        for (j=0; j<NY; j++){
-            if (xy_in[i][j]){
-                /* discretized Laplacian */
-		iplus = (i+1) % NX;
-		iminus = (i-1) % NX;
-		if (iminus < 0) iminus += NX;
-                jplus = (j+1) % NY;
-                jminus = (j-1) % NY;
-                if (jminus < 0) jminus += NY;
-                delta = phi[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
-
-                x = phi[i][j];
-		y = psi[i][j];
-
-                /* evolve phi */
-                phi[i][j] = -y + 2*x + courant2*delta - KAPPA*x - GAMMA*(x-y);
-
-                /* Old versions of the simulation used this: */
-//                 phi[i][j] = (-psi[i][j] + 2*phi[i][j] + courant2*delta)*damping;
-//                 where damping = 1.0 - 0.0001;
-
-                psi[i][j] = x;
-
-                if (FLOOR)
-                {
-                    if (phi[i][j] > VMAX) phi[i][j] = VMAX;
-                    if (phi[i][j] < -VMAX) phi[i][j] = -VMAX;
-                    if (psi[i][j] > VMAX) psi[i][j] = VMAX;
-                    if (psi[i][j] < -VMAX) psi[i][j] = -VMAX;
-                }
-            }
-        }
-    }
-//     printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
-}
 
 
 double compute_variance(phi, psi, xy_in)
@@ -448,6 +472,7 @@ double compute_variance(phi, psi, xy_in)
 }
 
 
+
 void animation()
 {
     double time, scale;
@@ -464,12 +489,18 @@ void animation()
         psi_tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
+    
+    /* initialise positions and radii of circles */
+    if (B_DOMAIN == D_CIRCLES) init_circle_config();
 
     courant2 = COURANT*COURANT;
+    courantb2 = COURANTB*COURANTB;
 
     /* initialize wave with a drop at one point, zero elsewhere */
-//     init_wave_flat(phi, psi, xy_in);
-    init_wave(0.0, 0.0, phi, psi, xy_in);
+    init_planar_wave(XMIN + 0.01, 0.0, phi, psi, xy_in);
+//     init_planar_wave(XMIN + 1.0, 0.0, phi, psi, xy_in);
+//     init_wave(-1.5, 0.0, phi, psi, xy_in);
+//     init_wave(0.0, 0.0, phi, psi, xy_in);
 
     /* add a drop at another point */
 //     add_drop_to_wave(1.0, 0.7, 0.0, phi, psi);
@@ -487,7 +518,7 @@ void animation()
 
     sleep(SLEEP1);
 
-    for (i=0; i<=NSTEPS; i++)
+    for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
     {
 	//printf("%d\n",i);
         /* compute the variance of the field to adjust color scheme */
@@ -499,16 +530,12 @@ void animation()
         }
         else scale = 1.0;
 
-//         /* TO BE ADAPTED */
-//         scale = 1.0;
-
         draw_wave(phi, psi, xy_in, scale, i);
         for (j=0; j<NVID; j++) 
         {
             evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
-//         for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
         
         draw_billiard();
 
@@ -517,7 +544,8 @@ void animation()
 
 	if (MOVIE)
         {
-            save_frame();
+            if (i >= INITIAL_TIME) save_frame();
+            else printf("Initial phase time %i of %i\n", i, INITIAL_TIME);
 
             /* it seems that saving too many files too fast can cause trouble with the file system */
             /* so this is to make a pause from time to time - parameter PAUSE may need adjusting   */