diff --git a/drop_billiard.c b/drop_billiard.c
index ab951b0..a45d841 100644
--- a/drop_billiard.c
+++ b/drop_billiard.c
@@ -42,6 +42,8 @@
 // #define YMIN -1.125
 // #define YMAX 1.125	/* 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 */
 
 #define B_DOMAIN 9      /* choice of domain shape */
@@ -55,6 +57,8 @@
 #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 LAMBDA 1.0	/* parameter controlling shape of billiard */
 #define LAMBDA 1.124950941	/* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */
@@ -66,6 +70,8 @@
 #define FOCI 1          /* set to 1 to draw focal points of ellipse */
 #define NPOLY 5             /* number of sides of polygon */
 #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 RESAMPLE 0      /* set to 1 if particles should be added when dispersion too large */
 
diff --git a/heat.c b/heat.c
new file mode 100644
index 0000000..375af8b
--- /dev/null
+++ b/heat.c
@@ -0,0 +1,812 @@
+/*********************************************************************************/
+/*                                                                               */
+/*  Animation of heat equation in a planar domain                                */
+/*                                                                               */
+/*  N. Berglund, May 2021                                                        */
+/*                                                                               */
+/*  Feel free to reuse, but if doing so it would be nice to drop a               */
+/*  line to nils.berglund@univ-orleans.fr - Thanks!                              */
+/*                                                                               */
+/*  compile with                                                                 */
+/*  gcc -o heat heat.c                                                           */
+/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp        */
+/*                                                                               */
+/*  To make a video, set MOVIE to 1 and create subfolder tif_heat                */
+/*  It may be possible to increase parameter PAUSE                               */
+/*                                                                               */
+/*  create movie using                                                           */
+/*  ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4                             */
+/*                                                                               */
+/*********************************************************************************/
+
+/*********************************************************************************/
+/*                                                                               */
+/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
+/* One could make it much faster by using a GPU                                  */
+/*                                                                               */
+/*********************************************************************************/
+
+#include <math.h>
+#include <string.h>
+#include <GL/glut.h>
+#include <GL/glu.h>
+#include <unistd.h>
+#include <sys/types.h>
+#include <tiffio.h>     /* Sam Leffler's libtiff library. */
+#include <omp.h>
+
+#define MOVIE 0         /* set to 1 to generate movie */
+
+/* General geometrical parameters */
+
+#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 */
+
+/* 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 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 */
+
+/* Choice of the billiard table */
+
+#define B_DOMAIN 25      /* 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 D_MENGER 15     /* Menger-Sierpinski carpet */ 
+#define D_JULIA_INT 16  /* interior of Julia set */ 
+
+/* 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 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 */
+#define MDEPTH 2            /* 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 */
+#define FOCI 1              /* set to 1 to draw focal points of ellipse */
+
+/* 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.00001
+#define DT 0.000004
+// #define DT 0.000002
+// #define DT 0.00000002
+// #define DT 0.000000005
+#define VISCOSITY 10.0
+#define T_OUT 2.0       /* outside temperature */
+#define T_IN 0.0        /* inside temperature */
+// #define T_OUT 0.0       /* outside temperature */
+// #define T_IN 2.0        /* inside temperature */
+#define SPEED 0.0       /* speed of drift to the right */
+
+/* Boundary conditions */
+
+#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) */
+
+/* Parameters for length and speed of simulation */
+
+#define NSTEPS 4500      /* number of frames of movie */
+#define NVID 50          /* number of iterations between images displayed on screen */
+// #define NVID 100          /* number of iterations between images displayed on screen */
+#define NSEG 100         /* number of segments of boundary */
+
+#define PAUSE 100       /* number of frames after which to pause */
+#define PSLEEP 1         /* sleep time during pause */
+#define SLEEP1  2        /* initial sleeping time */
+#define SLEEP2  1        /* final sleeping time */
+
+/* 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 */
+
+/* Field representation */
+
+#define FIELD_REP 0
+
+#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 FIELD_LINE_FACTOR 100 /* factor controlling precision when computing origin of field lines */
+
+/* Color schemes */
+
+#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 */
+#define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
+
+#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 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 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 "sub_wave.c"
+
+double courant2;  /* Courant parameter squared */
+double dx2;       /* spatial step size squared */
+double intstep;   /* integration step */
+double intstep1;  /* integration step used in absorbing boundary conditions */
+
+
+
+void init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
+/* initialise field with gaussian at position (x,y) */
+    double x, y, mean, amplitude, scalex, *phi[NX]; 
+    short int * xy_in[NX];
+
+{
+    int i, j, in;
+    double xy[2], dist2, module, phase, scale2;    
+
+    scale2 = scalex*scalex;
+    printf("Initialising field\n");
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
+
+            in = xy_in[i][j];
+            if (in == 1)
+            {
+                dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
+                module = amplitude*exp(-dist2/scale2);
+                if (module < 1.0e-15) module = 1.0e-15;
+
+                phi[i][j] = mean + module/scalex;
+            }   /* boundary temperatures */
+            else if (in >= 2) phi[i][j] = T_IN*pow(0.75, (double)(in-2));
+//             else if (in >= 2) phi[i][j] = T_IN*pow(1.0 - 0.5*(double)(in-2), (double)(in-2));
+//             else if (in >= 2) phi[i][j] = T_IN*(1.0 - (double)(in-2)/((double)MDEPTH))*(1.0 - (double)(in-2)/((double)MDEPTH));
+            else phi[i][j] = T_OUT;
+        }
+}
+
+void init_julia_set(phi, xy_in)
+/* change Julia set boundary condition */
+    double *phi[NX]; 
+    short int * xy_in[NX];
+
+{
+    int i, j, in;
+    double xy[2], dist2, module, phase, scale2;    
+
+//     printf("Changing Julia set\n");
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
+
+            in = xy_in[i][j];
+            if (in >= 2) phi[i][j] = T_IN;
+        }
+}
+
+
+/*********************/
+/* animation part    */
+/*********************/
+
+
+void compute_gradient(phi, nablax, nablay)
+/* compute the gradient of the field */
+double *phi[NX], *nablax[NX], *nablay[NX];
+{
+    int i, j, iplus, iminus, jplus, jminus; 
+    double dx;
+    
+    dx = (XMAX-XMIN)/((double)NX);
+    
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            iplus = i+1;  if (iplus == NX) iplus = NX-1;
+            iminus = i-1; if (iminus == -1) iminus = 0;
+            jplus = j+1;  if (jplus == NX) jplus = NY-1;
+            jminus = j-1; if (jminus == -1) jminus = 0;
+            nablax[i][j] = (phi[iplus][j] - phi[iminus][j])/dx;
+            nablay[i][j] = (phi[i][jplus] - phi[i][jminus])/dx;
+        }
+}
+
+void draw_field_line(x, y, xy_in, nablax, nablay, delta, nsteps)
+// void draw_field_line(x, y, nablax, nablay, delta, nsteps)
+/* draw a field line of the gradient, starting in (x,y) */
+double x, y, *nablax[NX], *nablay[NX], delta;
+int nsteps;
+short int *xy_in[NX];
+{
+    double x1, y1, x2, y2, pos[2], nabx, naby, norm2, norm;
+    int i = 0, ij[2], cont = 1;
+    
+    glColor3f(1.0, 1.0, 1.0);
+    glLineWidth(FIELD_LINE_WIDTH);
+    x1 = x;
+    y1 = y;
+    
+//     printf("Drawing field line \n");
+
+    glEnable(GL_LINE_SMOOTH);
+    glBegin(GL_LINE_STRIP);
+    xy_to_pos(x1, y1, pos);
+    glVertex2d(pos[0], pos[1]);
+    
+    i = 0;
+    while ((cont)&&(i < nsteps))
+    {
+        xy_to_ij(x1, y1, ij);
+        
+        if (ij[0] < 0) ij[0] = 0;
+        if (ij[0] > NX-1) ij[0] = NX-1;
+        if (ij[1] < 0) ij[1] = 0;
+        if (ij[1] > NY-1) ij[1] = NY-1;
+        
+        nabx = nablax[ij[0]][ij[1]];
+        naby = nablay[ij[0]][ij[1]];
+        
+        norm2 = nabx*nabx + naby*naby;
+        
+        if (norm2 > 1.0e-14)
+        {
+            /* avoid too large step size */
+            if (norm2 < 1.0e-9) norm2 = 1.0e-9;
+            norm = sqrt(norm2);
+            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;            
+        }
+        
+        if (!xy_in[ij[0]][ij[1]]) cont = 0;
+        
+        /* stop if the boundary is hit */
+//         if (xy_in[ij[0]][ij[1]] != 1) cont = 0;
+        
+//         printf("x1 = %.3lg \t y1 = %.3lg \n", x1, y1);
+                
+        xy_to_pos(x1, y1, pos);
+        glVertex2d(pos[0], pos[1]);
+        
+        i++;
+    }
+    glEnd();
+}
+
+void draw_wave(phi, xy_in, scale, time)
+/* draw the field */
+double *phi[NX], scale;
+short int *xy_in[NX];
+int time;
+{
+    int i, j, iplus, iminus, jplus, jminus, ij[2], counter = 0;
+    static int first = 1;
+    double rgb[3], xy[2], x1, y1, x2, y2, dx, value, angle, dangle, intens, deltaintens, sum = 0.0;
+    double *nablax[NX], *nablay[NX];
+    static double linex[N_FIELD_LINES*FIELD_LINE_FACTOR], liney[N_FIELD_LINES*FIELD_LINE_FACTOR], distance[N_FIELD_LINES*FIELD_LINE_FACTOR], integral[N_FIELD_LINES*FIELD_LINE_FACTOR + 1];
+
+    for (i=0; i<NX; i++) 
+    {
+        nablax[i] = (double *)malloc(NY*sizeof(double));
+        nablay[i] = (double *)malloc(NY*sizeof(double));
+    }
+    
+    /* compute the gradient */
+//     if (FIELD_REP > 0) 
+    compute_gradient(phi, nablax, nablay);
+    
+    /* compute the position of origins of field lines */
+    if ((first)&&(DRAW_FIELD_LINES))
+    {
+        first = 0;
+        
+        printf("computing linex\n");
+        
+        x1 = sqrt(3.58);
+        y1 = 0.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);
+            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);
+    }
+
+    dx = (XMAX-XMIN)/((double)NX);
+    glBegin(GL_QUADS);
+
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            if (FIELD_REP == F_INTENSITY) value = phi[i][j];
+            else if (FIELD_REP == F_GRADIENT)
+            {
+                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);
+                glColor3f(rgb[0], rgb[1], rgb[2]);
+            }
+            else glColor3f(0.0, 0.0, 0.0);
+
+            glVertex2i(i, j);
+            glVertex2i(i+1, j);
+            glVertex2i(i+1, j+1);
+            glVertex2i(i, j+1);
+        }
+    glEnd ();
+        
+    /* draw a field line */
+    if (DRAW_FIELD_LINES)
+    {
+        /* compute gradient norm along boundary and its integral */
+        for (i = 0; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
+        {
+            xy_to_ij(linex[i], liney[i], ij);
+            intens = module2(nablax[ij[0]][ij[1]], nablay[ij[0]][ij[1]])*distance[i];
+            if (i > 0) integral[i] = integral[i-1] + intens;
+            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++;
+        }
+        printf("%i lines\n", counter);
+    }
+
+    
+    for (i=0; i<NX; i++)
+    {
+        free(nablax[i]);
+        free(nablay[i]);
+    }
+}
+
+
+
+void 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];;
+    
+    for (i=0; i<NX; i++) newphi[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++){
+            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[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
+
+                x = phi[i][j];
+
+                /* evolve phi */
+                if (B_COND != BC_ABSORBING)
+                {
+                    newphi[i][j] = x + intstep*(delta1 - SPEED*(phi[iplus][j] - phi[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))
+                    {
+                        newphi[i][j] = x - intstep*delta2;
+                    }
+                     /* right border */
+                    else if (i==NX-1) 
+                    {
+                        newphi[i][j] = x - intstep1*(x - phi[i-1][j]);
+                    }
+                    /* upper border */
+                    else if (j==NY-1) 
+                    {
+                        newphi[i][j] = x - intstep1*(x - phi[i][j-1]);
+                    }
+                    /* left border */
+                    else if (i==0) 
+                    {
+                        newphi[i][j] = x - intstep1*(x - phi[1][j]);
+                    }
+                   /* lower border */
+                    else if (j==0) 
+                    {
+                        newphi[i][j] = x - intstep1*(x - phi[i][1]);
+                    }
+                }
+
+
+                if (FLOOR)
+                {
+                    if (newphi[i][j] > VMAX) phi[i][j] = VMAX;
+                    if (newphi[i][j] < -VMAX) phi[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];
+        }
+    }
+    
+    for (i=0; i<NX; i++)
+    {
+        free(newphi[i]);
+    }
+
+//     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];
+{
+    int i, j, n = 0;
+    double variance = 0.0;
+
+    for (i=1; i<NX; i++)
+        for (j=1; j<NY; j++)
+        {
+            if (xy_in[i][j])
+            {
+                n++;
+                variance += phi[i][j]*phi[i][j];
+            }
+        }
+    if (n==0) n=1;
+    return(variance/(double)n);
+}
+
+void renormalise_field(phi, xy_in, variance)
+/* renormalise variance of field */
+double *phi[NX], variance; 
+short int * xy_in[NX];
+{
+    int i, j;
+    double stdv;
+    
+    stdv = sqrt(variance);
+
+    for (i=1; i<NX; i++)
+        for (j=1; j<NY; j++)
+        {
+            if (xy_in[i][j])
+            {
+                phi[i][j] = phi[i][j]/stdv;
+            }
+        }
+}
+
+void print_level(level)
+int level;
+{
+    double pos[2];
+    char message[50];
+    
+    glColor3f(1.0, 1.0, 1.0);
+    sprintf(message, "Level %i", level);
+    xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
+    write_text(pos[0], pos[1], message);
+}
+
+
+void print_Julia_parameters()
+{
+    double pos[2];
+    char message[50];
+    
+    glColor3f(1.0, 1.0, 1.0);
+    if (julia_y >= 0.0) sprintf(message, "c = %.5f + %.5f i", julia_x, julia_y);
+    else sprintf(message, "c = %.5f %.5f i", julia_x, julia_y);
+    xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
+    write_text(pos[0], pos[1], message);
+}
+
+void set_Julia_parameters(time, phi, xy_in)
+int time;
+double *phi[NX];
+short int *xy_in[NX];
+{
+    double jangle, cosj, sinj, radius = 0.15;
+
+    jangle = (double)time*DPI/(double)NSTEPS;
+//     jangle = (double)time*0.001;
+//     jangle = (double)time*0.0001;
+
+    cosj = cos(jangle);
+    sinj = sin(jangle);
+    julia_x = -0.9 + radius*cosj;
+    julia_y = radius*sinj;
+    init_julia_set(phi, xy_in);
+    
+    printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
+}
+
+void set_Julia_parameters_cardioid(time, phi, xy_in)
+int time;
+double *phi[NX];
+short int *xy_in[NX];
+{
+    double jangle, cosj, sinj, yshift;
+
+    jangle = pow(1.05 + (double)time*0.00003, 0.333);
+    yshift = 0.02*sin((double)time*PID*0.002);
+//     jangle = pow(1.0 + (double)time*0.00003, 0.333);
+//     jangle = pow(0.05 + (double)time*0.00003, 0.333);
+//     jangle = pow(0.1 + (double)time*0.00001, 0.333);
+//     yshift = 0.04*sin((double)time*PID*0.002);
+
+    cosj = cos(jangle);
+    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);
+    
+    printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
+}
+
+void animation()
+{
+    double time, scale, dx, var, jangle, cosj, sinj;
+    double *phi[NX];
+    short int *xy_in[NX];
+    int i, j, s;
+
+    /* Since NX and NY are big, it seemed wiser to use some memory allocation here */
+    for (i=0; i<NX; i++)
+    {
+        phi[i] = (double *)malloc(NY*sizeof(double));
+        xy_in[i] = (short int *)malloc(NY*sizeof(short int));
+    }
+
+    dx = (XMAX-XMIN)/((double)NX);
+    intstep = DT/(dx*dx*VISCOSITY);
+    intstep1 = DT/(dx*VISCOSITY);
+    
+//     julia_x = 0.1; 
+//     julia_y = 0.6; 
+    
+    set_Julia_parameters(0, phi, xy_in);
+    
+    printf("Integration step %.3lg\n", intstep);
+
+    /* initialize wave wave function */
+    init_gaussian(-1.0, 0.0, 0.1, 0.0, 0.01, phi, xy_in);
+//     init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
+    
+    if (SCALE)
+    {
+        var = compute_variance(phi, xy_in);
+        scale = sqrt(1.0 + var);
+        renormalise_field(phi, xy_in, var);
+    }
+
+    blank();
+    glColor3f(0.0, 0.0, 0.0);
+    
+
+    glutSwapBuffers();
+    
+    draw_wave(phi, xy_in, 1.0, 0);
+    draw_billiard();
+    print_Julia_parameters(i);
+    
+//     print_level(MDEPTH);
+
+    glutSwapBuffers();
+
+    sleep(SLEEP1);
+    if (MOVIE) for (i=0; i<SLEEP1*25; i++) save_frame();
+
+    for (i=0; i<=NSTEPS; i++)
+    {
+        /* compute the variance of the field to adjust color scheme */
+        /* the color depends on the field divided by sqrt(1 + variance) */
+        if (SCALE)
+        {
+            var = compute_variance(phi, xy_in);
+            scale = sqrt(1.0 + var);
+//             printf("Norm: %5lg\t Scaling factor: %5lg\n", var, scale);
+            renormalise_field(phi, xy_in, var);
+        }
+        else scale = 1.0;
+        
+        draw_wave(phi, xy_in, scale, i);
+        
+        for (j=0; j<NVID; j++) evolve_wave(phi, xy_in);
+
+        draw_billiard();
+        
+//         print_level(MDEPTH);
+        print_Julia_parameters(i);
+
+	glutSwapBuffers();
+        
+        /* modify Julia set */
+        set_Julia_parameters(i, phi, xy_in);
+
+	if (MOVIE)
+        {
+            save_frame();
+
+            /* 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   */
+            if (i % PAUSE == PAUSE - 1)
+            {
+                printf("Making a short pause\n");
+                sleep(PSLEEP);
+                s = system("mv wave*.tif tif_heat/");
+            }
+        }
+
+    }
+
+    if (MOVIE)
+    {
+        for (i=0; i<20; i++) save_frame();
+        s = system("mv wave*.tif tif_heat/");
+    }
+    for (i=0; i<NX; i++)
+    {
+        free(phi[i]);
+    }
+
+}
+
+
+void display(void)
+{
+    glPushMatrix();
+
+    blank();
+    glutSwapBuffers();
+    blank();
+    glutSwapBuffers();
+
+    animation();
+    sleep(SLEEP2);
+
+    glPopMatrix();
+
+    glutDestroyWindow(glutGetWindow());
+
+}
+
+
+int main(int argc, char** argv)
+{
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
+    glutInitWindowSize(WINWIDTH,WINHEIGHT);
+    glutCreateWindow("Heat equation in a planar domain");
+
+    init();
+
+    glutDisplayFunc(display);
+
+    glutMainLoop();
+
+    return 0;
+}
+
diff --git a/particle_billiard.c b/particle_billiard.c
index 1806733..c884e99 100644
--- a/particle_billiard.c
+++ b/particle_billiard.c
@@ -39,6 +39,9 @@
 #define XMAX 2.0	/* x interval */
 #define YMIN -1.125
 #define YMAX 1.125	/* 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 */
+
 // #define XMIN -1.8
 // #define XMAX 1.8	/* x interval */
 // #define YMIN -0.91
@@ -57,31 +60,39 @@
 #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 LAMBDA 5.0	/* parameter controlling shape of billiard */
-#define LAMBDA -1.949855824	/* 7-sided Reuleaux triangle */
+#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 */
+// #define LAMBDA -1.949855824	/* 7-sided Reuleaux triangle */
 // #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 FOCI 1          /* set to 1 to draw focal points of ellipse */
-#define NPOLY 7             /* number of sides of polygon */
-#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define NPOLY 6             /* number of sides of polygon */
+#define APOLY 2.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 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 10000     /* number of particles */
+#define NPART 20000     /* 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 4000     /* number of frames of movie */
-#define TIME 125       /* time between movie frames, for fluidity of real-time simulation */ 
-#define DPHI 0.00001    /* integration step */
+#define NSTEPS 6000     /* 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 NVID 150         /* number of iterations between images displayed on screen */
 
 /* Decreasing TIME accelerates the animation and the movie                               */
@@ -92,18 +103,19 @@
 
 /* Colors and other graphical parameters */
 
-#define NCOLORS -7      /* number of colors */
-#define COLORSHIFT 220    /* hue of initial color */ 
-#define NSEG 100        /* number of segments of boundary */
-#define LENGTH 0.01     /* length of velocity vectors */
-#define BILLIARD_WIDTH 4    /* width of billiard */
-#define PARTICLE_WIDTH 2       /* width of particles */
+#define NCOLORS -10      /* number of colors */
+#define COLORSHIFT 200     /* hue of initial color */ 
+#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 PARTICLE_WIDTH 2    /* width of particles */
 #define FRONT_WIDTH 3       /* width of wave front */
 
 #define BLACK 1             /* set to 1 for black background */
 #define COLOR_OUTSIDE 0     /* set to 1 for colored outside */ 
 #define OUTER_COLOR 270.0   /* color outside billiard */
-#define PAINT_INT 1         /* set to 1 to paint interior in other color (for polygon) */
+#define PAINT_INT 0         /* set to 1 to paint interior in other color (for polygon/Reuleaux) */
 
 
 #define PAUSE 1000       /* number of frames after which to pause */
@@ -230,6 +242,7 @@ double *configs[NPARTMAX];
 
     glutSwapBuffers(); 
     blank();
+    if (PAINT_INT) paint_billiard_interior();
       
     glLineWidth(PARTICLE_WIDTH);
     
@@ -262,40 +275,40 @@ double *configs[NPARTMAX];
             glColor3f(rgb[0], rgb[1], rgb[2]);
         
             glBegin(GL_LINE_STRIP);
-            glVertex2d(x1, y1);
-            glVertex2d(x2, y2);
+            glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*y1);
+            glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*y2);
             glEnd ();
         
             /* taking care of boundary conditions - only needed for periodic boundary conditions */
