diff --git a/global_perc.c b/global_perc.c
new file mode 100644
index 0000000..fc968ca
--- /dev/null
+++ b/global_perc.c
@@ -0,0 +1,67 @@
+/* Global variables and parameters for wave_billiard, heat and schrodinger */
+
+/* Basic math */
+
+#define PI 	3.141592654
+#define DPI 	6.283185307
+#define PID 	1.570796327
+
+/* Lattice types */
+
+#define BC_SQUARE_DIRICHLET 0   /* square lattice, Dirichlet boundary conditions */
+#define BC_SQUARE_PERIODIC 1    /* square lattice, periodic boundary conditions */
+#define BC_SQUARE_BOND_DIRICHLET 2  /* square lattice, bond percolation, Dirichlet b.c. */
+#define BC_HEX_SITE_DIRICHLET 10    /* hexagonal lattice, site percolation, Dirichlet b.c. */
+#define BC_HEX_BOND_DIRICHLET 12    /* hexagonal lattice, bond percolation, Dirichlet b.c. */
+#define BC_TRIANGLE_SITE_DIRICHLET 20   /* triangular lattice, site percolation, Dirichlet b.c. */
+
+/* Plot types */
+
+#define PLOT_SQUARES 0      /* plot squares */
+#define PLOT_SQUARE_BONDS 1 /* plot edges of square lattice */
+#define PLOT_HEX 2          /* plot hexagons */
+#define PLOT_TRIANGLE 3     /* plot triangles */
+#define PLOT_HEX_BONDS 4    /* plot edges of hexagonal lattice */
+
+/* Color schemes */
+
+#define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
+#define C_HUE 1          /* color scheme modifies hue */
+#define C_PHASE 2        /* color scheme shows phase */
+#define C_ONEDIM 3       /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */
+#define C_ONEDIM_LINEAR 4   /* use preset 1d color scheme with linear scale */
+
+/* Color palettes */
+
+#define COL_JET 0       /* JET color palette */
+#define COL_HSLUV 1     /* HSLUV color palette (perceptually uniform) */
+#define COL_GRAY 2      /* grayscale */
+
+#define COL_TURBO 10     /* TURBO color palette (by Anton Mikhailov) */
+#define COL_VIRIDIS 11   /* Viridis color palette */
+#define COL_MAGMA 12     /* Magma color palette */
+#define COL_INFERNO 13   /* Inferno color palette */
+#define COL_PLASMA 14    /* Plasma color palette */
+#define COL_CIVIDIS 15   /* Cividis color palette */
+#define COL_PARULA 16    /* Parula color palette */
+#define COL_TWILIGHT 17  /* Twilight color palette */
+#define COL_TWILIGHT_SHIFTED 18  /* Shifted twilight color palette */
+
+#define COL_TURBO_CYCLIC 101    /* TURBO color palette (by Anton Mikhailov) corrected to be cyclic, beta */
+
+
+typedef struct
+{
+    double proba;           /* probability of being open */
+    short int nneighb;      /* number of neighbours */
+    int nghb[6];            /* indices of neighbours */
+    short int open;         /* vertex is open, flooded */
+    short int flooded;      /* vertex is open, flooded */
+    short int active;       /* used in cluster algorithm */
+    int cluster;            /* number of cluster, for search of all clusters */
+    int previous_cluster;   /* number of cluster of previous p, to limit number of color changes */
+    int tested;             /* 1 if the site has been tested for belonging to a cluster */
+} t_perco;
+
+
+
diff --git a/percolation.c b/percolation.c
new file mode 100644
index 0000000..aaa473c
--- /dev/null
+++ b/percolation.c
@@ -0,0 +1,288 @@
+/*********************************************************************************/
+/*                                                                               */
+/*  Simulation of percolation in 2D                                              */
+/*                                                                               */
+/*  N. Berglund, July 2022                                                       */
+/*                                                                               */
+/*  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 percolation percolation.c                                             */
+/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp        */
+/*                                                                               */
+/*  OMP acceleration may be more effective after executing                       */
+/*  export OMP_NUM_THREADS=2 in the shell before running the program             */
+/*                                                                               */
+/*  To make a video, set MOVIE to 1 and create subfolder tif_perc                */
+/*  It may be possible to increase parameter PAUSE                               */
+/*                                                                               */
+/*  create movie using                                                           */
+/*  ffmpeg -i perc.%05d.tif -vcodec libx264 perc.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>
+#include <time.h>
+
+#define MOVIE 0         /* set to 1 to generate movie */
+
+/* General geometrical parameters */
+
+// #define WINWIDTH 	1920  /* window width */
+// #define WINHEIGHT 	1000  /* window height */
+// #define NX 1920          /* number of grid points on x axis */
+// #define NY 992           /* number of grid points on y axis */
+// 
+// #define XMIN -2.0
+// #define XMAX 2.0	/* x interval  */
+// #define YMIN -1.041666667
+// #define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
+
+#define HIGHRES 0        /* set to 1 if resolution of grid is double that of displayed image */
+
+#define WINWIDTH 	1280  /* window width */
+#define WINHEIGHT 	720   /* window height */
+
+#define NX 1280          /* number of grid points on x axis */
+#define NY 720          /* number of grid points on y axis */
+
+#define XMIN -2.0
+#define XMAX 2.0	/* x interval  */
+#define YMIN -1.125
+#define YMAX 1.125	/* y interval for 9/16 aspect ratio */
+
+/* Boundary conditions, see list in global_pdes.c  */
+
+#define LATTICE 10
+
+#define FLOOD_LEFT_BOUNDARY 1   /* set to 1 to flood cells on left boundary */
+#define FIND_ALL_CLUSTERS 0     /* set to 1 to find all open clusters */
+
+#define PLOT_CLUSTER_SIZE 1     /* set to 1 to add a plot for the size of the percolation cluster */
+#define PLOT_CLUSTER_NUMBER 0   /* set to 1 to add a graph of the number of clusters */
+#define PLOT_CLUSTER_HISTOGRAM 0    /* set to 1 to add a histogram of the number of clusters */
+#define PRINT_LARGEST_CLUSTER_SIZE 0    /* set to 1 to print size of largest cluster */
+
+#define MAX_CLUSTER_NUMBER 6    /* vertical scale of the cluster number plot */
+#define HISTO_BINS 30           /* number of bins in histrogram */
+
+#define NSTEPS 1270        /* number of frames of movie */
+
+#define DEBUG 0         /* set to 1 for some debugging features */
+#define DEBUG_SLEEP_TIME 1  /* sleep time between frames when debugging */
+#define TEST_GRAPH 0   /* set to 1 to test graph connectivity matrix */
+#define TEST_START 2210    /* start position of connectivity test */
+
+#define PAUSE 200       /* number of frames after which to pause */
+#define PSLEEP 2         /* sleep time during pause */
+#define SLEEP1  1        /* initial sleeping time */
+#define SLEEP2  1        /* final sleeping time */
+#define MID_FRAMES 20    /* number of still frames between parts of two-part movie */
+#define END_FRAMES 100    /* number of still frames at end of movie */
+#define FADE 1           /* set to 1 to fade at end of movie */
+
+/* Color schemes */
+
+#define COLOR_PALETTE 15     /* Color palette, see list in global_pdes.c  */
+
+#define BLACK 1          /* background */
+
+#define COLOR_CLUSTERS_BY_SIZE 0    /* set to 1 to link cluster color to their size */
+
+#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 */
+
+#define HUE_CLOSED 330.0   /* color hue of closed cells */
+#define HUE_OPEN 200.0     /* color hue of open (dry) cells */
+#define HUE_FLOODED 45.0   /* color hue of open flooded cells */
+#define HUE_GRAPH 150.0    /* color hue in graph of cluster size */
+
+#define CLUSTER_HUEMIN 0.0     /* minimal color hue of clusters */
+#define CLUSTER_HUEMAX 250.0    /* maximal color hue of clusters */
+#define N_CLUSTER_COLORS 20     /* number of different colors of clusters */
+
+#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 180.0    /* mean value of hue for color scheme C_HUE */
+#define HUEAMP -180.0      /* amplitude of variation of hue for color scheme C_HUE */
+
+#define ADD_PLOT ((PLOT_CLUSTER_SIZE)||(PLOT_CLUSTER_NUMBER)||(PLOT_CLUSTER_HISTOGRAM))
+#define FIND_CLUSTER_SIZES ((COLOR_CLUSTERS_BY_SIZE)||(PLOT_CLUSTER_HISTOGRAM))
+
+#include "global_perc.c"        /* constants and global variables */
+#include "sub_perco.c"          
+
+
+
+void animation(int size)
+{
+    int i, j, k, s, nx, ny, ncells, nmaxcells, maxsize, nopen, nflooded, nstack, nclusters, maxclustersize = 0, maxclusterlabel;
+    int *plot_cluster_number, *cluster_sizes;
+    double p, *plot_cluster_size;
+    t_perco *cell;
+    t_perco **pstack;
+    
+    compute_nxny(size, &nx, &ny);
+    
+    nmaxcells = cell_number(NX, NY);
+    
+    cell = (t_perco *)malloc(nmaxcells*sizeof(t_perco));
+    if (PLOT_CLUSTER_SIZE) plot_cluster_size = (double *)malloc(NSTEPS*sizeof(double));
+    if (PLOT_CLUSTER_NUMBER) plot_cluster_number = (int *)malloc(NSTEPS*sizeof(double));
+    if (FIND_CLUSTER_SIZES) cluster_sizes = (int *)malloc(2*nmaxcells*sizeof(int));
+    
+    ncells = init_cell_lattice(cell, nx, ny);
+    printf("nx = %i, ny = %i, ncells = %i\n", nx, ny, ncells);
+    pstack = (t_perco* *)malloc(ncells*sizeof(struct t_perco *));
+    
+    if (TEST_GRAPH) test_neighbours(TEST_START, cell, nx, ny, size, ncells);
+    
+    init_cell_probabilities(cell, ncells);
+    
+    for (i=0; i<NSTEPS; i++)
+    {
+        p = p_schedule(i);
+        printf("\ni = %i, p = %.4lg\n", i, p);
+        
+        init_cell_state(cell, p, ncells, (i == 0));
+        
+        if (FLOOD_LEFT_BOUNDARY) nstack = init_flooded_cells(cell, nx, ny, pstack);
+        nopen = count_open_cells(cell, ncells);
+        
+        if (FLOOD_LEFT_BOUNDARY) nflooded = find_percolation_cluster(cell, ncells, pstack, nstack);
+    
+        if (FIND_ALL_CLUSTERS) 
+        {
+            nclusters = find_all_clusters(cell, ncells, (i == 0), &maxclusterlabel);
+            printf("Found %i clusters\n", nclusters);
+        }
+        
+        if (FIND_CLUSTER_SIZES) 
+        {
+            maxclustersize = find_cluster_sizes(cell, ncells, cluster_sizes, &maxclusterlabel); 
+            printf("Max cluster size %i, max cluster label %i, ncells %i\n", maxclustersize, maxclusterlabel, ncells);
+        }
+        
+//         print_cluster_sizes(cell, ncells, cluster_sizes);
+        
+        draw_configuration(cell, cluster_sizes, nx, ny, size, ncells);
+        print_p(p);
+        if (PRINT_LARGEST_CLUSTER_SIZE) print_largest_cluster_size(maxclustersize);
+        
+        printf("%i open cells, %i flooded cells out of %i cells\n", nopen, nflooded, ncells);
+        
+        if (PLOT_CLUSTER_SIZE) 
+        {
+            plot_cluster_size[i] = (double)(nflooded)/(double)(nopen);
+            draw_size_plot(plot_cluster_size, i,  pcritical(LATTICE));
+        }
+        
+        if (PLOT_CLUSTER_NUMBER) 
+        {
+            plot_cluster_number[i] = nclusters;
+            draw_cluster_number_plot(plot_cluster_number, ncells/MAX_CLUSTER_NUMBER, i,  pcritical(LATTICE));
+            print_nclusters(nclusters);
+        }
+        
+        if (PLOT_CLUSTER_HISTOGRAM) draw_cluster_histogram(ncells, cluster_sizes, maxclustersize, maxclusterlabel);
+                
+        glutSwapBuffers();
+        
+        if (DEBUG) 
+        {
+            printf("\n\n");
+            sleep(DEBUG_SLEEP_TIME);
+        }
+        
+        if (MOVIE) save_frame_perc();
+    }
+    
+    if (MOVIE) 
+    {
+        for (i=0; i<MID_FRAMES; i++) save_frame_perc();
+        s = system("mv perc*.tif tif_perc/");
+    }
+    
+    free(cell);
+    free(pstack);
+    
+    if (PLOT_CLUSTER_SIZE) free(plot_cluster_size);
+    if (PLOT_CLUSTER_NUMBER) free(plot_cluster_number);
+    if (FIND_CLUSTER_SIZES) free(cluster_sizes);
+}
+
+
+
+
+void display(void)
+{
+    time_t rawtime;
+    struct tm * timeinfo;
+
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    
+    glPushMatrix();
+
+    blank();
+    glutSwapBuffers();
+    blank();
+    glutSwapBuffers();
+
+//     animation(128);
+//     animation(64);
+    animation(32);
+    animation(16);
+    animation(8);
+    animation(4);
+    animation(2);
+    animation(1);
+    
+    
+    sleep(SLEEP2);
+
+    glPopMatrix();
+
+    glutDestroyWindow(glutGetWindow());
+    
+    printf("Start local time and date: %s", asctime(timeinfo));
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    printf("Current local time and date: %s", asctime(timeinfo));
+}
+
+
+int main(int argc, char** argv)
+{
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
+    glutInitWindowSize(WINWIDTH,WINHEIGHT);
+    glutCreateWindow("Percolation in a planar domain");
+
+    init();
+
+    glutDisplayFunc(display);
+
+    glutMainLoop();
+
+    return 0;
+}
+
diff --git a/sub_perco.c b/sub_perco.c
new file mode 100644
index 0000000..f1f5484
--- /dev/null
+++ b/sub_perco.c
@@ -0,0 +1,2388 @@
+/*********************/
+/* Graphics routines */
+/*********************/
+
+#include "colors_waves.c"
+
+#define CLUSTER_SHIFT 10    /* shift in numbering of open clusters */
+
+int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression)
+{
+  TIFF *file;
+  GLubyte *image, *p;
+  int i;
+
+  file = TIFFOpen(filename, "w");
+  if (file == NULL)
+  {
+    return 1;
+  }
+
+  image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3);
+
+  /* OpenGL's default 4 byte pack alignment would leave extra bytes at the
+     end of each image row so that each full row contained a number of bytes
+     divisible by 4.  Ie, an RGB row with 3 pixels and 8-bit componets would
+     be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist
+     just to pad the row out to 12 bytes (12 is divisible by 4). To make sure
+     the rows are packed as tight as possible (no row padding), set the pack
+     alignment to 1. */
+
+  glPixelStorei(GL_PACK_ALIGNMENT, 1);
+
+  glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
+  TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
+  TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
+  TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT);
+  TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3);
+  TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
+  TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1);
+  TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description);
+  p = image;
+  for (i = height - 1; i >= 0; i--)
+  {
+//     if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0)
+    if (TIFFWriteScanline(file, p, i, 0) < 0)
+    {
+      free(image);
+      TIFFClose(file);
+      return 1;
+    }
+    p += width * sizeof(GLubyte) * 3;
+  }
+  TIFFClose(file);
+  return 0;
+}
+
+void init()		/* initialisation of window */
+{
+    glLineWidth(3);
+
+    glClearColor(0.0, 0.0, 0.0, 1.0);
+    glClear(GL_COLOR_BUFFER_BIT);
+
+//     glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0);
+    glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0);
+}
+
+void blank()
+{
+    if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0);
+    else glClearColor(1.0, 1.0, 1.0, 1.0);
+    glClear(GL_COLOR_BUFFER_BIT);
+}
+
+
+void write_text_fixedwidth( double x, double y, char *st)
+{
+    int l, i;
+
+    l=strlen( st ); // see how many characters are in text string.
+    glRasterPos2d( x, y); // location to start printing text
+    for( i=0; i < l; i++) // loop until i is greater then l
+    {
+//         glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen
+//    glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen
+   glutBitmapCharacter(GLUT_BITMAP_9_BY_15, st[i]); // Print a character on the screen
+    }
+} 
+
+void write_text( double x, double y, char *st)
+{
+    int l,i;
+
+    l=strlen( st ); // see how many characters are in text string.
+    glRasterPos2d( x, y); // location to start printing text
+    for( i=0; i < l; i++) // loop until i is greater then l
+    {
+        glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen
+//    glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen
+    }
+} 
+
+
+void save_frame_perc()
+{
+  static int counter = 0;
+  char *name="perc.", n2[100];
+  char format[6]=".%05i";
+
+    counter++;
+    strcpy(n2, name);
+    sprintf(strstr(n2,"."), format, counter);
+    strcat(n2, ".tif");
+    printf(" saving frame %s \n",n2);
+    writetiff(n2, "Wave equation in a planar domain", 0, 0,
+         WINWIDTH, WINHEIGHT, COMPRESSION_LZW);
+
+}
+
+/*********************/
+/* some basic math   */
+/*********************/
+
+
+double ipow(double x, int n)
+ {
+    double y;
+    int i;
+    
+    y = x;
+    for (i=1; i<n; i++) y *= x;
+    
+    return(y);
+ }
+ 
+
+int ipowi(int base, int n)
+ {
+    int p, i;
+    
+    if (n == 0) return(1);
+    else
+    {
+        p = base;
+        for (i=1; i<n; i++) p *= base;
+    
+        return(p);
+    }
+ }
+ 
+
+
+ 
+/*********************/
+/* drawing routines  */
+/*********************/
+
+/* The billiard boundary is drawn in (x,y) coordinates               */
+/* However for the grid points, we use integer coordinates (i,j)     */
+/* GL would allow to always work in (x,y) coordinates but using both */
+/* sets of coordinates decreases number of double computations when  */
+/* drawing the field                                                 */
+
+void xy_to_ij(double x, double y, int ij[2])
+/* convert (x,y) position to (i,j) in table representing wave */
+{
+    double x1, y1;
+
+    x1 = (x - XMIN)/(XMAX - XMIN);
+    y1 = (y - YMIN)/(YMAX - YMIN);
+
+    ij[0] = (int)(x1 * (double)NX);
+    ij[1] = (int)(y1 * (double)NY);
+}
+
+
+void xy_to_pos(double x, double y, double pos[2])
+/* convert (x,y) position to double-valued position in table representing wave */
+{
+    double x1, y1;
+
+    x1 = (x - XMIN)/(XMAX - XMIN);
+    y1 = (y - YMIN)/(YMAX - YMIN);
+
+    pos[0] = x1 * (double)NX;
+    pos[1] = y1 * (double)NY;
+}
+
+void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3])
+{
+    double pos[2];
+    
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+    glBegin(GL_QUADS);
+    xy_to_pos(x - dx, y - dy, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x + dx, y - dy, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x + dx, y + dy, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x - dx, y + dy, pos);
+    glVertex2d(pos[0], pos[1]);
+    glEnd();
+}
+
+
+void draw_line(double x1, double y1, double x2, double y2)
+{
+    double pos[2];
+    
+    glBegin(GL_LINE_STRIP);
+    xy_to_pos(x1, y1, pos);
+    glVertex2d(pos[0], pos[1]);
+    xy_to_pos(x2, y2, pos);
+    glVertex2d(pos[0], pos[1]);
+    glEnd();    
+}
+
+void draw_rectangle(double x1, double y1, double x2, double 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();    
+}
+
+void draw_colored_rectangle(double x1, double y1, double x2, double y2, double hue)
+{
+    double pos[2], rgb[3];
+    
+    hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE);
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+    
+    glBegin(GL_QUADS);
+    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();    
+
+    glColor3f(0.0, 0.0, 0.0);
+    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();    
+
+}
+
+int graphical_rep(int bcondition)
+/* return type of drawing, depending on lattice/boundary conditions */
+{
+    switch (bcondition) {
+        case (BC_SQUARE_DIRICHLET): return(PLOT_SQUARES);
+        case (BC_SQUARE_PERIODIC): return(PLOT_SQUARES);
+        case (BC_SQUARE_BOND_DIRICHLET): return(PLOT_SQUARE_BONDS);
+        case (BC_HEX_SITE_DIRICHLET): return(PLOT_HEX); 
+        case (BC_HEX_BOND_DIRICHLET): return(PLOT_HEX_BONDS); 
+        case (BC_TRIANGLE_SITE_DIRICHLET): return(PLOT_TRIANGLE); 
+        default: return(0);
+    }
+    
+}
+
+double pcritical(int bcondition)
+/* critical probability in terms of lattice and boundary condition */
+{
+    switch (LATTICE) {
+        case (BC_SQUARE_DIRICHLET): return(0.59274);
+        case (BC_SQUARE_PERIODIC): return(0.59274);
+        case (BC_SQUARE_BOND_DIRICHLET): return(0.5);
+        case (BC_HEX_BOND_DIRICHLET): return(1.0 - 2.0*sin(PI/18.0));
+        case (BC_TRIANGLE_SITE_DIRICHLET): return(0.6970402);
+        default: return(0.5);
+        
+    }
+    
+}
+
+int cellnb(int i, int j, int group, int nx, int ny)
+/* convert 2d coordinates to 1d */
+{
+    switch (LATTICE) {
+        case (BC_SQUARE_DIRICHLET):
+        {
+            return(i*ny+j);
+            break;
+        }
+        case (BC_SQUARE_PERIODIC):
+        {
+            return(i*ny+j);
+            break;
+        }
+        case (BC_SQUARE_BOND_DIRICHLET):
+        {
+            if (group == 0) return(i+nx*j);
+            else return (nx*(ny+1) + i*ny+j);
+            break;
+        }
+        case (BC_HEX_SITE_DIRICHLET):
+        {
+            return(i+nx*j);
+            break;
+        }
+        case (BC_HEX_BOND_DIRICHLET):
+        {
+            return (group*nx*(ny+1) + i+nx*j);
+            break;
+        }
+        case (BC_TRIANGLE_SITE_DIRICHLET):
+        {
+            return(i+2*nx*j);
+            break;
+        }
+    }
+}
+
+int cell_to_ij(int c, int *i, int *j, int nx, int ny)
+/* convert 1d coordinates to 2d, returns group */
+{
+    int group;
+    
+    switch (LATTICE) {
+        case (BC_SQUARE_DIRICHLET):
+        {
+            *i = c/ny;
+            *j = c - *i*ny;
+            return(0);
+        }
+        case (BC_SQUARE_PERIODIC):
+        {
+            *i = c/ny;
+            *j = c - *i*ny;
+            return(0);
+        }
+        case (BC_SQUARE_BOND_DIRICHLET):
+        {
+            if (c < nx*(ny+1)) 
+            {
+                *j = c/nx;
+                *i = c - *j*nx;
+                return(0);
+            }
+            else 
+            {
+                c -= nx*(ny+1);
+                *i = c/ny;
+                *j = c - *i*ny;
+                return(1);
+            }
+        }
+        case (BC_HEX_SITE_DIRICHLET):
+        {
+            *j = c/nx;
+            *i = c - *j*nx;
+            return(0);
+        }
+        case (BC_HEX_BOND_DIRICHLET):
+        {
+            group = c/(nx*(ny+1));
+            c -= group*(nx*(ny+1));
+            *j = c/nx;
+            *i = c - *j*nx;
+            return(group);
+        }
+        case (BC_TRIANGLE_SITE_DIRICHLET):
+        {
+            *j= c/(2*nx);
+            *i = c - 2*(*j)*nx;
+            return(0);
+        }
+    }
+}
+
+double p_schedule(int i)
+/* percolation probability p as a function of time */
+{
+    double time, pstar;
+    int power = 2, factor;
+    
+    pstar = pcritical(LATTICE);
+    
+    factor = ipow(2, power);
+    
+    time = (double)i/(double)(NSTEPS-1);
+    
+    if (time > 0.5) return(pstar + factor*(1.0 - pstar)*ipow(time - 0.5, power));
+    else return(pstar - factor*pstar*ipow(0.5 - time, power));
+}
+
+int in_plot_box(double x, double y)
+{
+    int pos[2];
+    static double xmin, xmax, ymin, ymax;
+    static int first = 1;
+    
+    if (first)
+    {
+        xy_to_ij(XMAX - 1.0, YMAX - 1.0, pos);
+        xmin = (double)pos[0];
+        ymin = (double)pos[1];
+        first = 0;
+    }
+    
+    return((x > xmin)&&(y > ymin));
+//     return((x + 0.5*(double)size > xmin)&&(y + 0.5*(double)size > ymin));
+}
+
+double size_ratio_color(int clustersize, int ncells)
+/* color of cell as function of the size of its cluster */
+{
+    double ratio, minratio = 1.0e-2, x;
+    
+    ratio = (double)clustersize/(double)ncells;
+    if (ratio > 1.0) ratio = 1.0;
+    else if (ratio < minratio) ratio = minratio;
+    x = log(ratio/minratio)/log(1.0/minratio);
+    return(CLUSTER_HUEMIN*x + CLUSTER_HUEMAX*(1.0 -x));
+    
+    /* other attempts that seem to bug */
+//     ratio = log((double)(clustersize+1))/log((double)(ncells+1));
+// //     ratio = ((double)clustersize)/((double)ncells);
+// //     ratio = sqrt((double)clustersize)/sqrt((double)ncells);
+//     if (ratio > 1.0) ratio = 1.0;
+//     else if (ratio < 0.0) ratio = 0.0;
+//     return(CLUSTER_HUEMAX*ratio + CLUSTER_HUEMIN*(1.0 -ratio));
+}
+
+void set_cell_color(t_perco cell, int *cluster_sizes, int fade, int max_cluster_size)
+/* compute color of cell */
+{
+    int k, color, csize;
+    double rgb[3], fade_factor = 0.15, hue;
+    
+    if (!cell.open) hsl_to_rgb_palette(HUE_CLOSED, 0.9, 0.5, rgb, COLOR_PALETTE);
+    else
+        {
+            if (!cell.flooded) hsl_to_rgb_palette(HUE_OPEN, 0.9, 0.5, rgb, COLOR_PALETTE);
+            else hsl_to_rgb_palette(HUE_FLOODED, 0.9, 0.5, rgb, COLOR_PALETTE);
+            
+            if ((FIND_ALL_CLUSTERS)&&(COLOR_CLUSTERS_BY_SIZE))
+            {
+                csize = cluster_sizes[cell.cluster]; 
+                hue = size_ratio_color(csize, max_cluster_size);
+                hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE);
+            }
+            else if ((FIND_ALL_CLUSTERS)&&(cell.cluster > 1))
+            {
+                color = cell.cluster%N_CLUSTER_COLORS;
+                hue = CLUSTER_HUEMIN + (CLUSTER_HUEMAX - CLUSTER_HUEMIN)*(double)color/(double)N_CLUSTER_COLORS;
+                hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE);
+            }
+        }
+        
+    if (fade) for (k=0; k<3; k++) rgb[k] = 1.0 - fade_factor + fade_factor*rgb[k];
+        
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+}
+
+
+double plot_coord(double x, double xmin, double xmax)
+{
+    return(xmin + x*(xmax - xmin));
+}
+
+void draw_size_plot(double plot_cluster_size[NSTEPS], int i, double pcrit)
+/* draw plot of cluster sizes in terms of p */
+{
+    int j;
+    char message[100];
+    static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax;
+    double pos[2], x1, y1, x2, y2, rgb[3], x;
+    static int first = 1;
+    
+    if (first)
+    {
+        xmin = XMAX - 0.95;
+        xmax = XMAX - 0.05;
+        ymin = YMAX - 0.95;
+        ymax = YMAX - 0.05;
+                
+        xmid = 0.5*(xmin + xmax);
+        ymid = 0.5*(ymin + ymax);
+        
+        dx = 0.5*(xmax - xmin);
+        dy = 0.5*(ymax - ymin);
+        
+        plotxmin = xmin + 0.05;
+        plotxmax = xmax - 0.1;
+        plotymin = ymin + 0.07;
+        plotymax = ymax - 0.15;
+        
+        first = 0;
+    }
+    
+    rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0;
+//     erase_area_rgb(xmid, ymid, dx, dy, rgb);
+    
+    glColor3f(0.0, 0.0, 0.0);
+    glLineWidth(2);
+    
+    /* axes and labels */
+    draw_line(plotxmin, plotymin, plotxmax + 0.05, plotymin);
+    draw_line(plotxmin, plotymin, plotxmin, plotymax + 0.1);
+    draw_line(plotxmin - 0.02, plotymax, plotxmin + 0.02, plotymax);
+    x = plot_coord(pcrit, plotxmin, plotxmax);
+    draw_line(x, plotymin, x, plotymax);
+    draw_line(plotxmax, plotymin - 0.02, plotxmax, plotymin + 0.02);
+    
+    xy_to_pos(plotxmax + 0.06, plotymin - 0.03, pos);
+    sprintf(message, "p");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmax - 0.015, plotymin - 0.06, pos);
+    sprintf(message, "1");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(x - 0.02, plotymin - 0.04, pos);
+    sprintf(message, "pc");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmin + 0.02, plotymax + 0.05, pos);
+    sprintf(message, "nflooded/nopen");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmin - 0.05, plotymax - 0.01, pos);
+    sprintf(message, "1");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    /* plot */
+    hsl_to_rgb_palette(HUE_GRAPH, 0.9, 0.5, rgb, COLOR_PALETTE);
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+    x1 = plotxmin;
+    y1 = plotymin;
+    for (j=0; j<i; j++)
+    {
+//         x2 = plot_coord((double)j/(double)NSTEPS, plotxmin, plotxmax);
+        x2 = plot_coord(p_schedule(j), plotxmin, plotxmax);
+        y2 = plot_coord(plot_cluster_size[j], plotymin, plotymax);
+        
+        draw_line(x1, y1, x2, y2);
+        x1 = x2;
+        y1 = y2;
+    }
+}
+
+void draw_cluster_number_plot(int plot_cluster_number[NSTEPS], int max_number, int i, double pcrit)
+/* draw plot of number of clusters in terms of p */
+{
+    int j;
+    char message[100];
+    static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax;
+    double pos[2], x1, y1, x2, y2, rgb[3], x, y;
+    static int first = 1;
+    
+    if (first)
+    {
+        xmin = XMAX - 0.95;
+        xmax = XMAX - 0.05;
+        ymin = YMAX - 0.95;
+        ymax = YMAX - 0.05;
+                
+        xmid = 0.5*(xmin + xmax);
+        ymid = 0.5*(ymin + ymax);
+        
+        dx = 0.5*(xmax - xmin);
+        dy = 0.5*(ymax - ymin);
+        
+        plotxmin = xmin + 0.05;
+        plotxmax = xmax - 0.1;
+        plotymin = ymin + 0.07;
+        plotymax = ymax - 0.15;
+        
+        first = 0;
+    }
+    
+    rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0;
+//     erase_area_rgb(xmid, ymid, dx, dy, rgb);
+    
+    glColor3f(0.0, 0.0, 0.0);
+    glLineWidth(2);
+    
+    /* axes and labels */
+    draw_line(plotxmin, plotymin, plotxmax + 0.05, plotymin);
+    draw_line(plotxmin, plotymin, plotxmin, plotymax + 0.1);
+    draw_line(plotxmin - 0.02, plotymax, plotxmin + 0.02, plotymax);
+    x = plot_coord(pcrit, plotxmin, plotxmax);
+    draw_line(x, plotymin, x, plotymax);
+    draw_line(plotxmax, plotymin - 0.02, plotxmax, plotymin + 0.02);
+    
+    xy_to_pos(plotxmax + 0.06, plotymin - 0.03, pos);
+    sprintf(message, "p");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmax - 0.015, plotymin - 0.06, pos);
+    sprintf(message, "1");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(x - 0.02, plotymin - 0.04, pos);
+    sprintf(message, "pc");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmin + 0.02, plotymax + 0.05, pos);
+    sprintf(message, "clusters/cells");
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    xy_to_pos(plotxmin + 0.05, plotymax - 0.01, pos);
+    sprintf(message, "%.2f", 1.0/(double)MAX_CLUSTER_NUMBER);
+    write_text_fixedwidth(pos[0], pos[1], message);
+    
+    /* plot */
+    hsl_to_rgb_palette(HUE_GRAPH, 0.9, 0.5, rgb, COLOR_PALETTE);
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+    x1 = plotxmin;
+    y1 = plotymin;
+    for (j=0; j<i; j++)
+    {
+//         x2 = plot_coord((double)j/(double)NSTEPS, plotxmin, plotxmax);
+        x2 = plot_coord(p_schedule(j), plotxmin, plotxmax);
+        y2 = plot_coord((double)plot_cluster_number[j]/(double)max_number, plotymin, plotymax);
+        
+        draw_line(x1, y1, x2, y2);
+        x1 = x2;
+        y1 = y2;
+    }
+}
+
+void draw_cluster_histogram(int ncells, int *cluster_sizes, int maxclustersize, int maxclusterlabel)
+/* draw histogram of cluster size distribution */
+{   
+    int i, bin, nbins, binwidth, *histo, maxheight = 0, csize, dh, n, max_x_axis, di;
+    static int maxbins = HISTO_BINS;
+    char message[100];
+    static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax, x, y;
+    double pos[2], x1, y1, x2, y2, hue, delta;
+    static int first = 1;
+    
+    if (first)
+    {
+        xmin = XMAX - 0.95;
+        xmax = XMAX - 0.05;
+        ymin = YMAX - 0.95;
+        ymax = YMAX - 0.05;
+                
+        xmid = 0.5*(xmin + xmax);
+        ymid = 0.5*(ymin + ymax);
+        
+        dx = 0.5*(xmax - xmin);
+        dy = 0.5*(ymax - ymin);
+        
+        plotxmin = xmin + 0.15;
+        plotxmax = xmax - 0.1;
+        plotymin = ymin + 0.07;
+        plotymax = ymax - 0.1;
+        
+        first = 0;
+    }
+        
+    if (maxclustersize > 0)
+    {
+    
+//     if (maxclustersize < maxbins) nbins = maxclustersize;
+//     else 
+        nbins = maxbins;
+    
+        histo = (int *)malloc(nbins*sizeof(int));
+    
+        for (bin = 0; bin < nbins; bin++) histo[bin] = 0;
+    
+//         binwidth = maxclustersize/nbins;
+//         if (binwidth == 0) binwidth = 1;
+        binwidth = maxclustersize/nbins + 1;
+    
+        /* compute histogram */
+        for (i=CLUSTER_SHIFT; i<maxclusterlabel; i++) if (cluster_sizes[i] > 0)
+        {
+            bin = (cluster_sizes[i]-1)/binwidth;
+            histo[bin]++;
+        }
+        for (bin=0; bin<maxbins; bin++) if (histo[bin] > maxheight) maxheight = histo[bin];
+        
+        /* draw histogram */
+        glColor3f(0.0, 0.0, 0.0);
+        glLineWidth(2);
+    
+        x1 = plotxmin;
+        y1 = plotymin;
+        max_x_axis = maxclustersize;
+        if (max_x_axis < HISTO_BINS) max_x_axis = HISTO_BINS;
+        
+        for (bin=0; bin < nbins; bin++)
+        {
+            csize = bin*binwidth + binwidth/2;
+            hue = size_ratio_color(csize, ncells);
+            
+            x2 = plot_coord((double)((bin+1)*binwidth)/(double)(max_x_axis+1), plotxmin, plotxmax);
+//             x2 = plot_coord((double)(bin+1)/(double)nbins, plotxmin, plotxmax);
+            
+            y = log((double)(histo[bin]+1))/log((double)(maxheight+1));
+            y2 = plot_coord(y, plotymin, plotymax);
+        
+            draw_colored_rectangle(x1, y1, x2, y2, hue);
+            x1 = x2;
+        }
+        
+        draw_line(plotxmin, plotymin, plotxmax + 0.05, plotymin);
+        draw_line(plotxmin, plotymin, plotxmin, plotymax + 0.1);
+        
+        /* graduation of y axis */
+        x = log((double)(maxheight+1))/log(10.0);
+        for (i=0; i<(int)x + 1; i++)
+        {
+            n = ipowi(10, i);
+            y = log((double)(n+1))/log((double)(maxheight+1));
+            y1 = plot_coord(y, plotymin, plotymax);
+            xy_to_pos(plotxmin - 0.1, y1, pos);
+            draw_line(plotxmin - 0.02, y1, plotxmin + 0.02, y1);
+            sprintf(message, "%i", n);
+            xy_to_pos(plotxmin - 0.07 - 0.025*(double)i, y1 - 0.01, pos);
+            write_text_fixedwidth(pos[0], pos[1], message);
+        }
+    
+        /* graduation of x axis */
+        x = log((double)(max_x_axis+1))/log(10.0);
+        n = ipowi(10, (int)x);
+        y = (double)n/10.0;
+
+        delta = plot_coord((double)n/((double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin;
+        if (delta > 0.13 + 0.01*x) di = 1;
+        else if (delta > 0.08 + 0.01*x) di = 2;
+        else di = 5;
+        for (i=di; i<10; i+=di)
+        {
+            x1 = plot_coord(y*(double)i*10.0/(double)(max_x_axis+1), plotxmin, plotxmax);
+            if (i*n < max_x_axis*11/10)
+            {
+                sprintf(message, "%i", i*n);
+                xy_to_pos(x1 + 0.005 - 0.012*x, plotymin - 0.07, pos);
+                write_text_fixedwidth(pos[0], pos[1], message);
+            }
+        }
+        
+        delta = plot_coord((double)n/(10.0*(double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin;
+        for (i=0; i<100; i++) if (i*n < 11*max_x_axis)
+        {
+            y = (double)(i*n)/10.0;
+            x1 = plot_coord(y/(double)(max_x_axis+1), plotxmin, plotxmax);
+            xy_to_pos(x1, plotymin , pos);
+            if (i%10 == 0) draw_line(x1, plotymin - 0.02, x1, plotymin + 0.02);
+            else if (delta > 0.02) draw_line(x1, plotymin - 0.01, x1, plotymin + 0.01);
+        }
+        
+        /* for debugging */
+        if (DEBUG) for (bin=0; bin < nbins; bin++) printf("Bin %i = %i\n", bin, histo[bin]);
+        
+        
+        free(histo);
+    }
+}
+
+void draw_cell(int c, int nx, int ny, int size)
+/* draw a single cell, for debugging puposes */
+{
+    int i, j, ishift, k, group;
+    double rgb[3], x, y, alpha, r, fade = 0, dsize;
+    static double h;
+    static int first = 1;
+    
+    dsize = (double)size;
+    
+    group = cell_to_ij(c, &i, &j, nx, ny);
+    switch (graphical_rep(LATTICE)) {
+        case (PLOT_SQUARES):
+        {
+            glBegin(GL_QUADS);
+    
+            glVertex2i(i*size, j*size);
+            glVertex2i((i+1)*size, j*size);
+            glVertex2i((i+1)*size, (j+1)*size);
+            glVertex2i(i*size, (j+1)*size);
+
+            glEnd ();
+            
+            break;
+        }
+        case (PLOT_SQUARE_BONDS):
+        {
+            glLineWidth(5);
+            glBegin(GL_LINES);
+            
+            /* horizontal segments */
+            if (group == 0)
+            {
+                glVertex2i(i*size, j*size);
+                glVertex2i((i+1)*size, j*size);
+            }
+            
+            /* vertical segments */
+            else 
+            {
+                glVertex2i(i*size, j*size);
+                glVertex2i(i*size, (j+1)*size);
+            }
+                
+            glEnd ();
+            break;
+        }
+        case (PLOT_HEX):
+        {
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+    
+            x = ((double)i + 0.75)*dsize;
+            if (j%2 == 1) x += 0.5*dsize;
+            y = h*dsize*((double)j + 1.0);
+            
+            glBegin(GL_TRIANGLE_FAN);
+            glVertex2d(x, y);
+            for (k=0; k<7; k++)
+            {
+                alpha = (1.0 + 2.0*(double)k)*PI/6.0;
+                glVertex2d(x + r*cos(alpha), y + r*sin(alpha));
+            }
+            glEnd();
+            break;
+        }
+        case (PLOT_HEX_BONDS):
+        {
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+            
+            x = ((double)i + 0.75)*dsize;
+            if (j%2 == 1) x -= 0.5*dsize;
+            y = h*dsize*((double)j + 1.0);
+            
+            glBegin(GL_LINES);
+            switch (group){
+                case (0):   /* vertical bonds */
+                {
+                    glVertex2d(x-0.5*dsize, y+0.5*r);
+                    glVertex2d(x-0.5*dsize, y-0.5*r);
+                    break;
+                }
+                case (1):   /* NE-SW bonds */
+                {
+                    glVertex2d(x, y-r);
+                    glVertex2d(x+0.5*dsize, y-0.5*r);
+                    break;
+                }
+                case (2):   /* NW-SE bonds */
+                {
+                    glVertex2d(x-0.5*dsize, y-0.5*r);
+                    glVertex2d(x, y-r);
+                    break;
+                }
+            }
+
+            glEnd();
+            break;
+        }
+        case (PLOT_TRIANGLE):
+        {
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+    
+            x = 0.5*((double)i + 1.25)*dsize;
+            y = h*dsize*((double)j);
+            
+            printf("Drawing cell %i = (%i, %i) at (%.0f, %.0f)\n", c, i, j, x, y);
+                    
+            glBegin(GL_TRIANGLES);
+            if ((i+j)%2 == 1)
+            {
+                glVertex2d(x-0.5*dsize, y);
+                glVertex2d(x+0.5*dsize, y);
+                glVertex2d(x, y+h*dsize);
+            }
+            else
+            {
+                glVertex2d(x-0.5*dsize, y+h*dsize);
+                glVertex2d(x+0.5*dsize, y+h*dsize);
+                glVertex2d(x, y);
+            }
+            glEnd();
+            break;
+        }
+    }
+}
+
+void test_neighbours(int start, t_perco *cell, int nx, int ny, int size, int ncells)
+/* for debugging puposes */
+{
+    int i, k;
+    
+    for (i=start; i<ncells; i++)
+    {
+        printf("Testing cell %i of %i - %i neighbours\n", i, ncells, cell[i].nneighb);
+        blank();
+        glColor3f(1.0, 0.0, 0.0);
+        draw_cell(i, nx, ny, size);
+        glColor3f(0.0, 0.0, 1.0);
+        for (k=0; k<cell[i].nneighb; k++) draw_cell(cell[i].nghb[k], nx, ny, size);
+        glutSwapBuffers();
+        sleep(1);
+    }
+}   
+
+
+
+void draw_configuration(t_perco *cell, int *cluster_sizes, int nx, int ny, int size, int max_cluster_size)
+/* draw the configuration */
+{
+    int i, j, ishift, k;
+    double rgb[3], x, y, alpha, r, fade = 0, dsize;
+    static double h, h1;
+    static int first = 1;
+    
+    dsize = (double)size;
+    
+    switch (graphical_rep(LATTICE)) {
+        case (PLOT_SQUARES):
+        {
+            blank();
+            
+            glBegin(GL_QUADS);
+    
+            for (i=0; i<nx; i++)
+                for (j=0; j<ny; j++)
+                {
+                    if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j)*dsize))) fade = 1;
+                    else fade = 0;
+                    
+                    set_cell_color(cell[i*ny+j], cluster_sizes, fade, max_cluster_size);
+            
+                    glVertex2i(i*size, j*size);
+                    glVertex2i((i+1)*size, j*size);
+                    glVertex2i((i+1)*size, (j+1)*size);
+                    glVertex2i(i*size, (j+1)*size);
+                }
+
+            glEnd ();
+            break;
+        }
+        
+        case (PLOT_SQUARE_BONDS):
+        {
+            ishift = nx*(ny+1);
+            
+            blank();
+       
+            if (ADD_PLOT)
+            {
+                rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0;
+                erase_area_rgb(XMAX - 0.5, YMAX - 0.5, 0.5, 0.5, rgb);
+            }
+            
+            if (size < 8) glLineWidth(1 + size/4);
+            else glLineWidth(3);
+            
+            glBegin(GL_LINES);
+            
+            /* horizontal segments */
+            for (i=0; i<nx; i++)
+                for (j=0; j<ny+1; j++)
+                {
+                    if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j)*dsize))) fade = 1;
+                    else fade = 0;
+                    
+                    set_cell_color(cell[i+nx*j], cluster_sizes, fade, max_cluster_size);
+                    
+                    glVertex2i(i*size, j*size);
+                    glVertex2i((i+1)*size, j*size);
+                }
+                
+            /* vertical segments */
+            for (i=0; i<nx+1; i++)
+                for (j=0; j<ny; j++)
+                {
+                    if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j+1)*dsize))) fade = 1;
+                    else fade = 0;
+                    
+                    set_cell_color(cell[ishift+ny*i+j], cluster_sizes, fade, max_cluster_size);
+                    
+                    glVertex2i(i*size, j*size);
+                    glVertex2i(i*size, (j+1)*size);
+                }
+                
+            glEnd ();
+            break;
+        }
+        
+        case (PLOT_HEX):
+        {
+            blank();
+            
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+    
+            for (j=0; j<ny; j++)
+                for (i=0; i<nx; i++)
+                {
+                    
+                    x = ((double)i + 0.75)*dsize;
+                    if (j%2 == 1) x += 0.5*dsize;
+                    y = h*dsize*((double)j + 1.0);
+                    
+                    if ((ADD_PLOT)&&(in_plot_box(x, y))) fade = 1;
+                    else fade = 0;
+                    
+                    set_cell_color(cell[j*nx+i], cluster_sizes, fade, max_cluster_size);
+                    
+                    glBegin(GL_TRIANGLE_FAN);
+                    
+                    glVertex2d(x, y);
+                    for (k=0; k<7; k++)
+                    {
+                        alpha = (1.0 + 2.0*(double)k)*PI/6.0;
+                        glVertex2d(x + r*cos(alpha), y + r*sin(alpha));
+                    }
+                    glEnd();
+                }
+            break;
+        }
+        case (PLOT_HEX_BONDS):
+        {
+            blank();
+            
+            ishift = nx*(ny+1);
+            
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+    
+            if (ADD_PLOT)
+            {
+                rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0;
+                erase_area_rgb(XMAX - 0.5, YMAX - 0.5, 0.5, 0.5, rgb);
+            }
+            
+            if (size < 8) glLineWidth(1 + size/4);
+            else glLineWidth(3);
+            
+            glBegin(GL_LINES);
+            for (i=0; i<nx; i++)
+                for (j=0; j<ny; j++)
+                {
+                    x = ((double)i + 0.5)*dsize;
+                    if (j%2 == 1) x -= 0.5*dsize;
+                    y = h*dsize*((double)j + 0.75);
+                    
+                    if ((ADD_PLOT)&&(in_plot_box(x, y))) fade = 1;
+                    else fade = 0;
+                    
+                    /* vertical bonds */
+                    if (j<ny-1)
+                    {
+                        set_cell_color(cell[i+nx*j], cluster_sizes, fade, max_cluster_size);
+                    
+                        glVertex2d(x-0.5*dsize, y+0.5*r);
+                        glVertex2d(x-0.5*dsize, y-0.5*r);
+                    }
+                    
+                    if ((ADD_PLOT)&&(in_plot_box(x, y-0.5*h*dsize))) fade = 1;
+                    else fade = 0;
+                    
+                    /* NW-SE bonds */
+                    set_cell_color(cell[2*ishift + i+nx*j], cluster_sizes, fade, max_cluster_size);
+                    
+                    glVertex2d(x-0.5*dsize, y-0.5*r);
+                    glVertex2d(x, y-r);
+                    
+                    if ((ADD_PLOT)&&(in_plot_box(x+0.5*dsize, y-0.5*h*dsize))) fade = 1;
+                    else fade = 0;
+                    
+                    /* NE-SW bonds */
+                    set_cell_color(cell[ishift + i+nx*j], cluster_sizes, fade, max_cluster_size);
+                    
+                    glVertex2d(x, y-r);
+                    glVertex2d(x+0.5*dsize, y-0.5*r);
+
+                    
+                }
+           
+            glEnd();
+            break;
+        }
+        case (PLOT_TRIANGLE):
+        {
+            blank();
+            
+            if (first)
+            {
+                h = 0.5*sqrt(3.0);
+                h1 = 0.5/sqrt(3.0);
+                first = 0;
+            }
+            r = (double)size*0.5/h;
+    
+            for (j=0; j<ny; j++)
+                for (i=0; i<2*nx; i++)
+                {
+                    
+                    x = 0.5*((double)i + 1.25)*dsize;
+                    y = h*dsize*((double)j);
+                    
+                    if ((ADD_PLOT)&&(in_plot_box(x, y))) fade = 1;
+                    else fade = 0;
+                    
+                    set_cell_color(cell[2*j*nx+i], cluster_sizes, fade, max_cluster_size);
+                    
+                    glBegin(GL_TRIANGLES);
+                    if ((i+j)%2 == 1)
+                    {
+                        glVertex2d(x-0.5*dsize, y);
+                        glVertex2d(x+0.5*dsize, y);
+                        glVertex2d(x, y+h*dsize);
+                    }
+                    else
+                    {
+                        glVertex2d(x-0.5*dsize, y+h*dsize);
+                        glVertex2d(x+0.5*dsize, y+h*dsize);
+                        glVertex2d(x, y);
+                    }
+                    glEnd();
+                }
+            break;
+        }
+    }
+}
+
+
+void print_p(double p)
+{
+    char message[100];
+    double y = YMAX - 0.13, pos[2], rgb[3] = {1.0, 1.0, 1.0};
+    static double xleftbox, xlefttext, xrightbox, xrighttext;
+    static int first = 1;
+    
+    if (first)
+    {
+        xleftbox = XMIN + 0.2;
+        xlefttext = xleftbox - 0.45;
+        if (PLOT_CLUSTER_HISTOGRAM) xrightbox = XMAX - 0.41;
+        else xrightbox = XMAX - 0.27;
+        xrighttext = xrightbox - 0.45;
+        first = 0;
+    }
+    
+    if (!PLOT_CLUSTER_SIZE)
+    {
+        if (NX > 1280) erase_area_rgb(xrightbox, y + 0.025, 0.15, 0.05, rgb);
+        else erase_area_rgb(xrightbox, y + 0.025, 0.22, 0.05, rgb);
+    }
+    glColor3f(0.0, 0.0, 0.0);
+    if (NX > 1280) xy_to_pos(xrighttext + 0.35, y, pos);
+    else xy_to_pos(xrighttext + 0.3, y, pos);
+    sprintf(message, "p = %.4f", p);
+    write_text(pos[0], pos[1], message);
+}
+
+void print_nclusters(int nclusters)
+{
+    char message[100];
+    double y = YMAX - 0.25, pos[2], rgb[3] = {1.0, 1.0, 1.0};
+    static double xleftbox, xlefttext, xrightbox, xrighttext;
+    static int first = 1;
+    
+    if (first)
+    {
+        xleftbox = XMIN + 0.2;
+        xlefttext = xleftbox - 0.45;
+        xrightbox = XMAX - 0.31;
+        xrighttext = xrightbox - 0.48;
+        first = 0;
+    }
+    
+    if (!PLOT_CLUSTER_SIZE)
+    {
+        if (NX > 1280) erase_area_rgb(xrightbox, y + 0.025, 0.18, 0.05, rgb);
+        else erase_area_rgb(xrightbox, y + 0.025, 0.25, 0.05, rgb);
+    }
+    glColor3f(0.0, 0.0, 0.0);
+    if (NX > 1280) xy_to_pos(xrighttext + 0.35, y, pos);
+    else xy_to_pos(xrighttext + 0.3, y, pos);
+    if (nclusters == 1) sprintf(message, "%i cluster", nclusters);
+    else sprintf(message, "%i clusters", nclusters);
+    write_text_fixedwidth(pos[0], pos[1], message);
+}
+
+void print_largest_cluster_size(int max_cluster_size)
+{
+    char message[100];
+    double y = YMAX - 0.25, pos[2], rgb[3] = {1.0, 1.0, 1.0};
+    static double xleftbox, xlefttext, xrightbox, xrighttext;
+    static int first = 1;
+    
+    if (first)
+    {
+        xleftbox = XMIN + 0.2;
+        xlefttext = xleftbox - 0.45;
+        xrightbox = XMAX - 0.41;
+        xrighttext = xrightbox - 0.48;
+        first = 0;
+    }
+    
+    if (!PLOT_CLUSTER_SIZE)
+    {
+        if (NX > 1280) erase_area_rgb(xrightbox, y + 0.025, 0.18, 0.05, rgb);
+        else erase_area_rgb(xrightbox, y + 0.025, 0.25, 0.05, rgb);
+    }
+    glColor3f(0.0, 0.0, 0.0);
+    if (NX > 1280) xy_to_pos(xrighttext + 0.35, y, pos);
+    else xy_to_pos(xrighttext + 0.3, y, pos);
+    sprintf(message, "max size %i", max_cluster_size);
+    write_text_fixedwidth(pos[0], pos[1], message);
+}
+
+
+/*********************/
+/* animation part    */
+/*********************/
+
+int cell_number(int nx, int ny)
+/* total number of cells in graph */
+{
+    switch (LATTICE) {
+        case (BC_SQUARE_DIRICHLET): return(nx*ny);
+        case (BC_SQUARE_PERIODIC): return(nx*ny);
+        case (BC_SQUARE_BOND_DIRICHLET): return(2*nx*ny + nx + ny);
+        case (BC_HEX_SITE_DIRICHLET): return((int)((double)(nx*ny)*2.0/sqrt(3.0))); /* hex lattice requires more vertical space! */
+        case (BC_HEX_BOND_DIRICHLET): return(3*(int)((double)((nx+2)*(ny+2))*2.0/sqrt(3.0))); 
+        /* hex lattice requires more vertical space! */
+        case (BC_TRIANGLE_SITE_DIRICHLET): return((int)((double)(2*nx*ny)*2.0/sqrt(3.0))); /* triangle lattice requires more vertical space! */
+    }
+}
+
+
+void compute_nxny(int size, int *nx, int *ny)
+/* compute the number of rows and columns depending on lattice */
+{
+    switch (LATTICE) {
+        case (BC_HEX_SITE_DIRICHLET):
+        {
+            *nx = NX/size - 1;
+            *ny = (int)((double)NY*2.0/((double)size*sqrt(3.0)));
+            break;
+        }
+        case (BC_HEX_BOND_DIRICHLET):
+        {
+            *nx = NX/size + 1;
+            *ny = (int)((double)NY*2.0/((double)size*sqrt(3.0))) + 1;
+            break;
+        }
+        case (BC_TRIANGLE_SITE_DIRICHLET):
+        {
+            *nx = NX/size - 1;
+            *ny = (int)((double)NY*2.0/((double)size*sqrt(3.0))) + 1;
+            break;
+        }
+        default:
+        {
+            *nx = NX/size;
+            *ny = NY/size;
+        }
+    }
+    
+}
+
+int init_cell_lattice(t_perco *cell, int nx, int ny)
+/* initialize the graph of connected cells - returns the number of cells */
+{
+    int i, j, iplus, iminus, ishift, n;
+    
+    printf("Initializing cell lattice ...");
+    
+    switch (LATTICE) {
+        case (BC_SQUARE_DIRICHLET):
+        {
+            /* neighbours in the bulk */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=1; j<ny-1; j++){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j+1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j-1, 0, nx, ny);
+                }
+            }
+    
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++)
+            {
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j-1, 0, nx, ny);
+            }
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++)
+            {
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx-2, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j-1, 0, nx, ny);
+            }
+            
+            /* top boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, ny-1, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i+1, ny-1, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, ny-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, ny-2, 0, nx, ny);
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, 0, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i+1, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, 0, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, 1, 0, nx, ny);
+            }
+            
+            /* corners */
+            n = cellnb(0, 0, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(1, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, 1, 0, nx, ny);
+            
+            n = cellnb(0, ny-1, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(0, ny-2, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(1, ny-1, 0, nx, ny);            
+            
+            n = cellnb(nx-1, 0, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(nx-1, 1, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-2, 0, 0, nx, ny);
+            
+            n = cellnb(nx-1, ny-1, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(nx-1, ny-2, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-2, ny-1, 0, nx, ny);
+            
+            return(nx*ny);
+            break;
+        }
+        case (BC_SQUARE_PERIODIC):
+        {
+            /* neighbours in the bulk */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=1; j<ny-1; j++){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j+1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j-1, 0, nx, ny);
+                }
+            }
+    
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++)
+            {
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j+1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(0, j+1, 0, nx, ny);
+            }
+            
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++)
+            {
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(nx-2, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j+1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(nx-1, j-1, 0, nx, ny);
+            }
+            
+            /* top boundary */
+            #pragma omp parallel for private(i,iplus,iminus)
+            for (i=0; i<nx; i++)
+            {
+                n = cellnb(i, ny-1, 0, nx, ny);
+                iplus = (i+1) % nx;
+                iminus = (i-1) % nx;
+                if (iminus < 0) iminus += nx;
+                    
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(iplus, ny-1, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(iminus, ny-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, ny-2, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i, 0, 0, nx, ny);
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i,iplus,iminus)
+            for (i=0; i<nx; i++)
+            {
+                n = cellnb(i, 0, 0, nx, ny);
+                iplus = (i+1) % nx;
+                iminus = (i-1) % nx;
+                if (iminus < 0) iminus += nx;
+                    
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(iplus, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(iminus, 0, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, 1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i, ny-1, 0, nx, ny);
+            }
+            
+            return(nx*ny);
+            break;
+        }
+        case (BC_SQUARE_BOND_DIRICHLET): 
+        {
+            ishift = nx*(ny+1); 
+            
+            /* horizontal bonds */
+            /* neighbours in the bulk  */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=1; j<ny; j++){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 6;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j-1, 1, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j, 1, nx, ny);
+                    cell[n].nghb[4] = cellnb(i+1, j-1, 1, nx, ny);
+                    cell[n].nghb[5] = cellnb(i+1, j, 1, nx, ny);
+                }
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, 0, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(i+1, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, 0, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, 0, 1, nx, ny);
+                cell[n].nghb[3] = cellnb(i+1, 0, 1, nx, ny);
+            }
+
+            /* top boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, ny, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(i+1, ny, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, ny, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, ny-1, 1, nx, ny);
+                cell[n].nghb[3] = cellnb(i+1, ny-1, 1, nx, ny);
+            }
+            
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny; j++){
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j-1, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j, 1, nx, ny);
+                cell[n].nghb[3] = cellnb(1, j-1, 1, nx, ny);
+                cell[n].nghb[4] = cellnb(1, j, 1, nx, ny);
+            }
+            
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny; j++){
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(nx-2, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j-1, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j, 1, nx, ny);
+                cell[n].nghb[3] = cellnb(nx, j-1, 1, nx, ny);
+                cell[n].nghb[4] = cellnb(nx, j, 1, nx, ny);
+            }
+            
+            /* corners */
+            n = cellnb(0, 0, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(1, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, 0, 1, nx, ny);
+            cell[n].nghb[2] = cellnb(1, 0, 1, nx, ny);
+            
+            n = cellnb(nx-1, 0, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx-2, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, 0, 1, nx, ny);
+            cell[n].nghb[2] = cellnb(nx, 0, 1, nx, ny);
+            
+            n = cellnb(0, ny, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(1, ny, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, ny-1, 1, nx, ny);
+            cell[n].nghb[2] = cellnb(1, ny-1, 1, nx, ny);
+            
+            n = cellnb(nx-1, ny, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx-2, ny, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, ny-1, 1, nx, ny);
+            cell[n].nghb[2] = cellnb(nx, ny-1, 1, nx, ny);
+
+            /* vertical bonds */
+            /* neighbours in the bulk */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx; i++){
+                for (j=1; j<ny-1; j++){
+                    n = cellnb(i, j, 1, nx, ny);
+                    cell[n].nneighb = 6;
+                    cell[n].nghb[0] = cellnb(i, j-1, 1, nx, ny);
+                    cell[n].nghb[1] = cellnb(i, j+1, 1, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[4] = cellnb(i, j+1, 0, nx, ny);
+                    cell[n].nghb[5] = cellnb(i-1, j+1, 0, nx, ny);
+                }
+            }
+            
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++){
+                n = cellnb(0, j, 1, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(0, j-1, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j+1, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(0, j+1, 0, nx, ny);
+            }
+            
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++){
+                n = cellnb(nx, j, 1, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(nx, j-1, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(nx, j+1, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(nx-1, j+1, 0, nx, ny);
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx; i++){
+                n = cellnb(i, 0, 1, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(i, 1, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, 0, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, 0, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i-1, 1, 0, nx, ny);
+                cell[n].nghb[4] = cellnb(i, 1, 0, nx, ny);
+            }
+            
+            /* top boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx; i++){
+                n = cellnb(i, ny-1, 1, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(i, ny-2, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, ny-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i, ny-1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i-1, ny, 0, nx, ny);
+                cell[n].nghb[4] = cellnb(i, ny, 0, nx, ny);
+            }
+            
+            /* corners */
+            n = cellnb(0, 0, 1, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(0, 1, 1, nx, ny);
+            cell[n].nghb[1] = cellnb(0, 0, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(0, 1, 0, nx, ny);
+            
+            n = cellnb(0, ny-1, 1, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(0, ny-2, 1, nx, ny);
+            cell[n].nghb[1] = cellnb(0, ny-1, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(0, ny, 0, nx, ny);
+
+            n = cellnb(nx, 0, 1, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx, 1, 1, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, 0, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(nx-1, 1, 0, nx, ny);
+                    
+            n = cellnb(nx, ny-1, 1, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx, ny-2, 1, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, ny, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(nx-1, ny-1, 0, nx, ny);
+            
+            return(nx + ny + 2*nx*ny);
+            break;
+        }
+        case (BC_HEX_SITE_DIRICHLET):
+        {
+            /* neighbours in the bulk  */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=1; j<ny-1; j+=2){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 6;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i+1, j+1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j+1, 0, nx, ny);
+                    cell[n].nghb[4] = cellnb(i+1, j-1,0, nx, ny);
+                    cell[n].nghb[5] = cellnb(i, j-1, 0, nx, ny);
+                }
+            }
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=2; j<ny-1; j+=2){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 6;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, j, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i-1, j+1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j+1, 0, nx, ny);
+                    cell[n].nghb[4] = cellnb(i-1, j-1,0, nx, ny);
+                    cell[n].nghb[5] = cellnb(i, j-1, 0, nx, ny);
+                }
+            }
+            
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j+=2){
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(1, j+1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(0, j-1, 0, nx, ny);
+                cell[n].nghb[4] = cellnb(1, j-1, 0, nx, ny);
+            }
+            #pragma omp parallel for private(j)
+            for (j=2; j<ny-1; j+=2){
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j-1, 0, nx, ny);
+            }
+            
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j+=2){
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx-2, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j-1, 0, nx, ny);
+            }
+            #pragma omp parallel for private(j)
+            for (j=2; j<ny-1; j+=2){
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 5;
+                cell[n].nghb[0] = cellnb(nx-2, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-2, j+1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j+1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(nx-2, j-1, 0, nx, ny);
+                cell[n].nghb[4] = cellnb(nx-1, j-1, 0, nx, ny);
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++){
+                n = cellnb(i, 0, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(i+1, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, 0, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i-1, 1, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i, 1, 0, nx, ny);
+            }
+            
+            /* top boundary */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++){
+                n = cellnb(i, ny-1, 0, nx, ny);
+                cell[n].nneighb = 4;
+                cell[n].nghb[0] = cellnb(i+1, ny-1, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, ny-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(i-1, ny-2, 0, nx, ny);
+                cell[n].nghb[3] = cellnb(i, ny-2, 0, nx, ny);
+            }
+            
+            /* corners */
+            n = cellnb(0, 0, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(1, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, 1, 0, nx, ny);
+            
+            n = cellnb(nx-1, 0, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx-2, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, 1, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(nx-2, 1, 0, nx, ny);
+            
+            n = cellnb(0, ny-1, 0, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(1, ny-1, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, ny-2, 0, nx, ny);
+            
+            n = cellnb(nx-1, ny-1, 0, nx, ny);
+            cell[n].nneighb = 3;
+            cell[n].nghb[0] = cellnb(nx-2, ny-1, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(nx-1, ny-2, 0, nx, ny);
+            cell[n].nghb[2] = cellnb(nx-2, ny-2, 0, nx, ny);
+            
+            return(nx*ny);
+            break;
+        }
+        case (BC_HEX_BOND_DIRICHLET):
+        {
+            /* neighbours in the bulk */
+            /* vertical bonds */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx; i++){
+                for (j=1; j<ny; j+=2){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i-1, j, 1, nx, ny);
+                    cell[n].nghb[1] = cellnb(i, j, 2, nx, ny);
+                    cell[n].nghb[2] = cellnb(i-1, j+1, 1, nx, ny);
+                    cell[n].nghb[3] = cellnb(i-1, j+1, 2, nx, ny);
+                }
+            }
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx; i++){
+                for (j=0; j<ny; j+=2){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i-1, j, 1, nx, ny);
+                    cell[n].nghb[1] = cellnb(i, j, 2, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j+1, 1, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j+1, 2, nx, ny);
+                }
+            }
+            /* NE-SW bonds */
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<nx-1; i++){
+                for (j=1; j<ny-1; j+=2){
+                    n = cellnb(i, j, 1, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i+1, j, 2, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j-1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j, 2, nx, ny);
+                }
+            }
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<nx-1; i++){
+                for (j=2; j<ny-1; j+=2){
+                    n = cellnb(i, j, 1, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i+1, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i+1, j, 2, nx, ny);
+                    cell[n].nghb[2] = cellnb(i+1, j-1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i, j, 2, nx, ny);
+                }
+            }
+            
+            /* NW-SE bonds */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<nx-1; i++){
+                for (j=1; j<ny-1; j+=2){
+                    n = cellnb(i, j, 2, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i, j, 1, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, j-1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i-1, j, 1, nx, ny);
+                }
+            }
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<nx; i++){
+                for (j=2; j<ny-1; j+=2){
+                    n = cellnb(i, j, 2, nx, ny);
+                    cell[n].nneighb = 4;
+                    cell[n].nghb[0] = cellnb(i, j, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i, j, 1, nx, ny);
+                    cell[n].nghb[2] = cellnb(i+1, j-1, 0, nx, ny);
+                    cell[n].nghb[3] = cellnb(i-1, j, 1, nx, ny);
+                }
+            }
+            
+            /* left boundary */
+            /* vertical bonds */
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j+=2)
+            {
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 2;
+                cell[n].nghb[0] = cellnb(0, j, 2, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j+1, 1, nx, ny);
+            }
+            /* NE-SW bonds */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j+=2)
+            {
+                n = cellnb(0, j, 1, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(0, j-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(1, j, 2, nx, ny);
+            }
+            /* NW-SE bonds */
+            #pragma omp parallel for private(j)
+            for (j=2; j<ny-1; j+=2)
+            {
+                n = cellnb(0, j, 2, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(0, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(1, j-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(0, j, 1, nx, ny);
+            }
+            
+            /* right boundary */
+            /* vertical bonds */
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j+=2)
+            {
+                n = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 2;
+                cell[n].nghb[0] = cellnb(nx-1, j+1, 2, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-2, j, 1, nx, ny);
+            }
+            /* NE-SW bonds */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j+=2)
+            {
+                n = cellnb(nx-1, j, 1, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j-1, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-1, j, 2, nx, ny);
+            }
+            #pragma omp parallel for private(j)
+            for (j=2; j<ny-1; j+=2)
+            {
+                n = cellnb(nx-1, j, 1, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx, j, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx, j, 2, nx, ny);
+                cell[n].nghb[2] = cellnb(nx, j-1, 0, nx, ny);
+            }
+            /* NW-SE bonds */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j+=2)
+            {
+                n = cellnb(nx-1, j, 2, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx-1, j-1, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-2, j, 1, nx, ny);
+            }
+            #pragma omp parallel for private(j)
+            for (j=2; j<ny-1; j+=2)
+            {
+                n = cellnb(nx-1, j, 2, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(nx-1, j-1, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(nx-1, j, 0, nx, ny);
+                cell[n].nghb[2] = cellnb(nx-2, j, 1, nx, ny);
+            }
+            
+            /* bottom boundary */
+            /* NE-SW bonds */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, 0, 1, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i+1, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i, 0, 2, nx, ny);
+                cell[n].nghb[2] = cellnb(i+1, 0, 2, nx, ny);
+            }
+            /* NW-SE bonds */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, 0, 2, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i, 0, 0, nx, ny);
+                cell[n].nghb[1] = cellnb(i-1, 0, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(i, 0, 1, nx, ny);
+            }
+            
+            /* top boundary */
+            /* vertical bonds */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, ny-1, 0, nx, ny);
+                cell[n].nneighb = 2;
+                cell[n].nghb[0] = cellnb(i-1, ny-1, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(i, ny-1, 2, nx, ny);
+            }
+            /* NE-SW bonds */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, ny-1, 1, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i+1, ny-1, 2, nx, ny);
+                cell[n].nghb[1] = cellnb(i, ny-1, 2, nx, ny);
+                cell[n].nghb[2] = cellnb(i+1, ny-2, 0, nx, ny);
+            }
+            /* NW-SE bonds */
+            #pragma omp parallel for private(i)
+            for (i=1; i<nx-1; i++)
+            {
+                n = cellnb(i, ny-1, 2, nx, ny);
+                cell[n].nneighb = 3;
+                cell[n].nghb[0] = cellnb(i-1, ny-1, 1, nx, ny);
+                cell[n].nghb[1] = cellnb(i, ny-1, 1, nx, ny);
+                cell[n].nghb[2] = cellnb(i+1, ny-2, 0, nx, ny);
+            }
+            
+            /* corners */
+            n = cellnb(0, 0, 2, nx, ny);
+            cell[n].nneighb = 2;
+            cell[n].nghb[0] = cellnb(0, 0, 0, nx, ny);
+            cell[n].nghb[1] = cellnb(0, 0, 1, nx, ny);
+            
+            return(3*(nx+1)*(ny+1));    /* TO BE CHECKED */
+            break;
+        }
+        case (BC_TRIANGLE_SITE_DIRICHLET):
+        {
+            /* neighbours in the bulk */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<2*nx-1; i++){
+                for (j=1; j<ny-1; j++){
+                    n = cellnb(i, j, 0, nx, ny);
+                    cell[n].nneighb = 3;
+                    cell[n].nghb[0] = cellnb(i+1,j,0,nx,ny);
+                    cell[n].nghb[1] = cellnb(i-1,j,0,nx,ny);
+                    if ((i+j)%2 == 0) cell[n].nghb[2] = cellnb(i,j+1,0,nx,ny);
+                    else cell[n].nghb[2] = cellnb(i,j-1,0,nx,ny);
+                }
+            }
+            
+            /* left boundary */
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j++){
+                n = cellnb(0, j, 0, nx, ny);
+                cell[n].nneighb = 2;
+                cell[n].nghb[0] = cellnb(1, j, 0, nx, ny);
+                if (j%2 == 0) cell[n].nghb[1] = cellnb(0, j+1, 0, nx, ny);
+                else cell[n].nghb[1] = cellnb(0, j-1, 0, nx, ny);
+            }
+            
+            /* right boundary */
+            #pragma omp parallel for private(j)
+            for (j=1; j<ny-1; j++){
+                n = cellnb(2*nx-1, j, 0, nx, ny);
+                cell[n].nneighb = 2;
+                cell[n].nghb[0] = cellnb(2*nx-2, j, 0, nx, ny);
+                if (j%2 == 1) cell[n].nghb[1] = cellnb(2*nx-1, j+1, 0, nx, ny);
+                else cell[n].nghb[1] = cellnb(2*nx-1, j-1, 0, nx, ny);
+            }
+            
+            /* bottom boundary */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<2*nx-1; i++){
+                n = cellnb(i, 0, 0, nx, ny);
+                if (i%2 == 0)
+                {
+                    cell[n].nneighb = 3;
+                    cell[n].nghb[0] = cellnb(i-1, 0, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i+1, 0, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, 1, 0, nx, ny);
+                }
+                else 
+                {
+                    cell[n].nneighb = 2;
+                    cell[n].nghb[0] = cellnb(i+1, 0, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i-1, 0, 0, nx, ny);
+                }
+            }
+            
+            /* top boundary */
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<2*nx-1; i++){
+                n = cellnb(i, ny-1, 0, nx, ny);
+                if (i%2 == 1)
+                {
+                    cell[n].nneighb = 3;
+                    cell[n].nghb[0] = cellnb(i-1, ny-1, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i+1, ny-1, 0, nx, ny);
+                    cell[n].nghb[2] = cellnb(i, ny-2, 0, nx, ny);
+                }
+                else 
+                {
+                    cell[n].nneighb = 2;
+                    cell[n].nghb[0] = cellnb(i-1, ny-1, 0, nx, ny);
+                    cell[n].nghb[1] = cellnb(i+1, ny-1, 0, nx, ny);
+                }
+            }
+            
+            
+            /* corners (at the right) */
+            n = cellnb(2*nx-1, 0, 0, nx, ny);
+            cell[n].nneighb = 1;
+            cell[n].nghb[0] = cellnb(2*nx-2, 0, 0, nx, ny);
+            
+            n = cellnb(2*nx-1, ny-1, 0, nx, ny);
+            cell[n].nneighb = 1;
+            cell[n].nghb[0] = cellnb(2*nx-2, ny-1, 0, nx, ny);
+            
+            
+            return(2*nx*ny);
+            break;
+        }
+    }
+    printf("Done\n");
+}
+
+
+void init_cell_probabilities(t_perco *cell, int ncells)
+/* initialize the probabilities of cells being open */
+{
+    int i;
+    
+    printf("Initializing cell probabilities ...");
+    #pragma omp parallel for private(i)
+    for (i=0; i<ncells; i++)
+    {
+        cell[i].proba = (double)rand()/RAND_MAX;
+    }
+    printf("Done\n");
+}
+
+void init_cell_state(t_perco *cell, double p, int ncells, int first)
+/* initialize the probabilities of cells being open */
+{
+    int i;
+    
+    printf("Initializing cell state ...");
+    #pragma omp parallel for private(i)
+    for (i=0; i<ncells; i++)
+    {
+        if (cell[i].proba < p) 
+        {
+            cell[i].open = 1;
+            cell[i].active = 1;
+        }
+        else 
+        {
+            cell[i].open = 0;
+            cell[i].active = 0;
+        }
+        cell[i].flooded = 0;
+        if (first) cell[i].previous_cluster = 0;
+        else cell[i].previous_cluster = cell[i].cluster; 
+        cell[i].cluster = 0;
+        cell[i].tested = 0;
+    }
+    printf("Done\n");
+}
+
+int init_flooded_cells(t_perco *cell, int nx, int ny, t_perco* *pstack)
+/* make the left row of cells flooded, returns number of flooded cells */
+{
+    int j, ishift;
+    
+//     printf("Initializing flooded cells ...");
+    switch (graphical_rep(LATTICE)) {
+        case (PLOT_SQUARES): 
+        {
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j++)
+            {
+                cell[j].open = 1;
+                cell[j].flooded = 1;
+                cell[j].cluster = 1;
+                cell[j].previous_cluster = 1;
+                cell[j].active = 1;
+                pstack[j] = &cell[j];
+            }
+            return(ny);
+        }
+        case (PLOT_SQUARE_BONDS): 
+        {
+            ishift = nx*(ny+1);
+            
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j++)
+            {
+                cell[ishift + j].open = 1;
+                cell[ishift + j].flooded = 1;
+                cell[ishift + j].cluster = 1;
+                cell[ishift + j].previous_cluster = 1;
+                cell[ishift + j].active = 1;
+                pstack[j] = &cell[ishift + j];
+            }
+            return(ny);
+        }
+        case (PLOT_HEX): 
+        {
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j++)
+            {
+                cell[j*nx].open = 1;
+                cell[j*nx].flooded = 1;
+                cell[j*nx].cluster = 1;
+                cell[j*nx].previous_cluster = 1;
+                cell[j*nx].active = 1;
+                pstack[j] = &cell[j*nx];
+            }
+            return(ny);
+        }
+        case (PLOT_HEX_BONDS): 
+        {
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j+=2)
+            {
+                cell[j*nx].open = 1;
+                cell[j*nx].flooded = 1;
+                cell[j*nx].cluster = 1;
+                cell[j*nx].previous_cluster = 1;
+                cell[j*nx].active = 1;
+                pstack[j] = &cell[j*nx];
+            }
+            return(ny);
+        }
+        case (PLOT_TRIANGLE): 
+        {
+            #pragma omp parallel for private(j)
+            for (j=0; j<ny; j++)
+            {
+                cell[2*j*nx].open = 1;
+                cell[2*j*nx].flooded = 1;
+                cell[2*j*nx].cluster = 1;
+                cell[2*j*nx].previous_cluster = 1;
+                cell[2*j*nx].active = 1;
+                pstack[j] = &cell[2*j*nx];
+            }
+            return(ny);
+        }
+    }
+//     printf("Done\n");
+}
+
+
+int count_open_cells(t_perco *cell, int ncells)
+/* count the number of active cells */
+{
+    int n = 0, i;
+    
+    for (i=0; i<ncells; i++) if (cell[i].open) n++;
+    return(n);
+}
+
+int count_flooded_cells(t_perco *cell, int ncells)
+/* count the number of active cells */
+{
+    int n = 0, i;
+    
+    for (i=0; i<ncells; i++) if ((cell[i].open)&&(cell[i].flooded)) n++;
+    return(n);
+}
+
+int count_active_cells(t_perco *cell, int ncells)
+/* count the number of active cells */
+{
+    int n = 0, i;
+    
+    for (i=0; i<ncells; i++) 
+        if ((cell[i].active)&&(cell[i].flooded)) n++;
+        
+    return(n);
+}
+
+
+int find_open_cells(t_perco *cell, int ncells, int position, t_perco* *pstack, int *stacksize, int cluster)
+/* look for open neighbours of start_cell, returns difference in open cells */
+{
+    int k, nopen = 0, n;
+    
+//     printf("Position %i, open %i\n", position, pstack[position]->open);
+    
+    if (!pstack[position]->open) return(-1);
+    
+    pstack[position]->cluster = cluster;
+    
+    for (k=0; k<pstack[position]->nneighb; k++)
+    {
+        n = pstack[position]->nghb[k];
+        if ((cell[n].open)&&(cell[n].cluster != cluster))
+        {
+            cell[n].flooded = 1;
+            cell[n].tested = 1;
+            cell[n].cluster = cluster;
+            cell[n].active = 1;
+            nopen++;
+            if (*stacksize < ncells)
+            {
+                (*stacksize)++;
+                pstack[*stacksize-1] = &cell[n];
+            }
+        }
+    }
+
+    if (nopen == 0) 
+    {
+        pstack[position]->active = 0;
+        return(-1);
+    }
+    else return(nopen);
+}
+
+
+int find_percolation_cluster(t_perco *cell, int ncells, t_perco* *pstack, int nstack)
+/* new version of cluster search algorithm, using stacks; returns number of flooded cells */
+{
+    int position = 0, i, stacksize = nstack, nactive = nstack; 
+    
+    while (nactive > 0) 
+    {
+        /* find an active cell in stack */
+        while (!pstack[position]->active) position++;
+        if (position == stacksize) position = 0;
+        
+        nactive += find_open_cells(cell, ncells, position, pstack, &stacksize, 1);
+    }
+    printf("Stack size %i\n", stacksize);
+    return(stacksize);
+}
+
+
+int find_all_clusters(t_perco *cell, int ncells, int reset, int *max_cluster_label)
+/* find all clusters of open cells; returns number of clusters */
+{
+    int nclusters = 0, i, j, nactive, stacksize, position, newcluster = 1, maxlabel = 0;
+    t_perco **cstack;
+    static int cluster;
+    
+    cstack = (t_perco* *)malloc(ncells*sizeof(struct t_perco *));
+    
+    if (reset) cluster = CLUSTER_SHIFT;
+    else cluster = (*max_cluster_label) + CLUSTER_SHIFT;   /* give higher labels than before to new clusters */
+
+    for (i=0; i<ncells; i++) if ((cell[i].open)&&(cell[i].cluster != 1)&&(!cell[i].tested))
+    {
+        position = 0;
+        stacksize = 1;
+        nactive = 1;
+        cstack[0] = &cell[i];
+        
+        /* choice of new color */ 
+        if (cell[i].previous_cluster >= CLUSTER_SHIFT) newcluster = cell[i].previous_cluster;
+        else 
+        {
+            newcluster = cluster;
+            cluster++;
+        }
+        cell[i].cluster = newcluster;
+        cell[i].active = 1;
+        cell[i].tested = 1;
+        
+        while (nactive > 0) 
+        {
+            /* find an active cell in stack */
+            while (!cstack[position]->active) position++;
+            if (position == stacksize) position = 0;
+        
+            nactive += find_open_cells(cell, ncells, position, cstack, &stacksize, newcluster);
+        }
+//         printf("Found cluster %i with %i active cells\n", cluster, stacksize);
+
+        nclusters++;
+    }
+    
+    free(cstack);
+    
+    /* determine maximal cluster label */
+//     #pragma omp parallel for private(i)
+    for (i=0; i<ncells; i++) if ((cell[i].open)&&(cell[i].cluster > maxlabel)) maxlabel = cell[i].cluster;
+    
+    printf("Cluster = %i, maxlabel = %i\n", cluster, maxlabel);
+    
+    *max_cluster_label = maxlabel;
+    return(nclusters);
+}
+
+void print_cluster_sizes(t_perco *cell, int ncells, int *cluster_sizes)
+/* for debugging purposes */
+{
+    int j;
+    
+    for (j=0; j<ncells; j++) if (cell[j].open == 1) 
+            printf("(cell %i: cluster %i: size %i)\n", j, cell[j].cluster, cluster_sizes[cell[j].cluster]);
+    sleep(1);
+    printf("\n\n");
+}
+
+int find_cluster_sizes(t_perco *cell, int ncells, int *cluster_sizes, int *maxclusterlabel)
+/* determine the sizes of all clusters; returns the size of largest cluster */
+{
+    int i, max = 0, maxlabel = 0;
+    
+//     if (maxclusterlabel > ncells) maxclusterlabel = ncells;
+    
+    /* determine label of larges cluster */
+    for (i=0; i<ncells; i++) if (cell[i].cluster > maxlabel) maxlabel = cell[i].cluster;
+    *maxclusterlabel = maxlabel + 1;
+
+    #pragma omp parallel for private(i)
+    for (i=0; i<*maxclusterlabel; i++) cluster_sizes[i] = 0;
+    
+    for (i=0; i<ncells; i++) if (cell[i].open) cluster_sizes[cell[i].cluster]++;
+    
+    if (DEBUG)
+    {
+        printf("Computed %i cluster sizes\n", *maxclusterlabel);
+        print_cluster_sizes(cell, ncells, cluster_sizes);
+    }
+    
+    for (i=0; i<*maxclusterlabel; i++) if (cluster_sizes[i] > max) max = cluster_sizes[i];
+    
+//     for (i=0; i<ncells; i++) if (cell[i].open) printf("%i ", cluster_sizes[cell[i].cluster]);
+    
+    printf("Max cluster label %i\n", *maxclusterlabel);
+    printf("Largest cluster has %i cells\n", max); 
+    return(max);
+}