-            if (x2 > XMAX)
+            if (SCALING_FACTOR*x2 > XMAX)
             {
                 glBegin(GL_LINE_STRIP);
-                glVertex2d(x1+XMIN-XMAX, y1);
-                glVertex2d(x2+XMIN-XMAX, y2);
+                glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
+                glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
                 glEnd ();
             }
         
-            if (x2 < XMIN)
+            if (SCALING_FACTOR*x2 < XMIN)
             {
                 glBegin(GL_LINE_STRIP);
-                glVertex2d(x1-XMIN+XMAX, y1);
-                glVertex2d(x2-XMIN+XMAX, y2);
+                glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
+                glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
                 glEnd ();
             }
 
-            if (y2 > YMAX)
+            if (SCALING_FACTOR*y2 > YMAX)
             {
                 glBegin(GL_LINE_STRIP);
-                glVertex2d(x1, y1+YMIN-YMAX);
-                glVertex2d(x2, y2+YMIN-YMAX);
+                glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
+                glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
                 glEnd ();
             }
 
-            if (y2 < YMIN)
+            if (SCALING_FACTOR*y2 < YMIN)
             {
                 glBegin(GL_LINE_STRIP);
-                glVertex2d(x1, y1+YMAX-YMIN);
-                glVertex2d(x2, y2+YMAX-YMIN);
+                glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
+                glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
                 glEnd ();
             }
         }
@@ -305,15 +318,15 @@ double *configs[NPARTMAX];
         {
             glLineWidth(1.0);
             glBegin(GL_LINES);
-            glVertex2d(configs[i][4], configs[i][5]);
-            glVertex2d(configs[i][6], configs[i][7]);
+            glVertex2d(SCALING_FACTOR*configs[i][4], SCALING_FACTOR*configs[i][5]);
+            glVertex2d(SCALING_FACTOR*configs[i][6], SCALING_FACTOR*configs[i][7]);
             glEnd ();
             glLineWidth(3.0);
         }
     
         if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
     }
-    draw_billiard();    
+    if (DRAW_BILLIARD) draw_billiard();    
 }
 
 
@@ -355,9 +368,9 @@ double *configs[NPARTMAX];
 
 void animation()
 {
-    double time, dt, alpha;
+    double time, dt, alpha, r;
     double *configs[NPARTMAX];
-    int i, j, resamp = 1, s;
+    int i, j, resamp = 1, s, i1, i2;
     int *color, *newcolor, *active;
     
     /* Since NPARTMAX can be big, it seemed wiser to use some memory allocation here */
@@ -368,19 +381,21 @@ void animation()
         configs[i] = (double *)malloc(8*sizeof(double));
       
     /* initialize system by putting particles in a given point with a range of velocities */
-    alpha = 3.0*PI/7.0;
-    init_sym_drop_config(-0.99, 0.0, -alpha, alpha, configs);
+    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.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);
-//     init_drop_config(0.0, 0.5, 0.0, DPI, configs);
 
 //  other possible initial conditions :
 //     init_line_config(-0.6, 0.2, -0.6, 0.7, 0.0, configs);
 //     init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs);
-//     init_line_config(0.0, 0.1, 0.0, 0.7, 0.0, configs);
+    init_line_config(0.0, -0.3, 0.0, 0.3, 0.0, configs);
   
     blank();  
     glColor3f(0.0, 0.0, 0.0);
-    draw_billiard();
+    if (DRAW_BILLIARD) draw_billiard();
   
     glutSwapBuffers();   
   
@@ -391,6 +406,24 @@ void animation()
         newcolor[i] = 0;
         active[i] = 1;
     }
+    
+    if (FLOWER_COLOR)   /* adapt color scheme to flower configuration (beta implementation) */
+    {
+//         i1 = (int)((double)NPART*0.2538);     /* the 0.27 is just a trial-and-error guess, to be improved */
+//         i1 = (int)((double)NPART*0.1971);     /* the 0.27 is just a trial-and-error guess, to be improved */
+        i1 = (int)((double)NPART*0.3015);     /* the 0.27 is just a trial-and-error guess, to be improved */
+        i2 = NPART-i1;
+        for (i=i1; i<i2; i++) 
+        {
+            color[i] += NCOLORS/3;
+            newcolor[i] = NCOLORS/3;
+        }
+        for (i=i2; i<NPART; i++) 
+        {
+            color[i] += 2*NCOLORS/3;
+            newcolor[i] = 2*NCOLORS/3;
+        }
+    }
   
     sleep(SLEEP1);
   
@@ -399,7 +432,7 @@ void animation()
         graph_movie(TIME, newcolor, configs, active);
         
         draw_config(newcolor, configs, active);
-        draw_billiard();
+        if (DRAW_BILLIARD) draw_billiard();
         for (j=0; j<NPARTMAX; j++) color[j] = newcolor[j];
 
         
diff --git a/schrodinger.c b/schrodinger.c
index a5109bf..edea4c0 100644
--- a/schrodinger.c
+++ b/schrodinger.c
@@ -53,6 +53,8 @@
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
+#define JULIA_SCALE 1.0 /* scaling for Julia sets */
+
 /* Choice of the billiard table */
 
 #define B_DOMAIN 3      /* choice of domain shape */
@@ -70,10 +72,26 @@
 #define D_GRATING 10    /* diffraction grating */
 #define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
 
+#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
+#define D_JULIA_INT 16  /* interior of Julia set */ 
+
+/* 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 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 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 */
 
 /* You can add more billiard tables by adapting the functions */
@@ -81,7 +99,8 @@
 
 /* Physical patameters of wave equation */
 
-#define DT 0.00000002
+#define DT 0.00000005
+// #define DT 0.00000002
 // #define DT 0.000000005
 #define HBAR 1.0
 
@@ -96,7 +115,7 @@
 /* Parameters for length and speed of simulation */
 
 #define NSTEPS 4500      /* number of frames of movie */
-#define NVID 750         /* number of iterations between images displayed on screen */
+#define NVID 250         /* 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 */
@@ -145,6 +164,8 @@
 #define DPI 	6.283185307
 #define PID 	1.570796327
 
+double julia_x = 0.0, julia_y = 0.0;    /* parameters for Julia sets */
+
 #include "sub_wave.c"
 
 double courant2;  /* Courant parameter squared */
@@ -154,7 +175,6 @@ double intstep1;  /* integration step used in absorbing boundary conditions */
 
 
 
-
 void init_coherent_state(x, y, px, py, scalex, phi, psi, xy_in)
 /* initialise field with coherent state of position (x,y) and momentum (px, py) */
 /* phi is real part, psi is imaginary part */
@@ -252,10 +272,18 @@ int time;
 void evolve_wave(phi, psi, 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[NX], *psi[NX]; 
+    short int *xy_in[NX];
 {
     int i, j, iplus, iminus, jplus, jminus;
-    double delta1, delta2, x, y;
+    double delta1, delta2, x, y, *newphi[NX], *newpsi[NX];
+    
+    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++){
@@ -288,54 +316,69 @@ void evolve_wave(phi, psi, xy_in)
                 /* evolve phi and psi */
                 if (B_COND != BC_ABSORBING)
                 {
-                    phi[i][j] = x - intstep*delta2;
-                    psi[i][j] = y + intstep*delta1;
+                    newphi[i][j] = x - intstep*delta2;
+                    newpsi[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))
                     {
-                        phi[i][j] = x - intstep*delta2;
-                        psi[i][j] = y + intstep*delta1;
+                        newphi[i][j] = x - intstep*delta2;
+                        newpsi[i][j] = y + intstep*delta1;
                     }
                      /* right border */
                     else if (i==NX-1) 
                     {
-                        phi[i][j] = x - intstep1*(y - psi[i-1][j]);
-                        psi[i][j] = y + intstep1*(x - phi[i-1][j]);
+                        newphi[i][j] = x - intstep1*(y - psi[i-1][j]);
+                        newpsi[i][j] = y + intstep1*(x - phi[i-1][j]);
                     }
                     /* upper border */
                     else if (j==NY-1) 
                     {
-                        phi[i][j] = x - intstep1*(y - psi[i][j-1]);
-                        psi[i][j] = y + intstep1*(x - phi[i][j-1]);
+                        newphi[i][j] = x - intstep1*(y - psi[i][j-1]);
+                        newpsi[i][j] = y + intstep1*(x - phi[i][j-1]);
                     }
                     /* left border */
                     else if (i==0) 
                     {
-                        phi[i][j] = x - intstep1*(y - psi[1][j]);
-                        psi[i][j] = y + intstep1*(x - phi[1][j]);
+                        newphi[i][j] = x - intstep1*(y - psi[1][j]);
+                        newpsi[i][j] = y + intstep1*(x - phi[1][j]);
                     }
                    /* lower border */
                     else if (j==0) 
                     {
-                        phi[i][j] = x - intstep1*(y - psi[i][1]);
-                        psi[i][j] = y + intstep1*(x - phi[i][1]);
+                        newphi[i][j] = x - intstep1*(y - psi[i][1]);
+                        newpsi[i][j] = y + intstep1*(x - phi[i][1]);
                     }
                 }
 
 
                 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;
+                    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;
                 }
             }
         }
     }
+    
+    
+    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]);
 }
 
@@ -409,12 +452,15 @@ void animation()
 //     init_coherent_state(-0.5, 0.0, 1.0, 1.0, 0.05, phi, psi, xy_in);
     
     
+    
     if (SCALE)
     {
         var = compute_variance(phi,psi, xy_in);
         scale = sqrt(1.0 + var);
         renormalise_field(phi, psi, xy_in, var);
     }
+    
+    
 
     blank();
     glColor3f(0.0, 0.0, 0.0);
@@ -439,12 +485,14 @@ void animation()
         else scale = 1.0;
 
         draw_wave(phi, psi, xy_in, scale, i);
-
+        
+//         printf("Wave drawn\n");
+        
         for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
-
+        
         draw_billiard();
-
-	glutSwapBuffers();
+        
+        glutSwapBuffers();
 
 	if (MOVIE)
         {
diff --git a/sub_part_billiard.c b/sub_part_billiard.c
index 851d113..c96e175 100644
--- a/sub_part_billiard.c
+++ b/sub_part_billiard.c
@@ -239,10 +239,45 @@ void write_text( double x, double y, char *st)
 } 
 
 
-void paint_billiard_interior()      /* points billiard interior, for use before draw_conf */
+ void compute_flower_parameters(omega, co, so, axis1, axis2, phimax)
+ /* compute parameters needed for the flower billiard in terms of LAMBDA and NPOLY */
+ double *omega, *co, *so, *axis1, *axis2, *phimax;
+ {
+    double omega2, co2, so2, r, a, gamma, axissquare1;
+    
+    /* various angles */
+    *omega = DPI/((double)NPOLY);
+    omega2 = PI/((double)NPOLY);
+    co2 = cos(omega2);
+    so2 = sin(omega2);
+    *co = cos(*omega);
+    *so = sin(*omega);
+//     *co = co2*co2 - so2*so2;
+//     *so = 2.0*co2*so2;
+    
+    /* distance of edge of ellipse to the origin */
+    r = LAMBDA*co2/(*co);
+    
+    a = (r*co2 - *co)*(r*co2 - *co);
+    gamma = 0.5*r*r - r*co2*(*co) + 0.5*cos(2.0*(*omega));
+    axissquare1 = gamma + sqrt(gamma*gamma + a*(*so)*(*so));
+    
+    /* semi-minor axis */
+    *axis1 = sqrt(axissquare1);
+    
+    /* semi-major axis */
+    *axis2 = sqrt(axissquare1 + (*so)*(*so));
+    
+    /* max angle in ellipse parametrization */
+    *phimax = asin(r*so2/(*axis2));
+ }
+
+
+ 
+void paint_billiard_interior()      /* paints billiard interior, for use before draw_conf */
 {
-    double x0, x, y, phi, r = 0.01, alpha, dphi, omega;
-    int i, k, c;
+    double x0, x, y, phi, r = 0.01, alpha, dphi, omega, beta2, x2, s, x1, y1, angle, co, so, axis1, axis2, phimax;
+    int i, j, k, c;
     
     glLineWidth(4);
     
@@ -273,6 +308,64 @@ void paint_billiard_interior()      /* points billiard interior, for use before
             }
             break;
         }
+        case (D_REULEAUX):
+        {
+            omega = DPI/((double)NPOLY);
+            beta2 = asin(sin(omega*0.5)/LAMBDA);
+            if (LAMBDA > 0.0) x2 = cos(omega*0.5) + sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
+            else x2 = cos(omega*0.5) - sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
+            
+            if (PAINT_INT)
+            {
+                if (BLACK) glColor3f(1.0, 1.0, 1.0);
+                else glColor3f(0.0, 0.0, 0.0);
+                
+                glBegin(GL_TRIANGLE_FAN);
+                glVertex2d(0.0, 0.0);
+                for (i=0; i<=NPOLY; i++)
+                {
+                    for (j=0; j<NSEG; j++)
+                {
+                    s = 2.0*(((double)j/(double)NSEG)-0.5)*beta2;
+                    x1 = x2 - LAMBDA*cos(s);
+                    y1 = LAMBDA*sin(s);
+                    angle = i*omega + APOLY*PID;
+                    x = cos(angle)*x1 - sin(angle)*y1;
+                    y = sin(angle)*x1 + cos(angle)*y1;
+                    glVertex2d(x, y);
+                }
+               }
+                glEnd();
+            }
+            break;
+        }
+        case (D_FLOWER): 
+        {
+            compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
+            if (PAINT_INT)
+            {
+                if (BLACK) glColor3f(1.0, 1.0, 1.0);
+                else glColor3f(0.0, 0.0, 0.0);
+                
+                glBegin(GL_TRIANGLE_FAN);
+                glVertex2d(0.0, 0.0);
+                for (i=0; i<=NPOLY; i++)
+                {
+                    for (j=0; j<NSEG; j++)
+                    {
+                        s = 2.0*(((double)j/(double)NSEG)-0.5)*phimax;
+                        x1 = co + axis1*cos(s);
+                        y1 = axis2*sin(s);
+                        angle = i*omega + APOLY*PID;
+                        x = cos(angle)*x1 - sin(angle)*y1;
+                        y = sin(angle)*x1 + cos(angle)*y1;
+                        glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
+                    }
+               }
+                glEnd();
+            }
+            break;
+        }
         default: 
         {
             
@@ -282,7 +375,8 @@ void paint_billiard_interior()      /* points billiard interior, for use before
 
 void draw_billiard()      /* draws the billiard boundary */
 {
-    double x0, x, y, phi, r = 0.01, alpha, dphi, omega, x1, y1, x2, beta2, angle, s;
+    double x0, x, y, phi, r = 0.01, alpha, dphi, omega, x1, y1, x2, beta2, angle, s, x2plus, x2minus;
+    double omega2, co, so, axis1, axis2, phimax;
     int i, j, k, c;
     
     if (PAINT_INT) glColor3f(0.5, 0.5, 0.5);
@@ -366,6 +460,20 @@ void draw_billiard()      /* draws the billiard boundary */
                 glVertex2d(x, y);
             }
             glEnd();
+            
+            if (DRAW_CONSTRUCTION_LINES)
+            {
+                glColor3f(0.5, 0.5, 0.5);
+                glBegin(GL_LINE_STRIP);
+                glVertex2d(-LAMBDA, -1.0);
+                glVertex2d(-LAMBDA, 1.0);
+                glEnd();
+                glBegin(GL_LINE_STRIP);
+                glVertex2d(LAMBDA, -1.0);
+                glVertex2d(LAMBDA, 1.0);
+                glEnd();
+            }
+            
             break; 
         }
         case D_SINAI:
@@ -580,6 +688,86 @@ void draw_billiard()      /* draws the billiard boundary */
             glEnd ();
             break;
         }
+        case (D_FLOWER):
+        {
+            compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
+
+            /* draw inner polygon and radial lines */
+            if (DRAW_CONSTRUCTION_LINES)
+            {
+                glColor3f(0.5, 0.5, 0.5);
+                glBegin(GL_LINE_LOOP);
+                for (i=0; i<=NPOLY; i++)
+                {
+                    x = cos(i*omega + APOLY*PID);
+                    y = sin(i*omega + APOLY*PID);
+                    glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
+                }
+                glEnd ();
+
+                r = LAMBDA*cos(0.5*omega)/co;
+                for (i=0; i<=NPOLY; i++)
+                {
+                    glBegin(GL_LINE_STRIP);
+                    glVertex2d(0.0, 0.0);
+                    x = r*cos(((double)i + 0.5)*omega + APOLY*PID);
+                    y = r*sin(((double)i + 0.5)*omega + APOLY*PID);
+                    glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
+                    glEnd ();
+                }
+            }
+            
+            /* draw billiard boundary */
+            if (!PAINT_INT) 
+            {
+                if (BLACK) glColor3f(1.0, 1.0, 1.0);
+                else glColor3f(0.0, 0.0, 0.0);
+            }
+            glBegin(GL_LINE_STRIP);
+            for (i=0; i<=NPOLY; i++)
+//             for (i=0; i<=1; i++)
+            {
+                for (j=0; j<NSEG; j++)
+                {
+//                     s = 2.0*(((double)j/(double)NSEG)-0.5)*PI;
+                    s = 2.0*(((double)j/(double)NSEG)-0.5)*phimax;
+                    x1 = co + axis1*cos(s);
+                    y1 = axis2*sin(s);
+                    angle = i*omega + APOLY*PID;
+                    x = cos(angle)*x1 - sin(angle)*y1;
+                    y = sin(angle)*x1 + cos(angle)*y1;
+                    glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
+                }
+            }
+            glEnd ();
+            
+            
+            break;
+        }
+        case (D_ALT_REU):
+        {
+            omega = DPI/((double)NPOLY);
+            beta2 = asin(sin(omega*0.5)/LAMBDA);
+            x2plus = cos(omega*0.5) + sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
+            x2minus = cos(omega*0.5) - sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
+            glBegin(GL_LINE_STRIP);
+            for (i=0; i<=NPOLY; i++)
+            {
+                for (j=0; j<NSEG; j++)
+                {
+                    s = 2.0*(((double)j/(double)NSEG)-0.5)*beta2;
+                    if (i%2==0) x1 = x2plus - LAMBDA*cos(s);
+                    else x1 = x2minus + LAMBDA*cos(s);
+                    y1 = LAMBDA*sin(s);
+                    angle = i*omega + APOLY*PID;
+                    x = cos(angle)*x1 - sin(angle)*y1;
+                    y = sin(angle)*x1 + cos(angle)*y1;
+                    glVertex2d(x, y);
+                }
+            }
+            glEnd ();
+            break;
+        }
         default: 
         {
             printf("Function draw_billiard not defined for this billiard \n");
@@ -1688,57 +1876,38 @@ int color[NPARTMAX];
             rangle = 2.0*(double)k*omega2 + APOLY*PID;
             theta = alpha - rangle;
             
-//             if ((intb))     /* check if condition is ok */
-//             if ((cos(theta) > 0.0)&&(intb))     /* check if condition is ok */
+            ca = cos(rangle);
+            sa = sin(rangle);
+                                
+            x = pos[0]*ca + pos[1]*sa;
+            y = -pos[0]*sa + pos[1]*ca;
+                
+            a = (x-x2)*cos(theta) + y*sin(theta);
+            b = (x-x2)*(x-x2) + y*y - LAMBDA*LAMBDA;
+                
+            if (a*a - b > margin)
             {
-                ca = cos(rangle);
-                sa = sin(rangle);
-                
-//                 printf("theta = %.5lg\n", theta);
-//                 printf("rangle = %.5lg x0 = %.5lg y0 = %.5lg \n", rangle, pos[0], pos[1]);
-                
-                x = pos[0]*ca + pos[1]*sa;
-                y = -pos[0]*sa + pos[1]*ca;
-                
-//                 printf("x = %.5lg\t y = %.5lg\n", x, y);
-            
-                a = (x-x2)*cos(theta) + y*sin(theta);
-                b = (x-x2)*(x-x2) + y*y - LAMBDA*LAMBDA;
-                
-//                 printf("a = %.5lg\t b = %.5lg\n", a, b);
-                
-                if (a*a - b > margin)
+                if (LAMBDA > 0.0) t = -a - sqrt(a*a - b);
+                else t = -a + sqrt(a*a - b);
+                    
+                xi = x + t*cos(theta);
+                yi = y + t*sin(theta);
+                    
+                if ((t > margin)&&(vabs(yi) <= sin(omega2))) 
                 {
-//                     t = vabs(a) - sqrt(a*a - b);
-                    if (LAMBDA > 0.0) t = -a - sqrt(a*a - b);
-                    else t = -a + sqrt(a*a - b);
-                    
-                    xi = x + t*cos(theta);
-                    yi = y + t*sin(theta);
-                    
-//                     printf("t = %.5lg\t xi = %.5lg\t yi = %.5lg\n", t, xi, yi);
-                
-                    if ((t > margin)&&(vabs(yi) <= sin(omega2))) 
-                    {
-                        cval[nt] = k;
-                        tval[nt] = t;
+                    cval[nt] = k;
+                    tval[nt] = t;
                         
-                        /* rotate back */
-                        x1[nt] = xi*ca - yi*sa;
-                        y1[nt] = xi*sa + yi*ca;
+                    /* rotate back */
+                    x1[nt] = xi*ca - yi*sa;
+                    y1[nt] = xi*sa + yi*ca;
                     
-//                         intb = 0;
-//                         c = k;
-                        tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/LAMBDA);
-                        tempconf[nt][1] = PID - asin(yi/LAMBDA) - theta;      
-//                         tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/vabs(LAMBDA));
-//                         tempconf[nt][1] = PID - asin(yi/vabs(LAMBDA)) - theta;      
-                        nt++;
-                    }
+                    tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/LAMBDA);
+                    tempconf[nt][1] = PID - asin(yi/LAMBDA) - theta;      
+                    nt++;
                 }
             }
         }
-//         printf("nt = %i\n", nt);
         
         /* find earliest intersection */
         tmin = tval[0];
@@ -1783,6 +1952,390 @@ int color[NPARTMAX];
 	
 	return(c);
  }
+
+/****************************************************************************************/
+/* Bunimovich flower billiard */
+/****************************************************************************************/
+
+ int pos_flower(conf, pos, alpha)
+ /* determine position on boundary of polygon */
+ /* conf[0] is arclength on boundary, it belongs to [0,2*NPOLY*phimax) */
+ double conf[2], pos[2], *alpha;
+
+ {
+	double s, theta, omega, co, so, axis1, axis2, phimax, s1, x, y, angle;
+        int c;
+
+	s = conf[0]; 
+        theta = conf[1];
+        
+        compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
+        
+        c = (int)(s/(2.0*phimax));         /* side of shape */
+        
+        s1 = s - (((double)c)*2.0 + 1.0)*phimax;
+        
+        x = co + axis1*cos(s1);
+        y = axis2*sin(s1);
+        
+        angle = ((double)c)*omega + PID*APOLY;
+//         angle = 2.0*((double)c)*omega + PID*APOLY;
+        
+        pos[0] = x*cos(angle) - y*sin(angle);
+        pos[1] = x*sin(angle) + y*cos(angle);
+        
+        *alpha = argument(-axis1*sin(s1), axis2*cos(s1)) + theta + angle;
+        
+//         printf("alpha = %.5lg\t", *alpha);
+        
+        return(c);
+ }
+ 
+ int vflower_xy(config, alpha, pos)
+ /* determine initial configuration for start at point pos = (x,y) */
+ double config[8], alpha, pos[2];
+
+ {
+	double s, theta, omega, omega2, s1, rangle, x, y, x1, y1, xi, yi, t;
+        double ca, sa, aa, bb, cc, margin = 1.0e-14, tmin;
+        double co, so, co2, so2, ct, st, phimax, phi, axis1, axis2;
+	int k, c, intb=1, intc, i, nt = 0, ntmin, sign;
+
+        compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
+        
+        for (k=0; k<NPOLY; k++) if (intb)
+        {
+            /* rotate position so that kth side is vertical */
+//             rangle = (double)(2*k)*omega + APOLY*PID;
+            rangle = (double)k*omega + APOLY*PID;
+            theta = alpha - rangle;
+            
+            ca = cos(rangle);
+            sa = sin(rangle);
+            
+            ct = cos(theta);
+            st = sin(theta);
+                                
+            x = pos[0]*ca + pos[1]*sa;
+            y = -pos[0]*sa + pos[1]*ca;
+                
+            /* find intersection with elliptical arc */
+            aa = ct*ct/(axis1*axis1) + st*st/(axis2*axis2);
+            bb = (x-co)*ct/(axis1*axis1) + y*st/(axis2*axis2);
+            cc = (x-co)*(x-co)/(axis1*axis1) + y*y/(axis2*axis2) - 1.0;
+                
+//             if (bb*bb - aa*cc > margin) 
+            if (bb*bb - aa*cc >= 0.0) 
+            {
+                t = (-bb + sqrt(bb*bb - aa*cc))/aa;
+                    
+                xi = x + t*cos(theta);
+                yi = y + t*sin(theta);
+                
+                if (yi >= 0.0) phi = argument((xi - co)/axis1, yi/axis2);
+                else phi = -argument((xi - co)/axis1, -yi/axis2);
+                
+                if ((t > margin)&&((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax))) 
+                {
+                    intb = 0;
+                    c = k;
+                        
+                    /* rotate back */
+                    x1 = xi*ca - yi*sa;
+                    y1 = xi*sa + yi*ca;
+                    
+                    config[0] = (double)(2*k + 1)*phimax + phi;
+                    config[1] = argument(-axis1*sin(phi), axis2*cos(phi)) - theta;      
+                }
+            }
+        }
+        
+//        if (nt == 0) printf("nt = %i\t ntmin = %i \tcmin = %i\n", nt, ntmin, c);
+
+        if (config[1] < 0.0) config[1] += DPI;
+
+        config[2] = 0.0;	/* running time */ 
+	config[3] = module2(x1-pos[0], y1-pos[1]);     /* distance to collision */
+	config[4] = pos[0];    /* start position */
+	config[5] = pos[1];
+	config[6] = x1;        /* position of collision */
+	config[7] = y1;
+	
+	return(c);
+ }
+
+  int old_vflower_xy(config, alpha, pos)
+ /* determine initial configuration for start at point pos = (x,y) */
+ double config[8], alpha, pos[2];
+
+ {
+	double s, theta, omega, omega2, s1, rangle, x, y, x1[2*NPOLY], y1[2*NPOLY], xi, yi, t;
+        double ca, sa, aa, bb, cc, margin = 1.0e-14, tmin, tval[2*NPOLY], tempconf[2*NPOLY][2];
+        double co, so, co2, so2, ct, st, phimax, phi, axis1, axis2;
+	int k, c, intb=1, intc, i, nt = 0, cval[2*NPOLY], ntmin, sign;
+
+        compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
+        
+        for (k=0; k<NPOLY; k++)
+        {
+            /* rotate position so that kth side is vertical */
+//             rangle = (double)(2*k)*omega + APOLY*PID;
+            rangle = (double)k*omega + APOLY*PID;
+            theta = alpha - rangle;
+            
+            ca = cos(rangle);
+            sa = sin(rangle);
+            
+            ct = cos(theta);
+            st = sin(theta);
+                                
+            x = pos[0]*ca + pos[1]*sa;
+            y = -pos[0]*sa + pos[1]*ca;
+                
+            /* find intersection with elliptical arc */
+            aa = ct*ct/(axis1*axis1) + st*st/(axis2*axis2);
+            bb = (x-co)*ct/(axis1*axis1) + y*st/(axis2*axis2);
+            cc = (x-co)*(x-co)/(axis1*axis1) + y*y/(axis2*axis2) - 1.0;
+                
+//             if (bb*bb - aa*cc > margin) 
+            if (bb*bb - aa*cc >= 0.0) 
+            {
+                t = (-bb + sqrt(bb*bb - aa*cc))/aa;
+                    
+                xi = x + t*cos(theta);
+                yi = y + t*sin(theta);
+                
+                if (yi >= 0.0) phi = argument((xi - co)/axis1, yi/axis2);
+                else phi = -argument((xi - co)/axis1, -yi/axis2);
+                
+//                 phi = argument((xi - co)/axis1, yi/axis2);
+//                 if (phi > PI) phi += -DPI;
+                    
+                if ((t > margin)&&((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax))) 
+//                 if (((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax))) 
+//                 if ((t > margin)) 
+                {
+                    cval[nt] = k;
+//                     cval[nt] = 2*k;
+                    tval[nt] = t;
+                        
+                    /* rotate back */
+                    x1[nt] = xi*ca - yi*sa;
+                    y1[nt] = xi*sa + yi*ca;
+                    
+                    tempconf[nt][0] = (double)(2*k + 1)*phimax + phi;
+                    tempconf[nt][1] = argument(-axis1*sin(phi), axis2*cos(phi)) - theta;      
+                    nt++;
+                }
+            }
+        }
+        
+        /* find earliest intersection */
+        tmin = tval[0];
+        ntmin = 0;
+        for (i=1; i<nt; i++) 
+            if (tval[i] < tmin) 
+            {
+                tmin = tval[i];
+                ntmin = i;
+            }
+            
+        config[0] = tempconf[ntmin][0];
+        config[1] = tempconf[ntmin][1];
+        c = cval[ntmin];
+ 
+       if (nt == 0) printf("nt = %i\t ntmin = %i \tcmin = %i\n", nt, ntmin, c);
+        
+
+        if (config[1] < 0.0) config[1] += DPI;
+
+        config[2] = 0.0;	/* running time */ 
+	config[3] = module2(x1[ntmin]-pos[0], y1[ntmin]-pos[1]);     /* distance to collision */
+	config[4] = pos[0];    /* start position */
+	config[5] = pos[1];
+	config[6] = x1[ntmin];        /* position of collision */
+	config[7] = y1[ntmin];
+	
+	return(c);
+ }
+ 
+ int vflower(config)
+ /* determine initial configuration when starting from boundary */
+ double config[8];
+
+ {
+	double pos[2], alpha;
+	int c;
+
+        c = pos_flower(config, pos, &alpha);
+        
+        vflower_xy(config, alpha, pos);
+	
+	return(c);
+ } 
+ 
+ /****************************************************************************************/
+/* Alternating between Reuleaux-type and star-shaped billiard */
+/****************************************************************************************/
+
+ int pos_alt_reuleaux(conf, pos, alpha)
+ /* determine position on boundary of polygon */
+ /* conf[0] is arclength on boundary */
+ double conf[2], pos[2], *alpha;
+
+ {
+	double s, theta, omega2, beta2, beta, s1, angle, x2plus, x2minus, x, y;
+        int c;
+
+	s = conf[0]; 
+        theta = conf[1];
+        
+        omega2 = PI/((double)NPOLY);
+        beta2 = asin(sin(omega2)/vabs(LAMBDA));
+        beta = beta2*2.0;
+
+        c = (int)(s/beta);         /* side of shape */
+        
+        s1 = s - ((double)c)*beta;
+        
+        x2plus = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
+        x2minus = cos(omega2) - sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
+        
+        if (c%2 == 0) x = x2plus - LAMBDA*cos(s1 - beta2);
+        else x = x2minus + LAMBDA*cos(s1 - beta2);
+        
+        if (c%2 == 0) y = LAMBDA*sin(s1 - beta2);
+        else y = -LAMBDA*sin(s1 - beta2);
+        
+        /* test, to be removed */
+//         x = x2plus - LAMBDA*cos(s1 - beta2);
+//         y = LAMBDA*sin(s1 - beta2);
+        
+        angle = 2.0*((double)c)*omega2 + PID*APOLY;
+        
+        pos[0] = x*cos(angle) - y*sin(angle);
+        pos[1] = x*sin(angle) + y*cos(angle);
+        
+        *alpha = PID - s1 + beta2 + theta + 2.0*(double)c*omega2 + APOLY*PID;
+        
+//         printf("alpha = %.5lg\t", *alpha);
+        
+        return(c);
+ }
+ 
+ int valt_reuleaux_xy(config, alpha, pos)
+ /* determine initial configuration for start at point pos = (x,y) */
+ double config[8], alpha, pos[2];
+
+ {
+	double s, theta, omega2, beta, s1, rangle, x, y, x1[NPOLY], y1[NPOLY], xi, yi, t, x2plus, x2minus, arcsine;
+        double ca, sa, a, b, margin = 1.0e-14, tmin, tval[NPOLY], tempconf[NPOLY][2];
+	int k, c, intb=1, intc, i, nt = 0, cval[NPOLY], ntmin;
+
+        /* dimensions/angles of polygon */
+        omega2 = PI/((double)NPOLY);
+        beta = 2.0*asin(sin(omega2)/vabs(LAMBDA));
+//         printf("beta = %.5lg\n", beta);
+        
+        x2plus = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
+        x2minus = cos(omega2) - sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
+//         printf("x2 = %.5lg\n", x2);
+        
+        for (k=0; k<NPOLY; k++)
+        {
+            /* rotate position so that kth side is vertical */
+            rangle = 2.0*(double)k*omega2 + APOLY*PID;
+            theta = alpha - rangle;
+            
+            ca = cos(rangle);
+            sa = sin(rangle);
+                                
+            x = pos[0]*ca + pos[1]*sa;
+            y = -pos[0]*sa + pos[1]*ca;
+                
+            if (k%2==0)
+            {
+                a = (x-x2plus)*cos(theta) + y*sin(theta);
+                b = (x-x2plus)*(x-x2plus) + y*y - LAMBDA*LAMBDA;
+            }
+            else 
+            {
+                a = (x-x2minus)*cos(theta) + y*sin(theta);
+                b = (x-x2minus)*(x-x2minus) + y*y - LAMBDA*LAMBDA;
+            }
+    
+            if (a*a - b > margin)
+            {
+                if (k%2==0) t = -a - sqrt(a*a - b);
+                else t = -a + sqrt(a*a - b);
+                    
+                xi = x + t*cos(theta);
+                yi = y + t*sin(theta);
+                    
+                if ((t > margin)&&(vabs(yi) <= sin(omega2))) 
+                {
+                    cval[nt] = k;
+                    tval[nt] = t;
+                        
+                    /* rotate back */
+                    x1[nt] = xi*ca - yi*sa;
+                    y1[nt] = xi*sa + yi*ca;
+                    
+                    if (k%2==0) arcsine = asin(yi/LAMBDA);
+                    else arcsine = -asin(yi/LAMBDA);
+                    
+                    tempconf[nt][0] = ((double)k + 0.5)*beta + arcsine;
+                    tempconf[nt][1] = PID - arcsine - theta;      
+                    nt++;
+                }
+            }
+        }
+        
+        /* find earliest intersection */
+        tmin = tval[0];
+        ntmin = 0;
+        for (i=1; i<nt; i++) 
+            if (tval[i] < tmin) 
+            {
+                tmin = tval[i];
+                ntmin = i;
+            }
+            
+        config[0] = tempconf[ntmin][0];
+        config[1] = tempconf[ntmin][1];
+        c = cval[ntmin];
+ 
+//         printf("nt = %i\t ntmin = %i \tcmin = %i\n", nt, ntmin, c);
+        
+
+        if (config[1] < 0.0) config[1] += DPI;
+
+        config[2] = 0.0;	/* running time */ 
+	config[3] = module2(x1[ntmin]-pos[0], y1[ntmin]-pos[1]);     /* distance to collision */
+	config[4] = pos[0];    /* start position */
+	config[5] = pos[1];
+	config[6] = x1[ntmin];        /* position of collision */
+	config[7] = y1[ntmin];
+	
+	return(c);
+ }
+
+ int valt_reuleaux(config)
+ /* determine initial configuration when starting from boundary */
+ double config[8];
+
+ {
+	double pos[2], alpha;
+	int c;
+
+        c = pos_alt_reuleaux(config, pos, &alpha);
+        
+        valt_reuleaux_xy(config, alpha, pos);
+	
+	return(c);
+ }
+
+
  
 /****************************************************************************************/
 /* general billiard */
@@ -1833,6 +2386,16 @@ int color[NPARTMAX];
             return(pos_reuleaux(conf, pos, &alpha));
             break;
         }
+        case (D_FLOWER):
+        {
+            return(pos_flower(conf, pos, &alpha));
+            break;
+        }
+        case (D_ALT_REU):
+        {
+            return(pos_alt_reuleaux(conf, pos, &alpha));
+            break;
+        }
         default: 
         {
             printf("Function pos_billiard not defined for this billiard \n");
@@ -1887,6 +2450,16 @@ int color[NPARTMAX];
             return(vreuleaux_xy(config, alpha, pos));
             break;
         }
+        case (D_FLOWER):
+        {
+            return(vflower_xy(config, alpha, pos));
+            break;
+        }
+        case (D_ALT_REU):
+        {
+            return(valt_reuleaux_xy(config, alpha, pos));
+            break;
+        }
         default: 
         {
             printf("Function vbilliard_xy not defined for this billiard \n");
@@ -1960,6 +2533,20 @@ int color[NPARTMAX];
             return(vreuleaux(config, alpha, pos));
             break;
         }
+        case (D_FLOWER):
+        {
+            c = pos_flower(config, pos, &alpha);
+        
+            return(vflower(config, alpha, pos));
+            break;
+        }
+        case (D_ALT_REU):
+        {
+            c = pos_alt_reuleaux(config, pos, &alpha);
+        
+            return(valt_reuleaux(config, alpha, pos));
+            break;
+        }
         default: 
         {
             printf("Function vbilliard not defined for this billiard \n");
@@ -1971,7 +2558,7 @@ int color[NPARTMAX];
  /* returns 1 if (x,y) represents a point in the billiard */
  double x, y;
  {
-    double l2, r2, omega, omega2, c, angle, x1, y1, x2, co, so;
+    double l2, r1, r2, omega, omega2, c, angle, x1, y1, x2, co, so, x2plus, x2minus;
     int condition, k;
  
     switch (B_DOMAIN) {
@@ -2021,8 +2608,9 @@ int color[NPARTMAX];
         case D_ANNULUS:
         {
             l2 = LAMBDA*LAMBDA;
-            r2 = x*x + y*y;
-            if ((r2 > l2)&&(r2 < 1.0)) return(1);
+            r1 = x*x + y*y;
+            r2 = (x-MU)*(x-MU) + y*y;
+            if ((r2 > l2)&&(r1 < 1.0)) return(1);
             else return(0);
             break;
         }
@@ -2056,22 +2644,42 @@ int color[NPARTMAX];
                 y1 = -x*sin(angle) + y*cos(angle);
                 if (LAMBDA > 0.0) condition = condition*((x1-x2)*(x1-x2) + y1*y1 > LAMBDA*LAMBDA);
                 else condition = condition*((x1-x2)*(x1-x2) + y1*y1 < LAMBDA*LAMBDA);
-                
-//                 if (!condition)
-//                 {
-//                     printf("x = %.5lg \t y = %.5lg \t x1 = %.5lg \t y1 = %.5lg \t angle = %.5lg \n", x, y, x1, y1, angle);
-//                     printf("k = %i \t condition = %i\n", k, condition);
-//                     sleep(1);
-//                 }
             }
             return(condition);
             break;            
         }
         /* D_REULEAUX : distance to all centers of arcs should be larger than LAMBDA */
+        case D_FLOWER:
+        {
+            /* TO DO */
+            return(1);
+            break;
+        }
+        case D_ALT_REU:
+        {
+            condition = 1;
+            omega2 = PI/((double)NPOLY);
+            co = cos(omega2);
+            so = sin(omega2);
+            x2plus = co + sqrt(LAMBDA*LAMBDA - so*so);
+            x2minus = co - sqrt(LAMBDA*LAMBDA - so*so);
+                        
+            for (k=0; k<NPOLY; k++)  
+            {
+                angle = 2.0*(double)k*omega2 + APOLY*PID;
+                
+                x1 = x*cos(angle) + y*sin(angle);
+                y1 = -x*sin(angle) + y*cos(angle);
+                if (k%2==0) condition = condition*((x1-x2plus)*(x1-x2plus) + y1*y1 > LAMBDA*LAMBDA);
+                else condition = condition*((x1-x2minus)*(x1-x2minus) + y1*y1 < LAMBDA*LAMBDA);
+            }
+            return(condition);
+            break;            
+        }
         default: 
         {
             printf("Function ij_in_billiard not defined for this billiard \n");
-            return(0);
+            return(1);
         }
     }
  }
diff --git a/sub_wave.c b/sub_wave.c
index 58c921d..187b31b 100644
--- a/sub_wave.c
+++ b/sub_wave.c
@@ -198,6 +198,7 @@ void write_text( double x, double y, char *st)
 
 
 
+
 /*********************/
 /* some basic math   */
 /*********************/
@@ -296,6 +297,22 @@ double xy[2];
     xy[1] = YMIN + ((double)j)*(YMAX-YMIN)/((double)NY);
 }
 
+void draw_rectangle(x1, y1, x2, y2)
+double x1, y1, x2, y2;
+{
+    double pos[2];
+    
+    glBegin(GL_LINE_LOOP);
+    xy_to_pos(x1, y1, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x2, y1, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x2, y2, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x1, y2, pos);
+    glVertex2d(pos[0], pos[1]);
+    glEnd();    
+}
 
 
 
@@ -303,7 +320,7 @@ int xy_in_billiard(x, y)
 /* returns 1 if (x,y) represents a point in the billiard */
 double x, y;
 {
-    double l2, r2, omega, c, angle, z;
+    double l2, r2, r2mu, omega, c, angle, z, x1, y1, u, v, u1, v1;
     int i, k, k1, k2, condition;
 
     switch (B_DOMAIN) {
@@ -399,6 +416,124 @@ 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_MENGER:       
+        {
+            x1 = 0.5*(x+1.0);
+            y1 = 0.5*(y+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_JULIA_INT:  
+        {
+            u = x/JULIA_SCALE;
+            v = y/JULIA_SCALE;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + julia_x;
+                v = 2.0*u*v + julia_y;
+                u = u1;
+                i++;
+            }
+            if (u*u + v*v < MANDELLIMIT) return(1);
+            else return(0);
+        }
+        case D_ANNULUS_HEATED:      /* returns 2 if in inner circle */ 
+        {
+            l2 = LAMBDA*LAMBDA;
+            r2 = x*x + y*y;
+            r2mu = (x-MU)*(x-MU) + y*y;
+            if ((r2mu > l2)&&(r2 < 1.0)) return(1);
+            else if (r2mu <= l2) return(2);
+            else return (0);
+        }
+        case D_MENGER_HEATED: 
+        {
+            if ((vabs(x) >= 1.0)||(vabs(y) >= 1.0)) return(0);
+            else
+            {
+                x1 = 0.5*(x+1.0);
+                y1 = 0.5*(y+1.0);
+                for (k=0; k<MDEPTH; k++)
+                {
+                    x1 = x1*(double)MRATIO;
+                    y1 = y1*(double)MRATIO;
+                    if ((((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(k+2);
+                }
+                return(1);
+            }
+        }
+        case D_MENGER_H_OPEN:       /* returns 2 if in inner circle */ 
+        {
+            x1 = 0.5*(x+1.0);
+            y1 = 0.5*(y+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(k+2);
+            }
+            return(1);
+        }
+        case D_MANDELBROT:  
+        {
+            u = 0.0;
+            v = 0.0;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + x;
+                v = 2.0*u*v + y;
+                u = u1;
+                i++;
+                /* old version used */
+                /* u1 = u*u - v*v - x; */
+                /* v = 2.0*u*v - y; */
+            }
+            if (u*u + v*v < MANDELLIMIT) return(0);
+            else if ((x-0.5)*(x-0.5)/3.0 + y*y/1.0 > 1.2) return(2);
+            else return(1);
+        }
+        case D_MANDELBROT_CIRCLE:  
+        {
+            u = 0.0;
+            v = 0.0;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + x;
+                v = 2.0*u*v + y;
+                u = u1;
+                i++;
+            }
+            if (u*u + v*v < MANDELLIMIT) return(0);
+            else if ((x-LAMBDA)*(x-LAMBDA) + (y-0.5)*(y-0.5) < MU*MU) return(2);
+            else return(1);
+        }
+        case D_JULIA:  
+        {
+            u = x/JULIA_SCALE;
+            v = y/JULIA_SCALE;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + julia_x;
+                v = 2.0*u*v + julia_y;
+                u = u1;
+                i++;
+//                 printf("x = %.5lg y = %.5lg i = %i r2 = %.5lg\n", x, y, i, u*u+v*v);
+            }
+//             printf("i = %i x = %.5lg y = %.5lg r2 = %.5lg\n", i, x, y, u*u+v*v);
+            if (u*u + v*v < MANDELLIMIT) return(0);
+            else if (x*x/3.0 + y*y/1.0 > 1.2) return(2);
+//             else if ((vabs(x) > XMAX - 0.01)||(vabs(y) > YMAX - 0.01)) return(2);
+            else return(1);
+        }
         default:
         {
             printf("Function ij_in_billiard not defined for this billiard \n");
@@ -423,8 +558,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;
-    int i, j, k1, k2;
+    double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
+    int i, j, k1, k2, mr2;
 
     if (BLACK) glColor3f(1.0, 1.0, 1.0);
     else glColor3f(0.0, 0.0, 0.0);
@@ -710,7 +845,195 @@ void draw_billiard()      /* draws the billiard boundary */
             glEnd ();
             break;
         }
-        default:
+        case (D_MENGER):
+        {
+            glLineWidth(3);
+//             draw_rectangle(XMIN, -1.0, XMAX, 1.0);
+            
+            /* level 1 */
+            if (MDEPTH > 0)
+            {
+                glLineWidth(2);
+                x = 1.0/((double)MRATIO);
+                draw_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_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_rectangle(x, y, x+l, y+l);
+                        }
+            }
+            
+            break;
+        }        
+         case (D_JULIA_INT):
+        {
+            /* Do nothing */
+            break;
+        }
+        case (D_ANNULUS_HEATED):
+        {
+            glBegin(GL_LINE_LOOP);
+            for (i=0; i<=NSEG; i++)
+            {
+                phi = (double)i*DPI/(double)NSEG;
+                x = MU + LAMBDA*cos(phi);
+                y = LAMBDA*sin(phi);
+                xy_to_pos(x, y, pos);
+                glVertex2d(pos[0], pos[1]);
+            }
+            glEnd ();
+            glBegin(GL_LINE_LOOP);
+            for (i=0; i<=NSEG; i++)
+            {
+                phi = (double)i*DPI/(double)NSEG;
+                x = cos(phi);
+                y = sin(phi);
+                xy_to_pos(x, y, pos);
+                glVertex2d(pos[0], pos[1]);
+            }
+            glEnd ();      
+            break;
+        }
+        case (D_MENGER_HEATED):
+        {
+            glLineWidth(3);
+            draw_rectangle(-1.0, -1.0, 1.0, 1.0);
+            
+            /* level 1 */
+            if (MDEPTH > 0)
+            {
+                glLineWidth(2);
+                x = 1.0/((double)MRATIO);
+                draw_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_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_rectangle(x, y, x+l, y+l);
+                        }
+            }
+            
+            break;
+        }
+        case (D_MENGER_H_OPEN):
+        {
+            glLineWidth(3);
+//             draw_rectangle(XMIN, -1.0, XMAX, 1.0);
+            
+            /* level 1 */
+            if (MDEPTH > 0)
+            {
+                glLineWidth(2);
+                x = 1.0/((double)MRATIO);
+                draw_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_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_rectangle(x, y, x+l, y+l);
+                        }
+            }
+            
+            break;
+        }        
+        case (D_MANDELBROT):
+        {
+            /* Do nothing */
+            break;
+        }
+        case (D_MANDELBROT_CIRCLE):
+        {
+            /* Do nothing */
+            break;
+        }
+        case (D_JULIA):
+        {
+            /* Do nothing */
+            break;
+        }
+       default:
         {
             printf("Function draw_billiard not defined for this billiard \n");
         }
diff --git a/wave_billiard.c b/wave_billiard.c
index b262cee..6ddfea5 100644
--- a/wave_billiard.c
+++ b/wave_billiard.c
@@ -46,8 +46,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                                 */
@@ -57,9 +59,11 @@
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
+#define JULIA_SCALE 1.1 /* scaling for Julia sets */
+
 /* Choice of the billiard table */
 
-#define B_DOMAIN 8      /* choice of domain shape */
+#define B_DOMAIN 16      /* choice of domain shape */
 
 #define D_RECTANGLE 0   /* rectangular domain */
 #define D_ELLIPSE 1     /* elliptical domain */
@@ -74,21 +78,39 @@
 #define D_GRATING 10    /* diffraction grating */
 #define D_EHRENFEST 11  /* Ehrenfest urn type geometry */
 
+#define D_MENGER 15     /* Menger-Sierpinski carpet */ 
+#define D_JULIA_INT 16  /* interior of Julia set */ 
+
+/* 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 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 */
 
 /* You can add more billiard tables by adapting the functions */
 /* xy_in_billiard and draw_billiard below */
 
-/* Physical patameters of wave equation */
+/* Physical parameters of wave equation */
 
+#define OMEGA 0.9           /* frequency of periodic excitation */
 #define COURANT 0.01       /* Courant number */
 #define GAMMA 0.0      /* damping factor in wave equation */
 // #define GAMMA 5.0e-10      /* damping factor in wave equation */
-#define KAPPA 5.0e-6       /* "elasticity" term enforcing oscillations */
+#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 */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
@@ -96,13 +118,22 @@
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
+/* Boundary conditions */
+
+#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 */
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 5000      /* number of frames of movie */
+#define NSTEPS 4000      /* number of frames of movie */
 #define NVID 25          /* number of iterations between images displayed on screen */
 #define NSEG 100         /* number of segments of boundary */
 
@@ -120,7 +151,7 @@
 #define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
 #define C_HUE 1          /* color scheme modifies hue */
 
-#define SCALE 1          /* set to 1 to adjust color scheme to variance of field */
+#define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
 
@@ -128,8 +159,8 @@
 #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 100.0    /* mean value of hue for color scheme C_HUE */
-#define HUEAMP 80.0      /* amplitude of variation of hue for color scheme C_HUE */
+#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 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 */
 
@@ -139,6 +170,10 @@
 #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 "sub_wave.c"
 
 double courant2;  /* Courant parameter squared */
@@ -162,6 +197,24 @@ void init_wave(x, y, phi, psi, xy_in)
         }
 }
 
+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];
+
+{
+    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);
+	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
+	    phi[i][j] = 0.0;
+            psi[i][j] = 0.0;
+        }
+}
+
 void add_drop_to_wave(factor, x, y, phi, psi)
 /* add drop at (x,y) to the field with given prefactor */
 double factor, x, y, *phi[NX], *psi[NX];
@@ -178,6 +231,33 @@ 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) */
+double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
+{
+    int i, j, ij1[2], ij2[2], imin, imax, jmin, jmax, d = 5;
+    double xy[2], dist2;
+    
+    xy_to_ij(x1, y1, ij1);
+    xy_to_ij(x2, y2, ij2);
+    imin = ij1[0] - d;     if (imin < 0) imin = 0;
+    imax = ij2[0] + d;     if (imax >= NX) imax = NX-1;
+    jmin = ij1[1] - d;     if (jmin < 0) jmin = 0;
+    jmax = ij2[1] + d;     if (jmax >= NY) jmax = NY-1;
+
+    for (i = imin; i < imax; i++)
+        for (j = jmin; j < jmax; j++)
+        {
+            ij_to_xy(i, j, xy);
+            dist2 = (xy[0]-x1)*(xy[0]-x1);      /* to be improved */
+//             dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
+//             if (dist2 < 0.01)
+            if (dist2 < 0.001)
+                phi[i][j] = amplitude*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01)*cos(t*OMEGA);
+//                 phi[i][j] += 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01)*cos(t*OMEGA);
+        }
+}
+
 
 
 
@@ -215,7 +295,92 @@ 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 value of field at time t, psi at time t-1 */
+    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 delta, x, y, c, cc;
+    
+    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 ((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;
+                }
+                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 ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
+                        phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(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);
+                    
+                    /* lower border */
+                    else if (j==0) 
+                        phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA*x - GAMMA*(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);
+                    
+                    /* left border */
+                    else if (i==0) 
+                        phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA*x - GAMMA*(x-y);
+                }
+                psi_out[i][j] = x;
+
+                if (FLOOR)
+                {
+                    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;
+                }
+            }
+        }
+    }
+//     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);
+}
+
+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];
@@ -286,7 +451,7 @@ double compute_variance(phi, psi, xy_in)
 void animation()
 {
     double time, scale;
-    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;
 
@@ -295,12 +460,15 @@ 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));
     }
 
     courant2 = COURANT*COURANT;
 
     /* 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);
 
     /* add a drop at another point */
@@ -329,12 +497,19 @@ void animation()
             scale = sqrt(1.0 + compute_variance(phi,psi, xy_in));
 //             printf("Scaling factor: %5lg\n", scale);
         }
-
         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, xy_in);
+        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();
 
 
@@ -365,6 +540,9 @@ void animation()
     {
         free(phi[i]);
         free(psi[i]);
+        free(phi_tmp[i]);
+        free(psi_tmp[i]);
+        free(xy_in[i]);
     }
 
 }