diff --git a/.gitignore b/.gitignore
index 564f67366f2ec7f7d26db2df33a0ff074479cc50..cc2ab72bf338c8744a5bba0de6a9695e6b6f4a31 100644
--- a/.gitignore
+++ b/.gitignore
@@ -85,3 +85,4 @@ scripts/stefan copy/
 scripts/stefan/
 test.py
 /Release/
+build*
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9bc67a547a7971b7e725bf4d2b5b92e97983128e..22b48d8d43f24a96c1abf5ff0d0daa089e54cba6 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -25,6 +25,7 @@ configure-linux:
   script:
     - mkdir -p build
     - cd build
+    - which git
     - cmake -DCMAKE_BUILD_TYPE=Debug -DBoost_NO_SYSTEM_PATHS=true -DBOOST_ROOT=~/boost_1_61_0  ..
     - echo "configure | ${CI_PROJECT_DIR}"
   stage: configure
@@ -48,7 +49,8 @@ make-linux:
     - cd build
     - make -j$nproc
     - echo "compile | ${CI_PROJECT_DIR}"
-
+  tags:
+    - linux
 
 after_script:
   - echo "End CI" # todo: run report script
diff --git a/.vs/dirent.h b/.vs/dirent.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcc2f556a42be6054e60fb82a8a3100f336db9b7
--- /dev/null
+++ b/.vs/dirent.h
@@ -0,0 +1,1157 @@
+/*
+ * Dirent interface for Microsoft Visual Studio
+ *
+ * Copyright (C) 2006-2012 Toni Ronkko
+ * This file is part of dirent.  Dirent may be freely distributed
+ * under the MIT license.  For all details and documentation, see
+ * https://github.com/tronkko/dirent
+ */
+#ifndef DIRENT_H
+#define DIRENT_H
+
+/*
+ * Include windows.h without Windows Sockets 1.1 to prevent conflicts with
+ * Windows Sockets 2.0.
+ */
+#ifndef WIN32_LEAN_AND_MEAN
+#   define WIN32_LEAN_AND_MEAN
+#endif
+#include <windows.h>
+
+#include <stdio.h>
+#include <stdarg.h>
+#include <wchar.h>
+#include <string.h>
+#include <stdlib.h>
+#include <malloc.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <errno.h>
+
+/* Indicates that d_type field is available in dirent structure */
+#define _DIRENT_HAVE_D_TYPE
+
+/* Indicates that d_namlen field is available in dirent structure */
+#define _DIRENT_HAVE_D_NAMLEN
+
+/* Entries missing from MSVC 6.0 */
+#if !defined(FILE_ATTRIBUTE_DEVICE)
+#   define FILE_ATTRIBUTE_DEVICE 0x40
+#endif
+
+/* File type and permission flags for stat(), general mask */
+#if !defined(S_IFMT)
+#   define S_IFMT _S_IFMT
+#endif
+
+/* Directory bit */
+#if !defined(S_IFDIR)
+#   define S_IFDIR _S_IFDIR
+#endif
+
+/* Character device bit */
+#if !defined(S_IFCHR)
+#   define S_IFCHR _S_IFCHR
+#endif
+
+/* Pipe bit */
+#if !defined(S_IFFIFO)
+#   define S_IFFIFO _S_IFFIFO
+#endif
+
+/* Regular file bit */
+#if !defined(S_IFREG)
+#   define S_IFREG _S_IFREG
+#endif
+
+/* Read permission */
+#if !defined(S_IREAD)
+#   define S_IREAD _S_IREAD
+#endif
+
+/* Write permission */
+#if !defined(S_IWRITE)
+#   define S_IWRITE _S_IWRITE
+#endif
+
+/* Execute permission */
+#if !defined(S_IEXEC)
+#   define S_IEXEC _S_IEXEC
+#endif
+
+/* Pipe */
+#if !defined(S_IFIFO)
+#   define S_IFIFO _S_IFIFO
+#endif
+
+/* Block device */
+#if !defined(S_IFBLK)
+#   define S_IFBLK 0
+#endif
+
+/* Link */
+#if !defined(S_IFLNK)
+#   define S_IFLNK 0
+#endif
+
+/* Socket */
+#if !defined(S_IFSOCK)
+#   define S_IFSOCK 0
+#endif
+
+/* Read user permission */
+#if !defined(S_IRUSR)
+#   define S_IRUSR S_IREAD
+#endif
+
+/* Write user permission */
+#if !defined(S_IWUSR)
+#   define S_IWUSR S_IWRITE
+#endif
+
+/* Execute user permission */
+#if !defined(S_IXUSR)
+#   define S_IXUSR 0
+#endif
+
+/* Read group permission */
+#if !defined(S_IRGRP)
+#   define S_IRGRP 0
+#endif
+
+/* Write group permission */
+#if !defined(S_IWGRP)
+#   define S_IWGRP 0
+#endif
+
+/* Execute group permission */
+#if !defined(S_IXGRP)
+#   define S_IXGRP 0
+#endif
+
+/* Read others permission */
+#if !defined(S_IROTH)
+#   define S_IROTH 0
+#endif
+
+/* Write others permission */
+#if !defined(S_IWOTH)
+#   define S_IWOTH 0
+#endif
+
+/* Execute others permission */
+#if !defined(S_IXOTH)
+#   define S_IXOTH 0
+#endif
+
+/* Maximum length of file name */
+#if !defined(PATH_MAX)
+#   define PATH_MAX MAX_PATH
+#endif
+#if !defined(FILENAME_MAX)
+#   define FILENAME_MAX MAX_PATH
+#endif
+#if !defined(NAME_MAX)
+#   define NAME_MAX FILENAME_MAX
+#endif
+
+/* File type flags for d_type */
+#define DT_UNKNOWN 0
+#define DT_REG S_IFREG
+#define DT_DIR S_IFDIR
+#define DT_FIFO S_IFIFO
+#define DT_SOCK S_IFSOCK
+#define DT_CHR S_IFCHR
+#define DT_BLK S_IFBLK
+#define DT_LNK S_IFLNK
+
+/* Macros for converting between st_mode and d_type */
+#define IFTODT(mode) ((mode) & S_IFMT)
+#define DTTOIF(type) (type)
+
+/*
+ * File type macros.  Note that block devices, sockets and links cannot be
+ * distinguished on Windows and the macros S_ISBLK, S_ISSOCK and S_ISLNK are
+ * only defined for compatibility.  These macros should always return false
+ * on Windows.
+ */
+#if !defined(S_ISFIFO)
+#   define S_ISFIFO(mode) (((mode) & S_IFMT) == S_IFIFO)
+#endif
+#if !defined(S_ISDIR)
+#   define S_ISDIR(mode) (((mode) & S_IFMT) == S_IFDIR)
+#endif
+#if !defined(S_ISREG)
+#   define S_ISREG(mode) (((mode) & S_IFMT) == S_IFREG)
+#endif
+#if !defined(S_ISLNK)
+#   define S_ISLNK(mode) (((mode) & S_IFMT) == S_IFLNK)
+#endif
+#if !defined(S_ISSOCK)
+#   define S_ISSOCK(mode) (((mode) & S_IFMT) == S_IFSOCK)
+#endif
+#if !defined(S_ISCHR)
+#   define S_ISCHR(mode) (((mode) & S_IFMT) == S_IFCHR)
+#endif
+#if !defined(S_ISBLK)
+#   define S_ISBLK(mode) (((mode) & S_IFMT) == S_IFBLK)
+#endif
+
+/* Return the exact length of the file name without zero terminator */
+#define _D_EXACT_NAMLEN(p) ((p)->d_namlen)
+
+/* Return the maximum size of a file name */
+#define _D_ALLOC_NAMLEN(p) ((PATH_MAX)+1)
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+/* Wide-character version */
+struct _wdirent {
+    /* Always zero */
+    long d_ino;
+
+    /* File position within stream */
+    long d_off;
+
+    /* Structure size */
+    unsigned short d_reclen;
+
+    /* Length of name without \0 */
+    size_t d_namlen;
+
+    /* File type */
+    int d_type;
+
+    /* File name */
+    wchar_t d_name[PATH_MAX+1];
+};
+typedef struct _wdirent _wdirent;
+
+struct _WDIR {
+    /* Current directory entry */
+    struct _wdirent ent;
+
+    /* Private file data */
+    WIN32_FIND_DATAW data;
+
+    /* True if data is valid */
+    int cached;
+
+    /* Win32 search handle */
+    HANDLE handle;
+
+    /* Initial directory name */
+    wchar_t *patt;
+};
+typedef struct _WDIR _WDIR;
+
+/* Multi-byte character version */
+struct dirent {
+    /* Always zero */
+    long d_ino;
+
+    /* File position within stream */
+    long d_off;
+
+    /* Structure size */
+    unsigned short d_reclen;
+
+    /* Length of name without \0 */
+    size_t d_namlen;
+
+    /* File type */
+    int d_type;
+
+    /* File name */
+    char d_name[PATH_MAX+1];
+};
+typedef struct dirent dirent;
+
+struct DIR {
+    struct dirent ent;
+    struct _WDIR *wdirp;
+};
+typedef struct DIR DIR;
+
+
+/* Dirent functions */
+static DIR *opendir (const char *dirname);
+static _WDIR *_wopendir (const wchar_t *dirname);
+
+static struct dirent *readdir (DIR *dirp);
+static struct _wdirent *_wreaddir (_WDIR *dirp);
+
+static int readdir_r(
+    DIR *dirp, struct dirent *entry, struct dirent **result);
+static int _wreaddir_r(
+    _WDIR *dirp, struct _wdirent *entry, struct _wdirent **result);
+
+static int closedir (DIR *dirp);
+static int _wclosedir (_WDIR *dirp);
+
+static void rewinddir (DIR* dirp);
+static void _wrewinddir (_WDIR* dirp);
+
+static int scandir (const char *dirname, struct dirent ***namelist,
+    int (*filter)(const struct dirent*),
+    int (*compare)(const void *, const void *));
+
+static int alphasort (const struct dirent **a, const struct dirent **b);
+
+static int versionsort (const struct dirent **a, const struct dirent **b);
+
+
+/* For compatibility with Symbian */
+#define wdirent _wdirent
+#define WDIR _WDIR
+#define wopendir _wopendir
+#define wreaddir _wreaddir
+#define wclosedir _wclosedir
+#define wrewinddir _wrewinddir
+
+
+/* Internal utility functions */
+static WIN32_FIND_DATAW *dirent_first (_WDIR *dirp);
+static WIN32_FIND_DATAW *dirent_next (_WDIR *dirp);
+
+static int dirent_mbstowcs_s(
+    size_t *pReturnValue,
+    wchar_t *wcstr,
+    size_t sizeInWords,
+    const char *mbstr,
+    size_t count);
+
+static int dirent_wcstombs_s(
+    size_t *pReturnValue,
+    char *mbstr,
+    size_t sizeInBytes,
+    const wchar_t *wcstr,
+    size_t count);
+
+static void dirent_set_errno (int error);
+
+
+/*
+ * Open directory stream DIRNAME for read and return a pointer to the
+ * internal working area that is used to retrieve individual directory
+ * entries.
+ */
+static _WDIR*
+_wopendir(
+    const wchar_t *dirname)
+{
+    _WDIR *dirp = NULL;
+    int error;
+
+    /* Must have directory name */
+    if (dirname == NULL  ||  dirname[0] == '\0') {
+        dirent_set_errno (ENOENT);
+        return NULL;
+    }
+
+    /* Allocate new _WDIR structure */
+    dirp = (_WDIR*) malloc (sizeof (struct _WDIR));
+    if (dirp != NULL) {
+        DWORD n;
+
+        /* Reset _WDIR structure */
+        dirp->handle = INVALID_HANDLE_VALUE;
+        dirp->patt = NULL;
+        dirp->cached = 0;
+
+        /* Compute the length of full path plus zero terminator
+         *
+         * Note that on WinRT there's no way to convert relative paths
+         * into absolute paths, so just assume it is an absolute path.
+         */
+#       if defined(WINAPI_FAMILY) && (WINAPI_FAMILY == WINAPI_FAMILY_PHONE_APP)
+            n = wcslen(dirname);
+#       else
+            n = GetFullPathNameW (dirname, 0, NULL, NULL);
+#       endif
+
+        /* Allocate room for absolute directory name and search pattern */
+        dirp->patt = (wchar_t*) malloc (sizeof (wchar_t) * n + 16);
+        if (dirp->patt) {
+
+            /*
+             * Convert relative directory name to an absolute one.  This
+             * allows rewinddir() to function correctly even when current
+             * working directory is changed between opendir() and rewinddir().
+             *
+             * Note that on WinRT there's no way to convert relative paths
+             * into absolute paths, so just assume it is an absolute path.
+             */
+#           if defined(WINAPI_FAMILY) && (WINAPI_FAMILY == WINAPI_FAMILY_PHONE_APP)
+                wcsncpy_s(dirp->patt, n+1, dirname, n);
+#           else
+                n = GetFullPathNameW (dirname, n, dirp->patt, NULL);
+#           endif
+            if (n > 0) {
+                wchar_t *p;
+
+                /* Append search pattern \* to the directory name */
+                p = dirp->patt + n;
+                if (dirp->patt < p) {
+                    switch (p[-1]) {
+                    case '\\':
+                    case '/':
+                    case ':':
+                        /* Directory ends in path separator, e.g. c:\temp\ */
+                        /*NOP*/;
+                        break;
+
+                    default:
+                        /* Directory name doesn't end in path separator */
+                        *p++ = '\\';
+                    }
+                }
+                *p++ = '*';
+                *p = '\0';
+
+                /* Open directory stream and retrieve the first entry */
+                if (dirent_first (dirp)) {
+                    /* Directory stream opened successfully */
+                    error = 0;
+                } else {
+                    /* Cannot retrieve first entry */
+                    error = 1;
+                    dirent_set_errno (ENOENT);
+                }
+
+            } else {
+                /* Cannot retrieve full path name */
+                dirent_set_errno (ENOENT);
+                error = 1;
+            }
+
+        } else {
+            /* Cannot allocate memory for search pattern */
+            error = 1;
+        }
+
+    } else {
+        /* Cannot allocate _WDIR structure */
+        error = 1;
+    }
+
+    /* Clean up in case of error */
+    if (error  &&  dirp) {
+        _wclosedir (dirp);
+        dirp = NULL;
+    }
+
+    return dirp;
+}
+
+/*
+ * Read next directory entry.
+ *
+ * Returns pointer to static directory entry which may be overwritten by
+ * subsequent calls to _wreaddir().
+ */
+static struct _wdirent*
+_wreaddir(
+    _WDIR *dirp)
+{
+    struct _wdirent *entry;
+
+    /*
+     * Read directory entry to buffer.  We can safely ignore the return value
+     * as entry will be set to NULL in case of error.
+     */
+    (void) _wreaddir_r (dirp, &dirp->ent, &entry);
+
+    /* Return pointer to statically allocated directory entry */
+    return entry;
+}
+
+/*
+ * Read next directory entry.
+ *
+ * Returns zero on success.  If end of directory stream is reached, then sets
+ * result to NULL and returns zero.
+ */
+static int
+_wreaddir_r(
+    _WDIR *dirp,
+    struct _wdirent *entry,
+    struct _wdirent **result)
+{
+    WIN32_FIND_DATAW *datap;
+
+    /* Read next directory entry */
+    datap = dirent_next (dirp);
+    if (datap) {
+        size_t n;
+        DWORD attr;
+
+        /*
+         * Copy file name as wide-character string.  If the file name is too
+         * long to fit in to the destination buffer, then truncate file name
+         * to PATH_MAX characters and zero-terminate the buffer.
+         */
+        n = 0;
+        while (n < PATH_MAX  &&  datap->cFileName[n] != 0) {
+            entry->d_name[n] = datap->cFileName[n];
+            n++;
+        }
+        entry->d_name[n] = 0;
+
+        /* Length of file name excluding zero terminator */
+        entry->d_namlen = n;
+
+        /* File type */
+        attr = datap->dwFileAttributes;
+        if ((attr & FILE_ATTRIBUTE_DEVICE) != 0) {
+            entry->d_type = DT_CHR;
+        } else if ((attr & FILE_ATTRIBUTE_DIRECTORY) != 0) {
+            entry->d_type = DT_DIR;
+        } else {
+            entry->d_type = DT_REG;
+        }
+
+        /* Reset dummy fields */
+        entry->d_ino = 0;
+        entry->d_off = 0;
+        entry->d_reclen = sizeof (struct _wdirent);
+
+        /* Set result address */
+        *result = entry;
+
+    } else {
+
+        /* Return NULL to indicate end of directory */
+        *result = NULL;
+
+    }
+
+    return /*OK*/0;
+}
+
+/*
+ * Close directory stream opened by opendir() function.  This invalidates the
+ * DIR structure as well as any directory entry read previously by
+ * _wreaddir().
+ */
+static int
+_wclosedir(
+    _WDIR *dirp)
+{
+    int ok;
+    if (dirp) {
+
+        /* Release search handle */
+        if (dirp->handle != INVALID_HANDLE_VALUE) {
+            FindClose (dirp->handle);
+            dirp->handle = INVALID_HANDLE_VALUE;
+        }
+
+        /* Release search pattern */
+        if (dirp->patt) {
+            free (dirp->patt);
+            dirp->patt = NULL;
+        }
+
+        /* Release directory structure */
+        free (dirp);
+        ok = /*success*/0;
+
+    } else {
+
+        /* Invalid directory stream */
+        dirent_set_errno (EBADF);
+        ok = /*failure*/-1;
+
+    }
+    return ok;
+}
+
+/*
+ * Rewind directory stream such that _wreaddir() returns the very first
+ * file name again.
+ */
+static void
+_wrewinddir(
+    _WDIR* dirp)
+{
+    if (dirp) {
+        /* Release existing search handle */
+        if (dirp->handle != INVALID_HANDLE_VALUE) {
+            FindClose (dirp->handle);
+        }
+
+        /* Open new search handle */
+        dirent_first (dirp);
+    }
+}
+
+/* Get first directory entry (internal) */
+static WIN32_FIND_DATAW*
+dirent_first(
+    _WDIR *dirp)
+{
+    WIN32_FIND_DATAW *datap;
+
+    /* Open directory and retrieve the first entry */
+    dirp->handle = FindFirstFileExW(
+        dirp->patt, FindExInfoStandard, &dirp->data,
+        FindExSearchNameMatch, NULL, 0);
+    if (dirp->handle != INVALID_HANDLE_VALUE) {
+
+        /* a directory entry is now waiting in memory */
+        datap = &dirp->data;
+        dirp->cached = 1;
+
+    } else {
+
+        /* Failed to re-open directory: no directory entry in memory */
+        dirp->cached = 0;
+        datap = NULL;
+
+    }
+    return datap;
+}
+
+/*
+ * Get next directory entry (internal).
+ *
+ * Returns 
+ */
+static WIN32_FIND_DATAW*
+dirent_next(
+    _WDIR *dirp)
+{
+    WIN32_FIND_DATAW *p;
+
+    /* Get next directory entry */
+    if (dirp->cached != 0) {
+
+        /* A valid directory entry already in memory */
+        p = &dirp->data;
+        dirp->cached = 0;
+
+    } else if (dirp->handle != INVALID_HANDLE_VALUE) {
+
+        /* Get the next directory entry from stream */
+        if (FindNextFileW (dirp->handle, &dirp->data) != FALSE) {
+            /* Got a file */
+            p = &dirp->data;
+        } else {
+            /* The very last entry has been processed or an error occurred */
+            FindClose (dirp->handle);
+            dirp->handle = INVALID_HANDLE_VALUE;
+            p = NULL;
+        }
+
+    } else {
+
+        /* End of directory stream reached */
+        p = NULL;
+
+    }
+
+    return p;
+}
+
+/*
+ * Open directory stream using plain old C-string.
+ */
+static DIR*
+opendir(
+    const char *dirname) 
+{
+    struct DIR *dirp;
+    int error;
+
+    /* Must have directory name */
+    if (dirname == NULL  ||  dirname[0] == '\0') {
+        dirent_set_errno (ENOENT);
+        return NULL;
+    }
+
+    /* Allocate memory for DIR structure */
+    dirp = (DIR*) malloc (sizeof (struct DIR));
+    if (dirp) {
+        wchar_t wname[PATH_MAX + 1];
+        size_t n;
+
+        /* Convert directory name to wide-character string */
+        error = dirent_mbstowcs_s(
+            &n, wname, PATH_MAX + 1, dirname, PATH_MAX + 1);
+        if (!error) {
+
+            /* Open directory stream using wide-character name */
+            dirp->wdirp = _wopendir (wname);
+            if (dirp->wdirp) {
+                /* Directory stream opened */
+                error = 0;
+            } else {
+                /* Failed to open directory stream */
+                error = 1;
+            }
+
+        } else {
+            /*
+             * Cannot convert file name to wide-character string.  This
+             * occurs if the string contains invalid multi-byte sequences or
+             * the output buffer is too small to contain the resulting
+             * string.
+             */
+            error = 1;
+        }
+
+    } else {
+        /* Cannot allocate DIR structure */
+        error = 1;
+    }
+
+    /* Clean up in case of error */
+    if (error  &&  dirp) {
+        free (dirp);
+        dirp = NULL;
+    }
+
+    return dirp;
+}
+
+/*
+ * Read next directory entry.
+ */
+static struct dirent*
+readdir(
+    DIR *dirp)
+{
+    struct dirent *entry;
+
+    /*
+     * Read directory entry to buffer.  We can safely ignore the return value
+     * as entry will be set to NULL in case of error.
+     */
+    (void) readdir_r (dirp, &dirp->ent, &entry);
+
+    /* Return pointer to statically allocated directory entry */
+    return entry;
+}
+
+/*
+ * Read next directory entry into called-allocated buffer.
+ *
+ * Returns zero on success.  If the end of directory stream is reached, then
+ * sets result to NULL and returns zero.
+ */
+static int
+readdir_r(
+    DIR *dirp,
+    struct dirent *entry,
+    struct dirent **result)
+{
+    WIN32_FIND_DATAW *datap;
+
+    /* Read next directory entry */
+    datap = dirent_next (dirp->wdirp);
+    if (datap) {
+        size_t n;
+        int error;
+
+        /* Attempt to convert file name to multi-byte string */
+        error = dirent_wcstombs_s(
+            &n, entry->d_name, PATH_MAX + 1, datap->cFileName, PATH_MAX + 1);
+
+        /*
+         * If the file name cannot be represented by a multi-byte string,
+         * then attempt to use old 8+3 file name.  This allows traditional
+         * Unix-code to access some file names despite of unicode
+         * characters, although file names may seem unfamiliar to the user.
+         *
+         * Be ware that the code below cannot come up with a short file
+         * name unless the file system provides one.  At least
+         * VirtualBox shared folders fail to do this.
+         */
+        if (error  &&  datap->cAlternateFileName[0] != '\0') {
+            error = dirent_wcstombs_s(
+                &n, entry->d_name, PATH_MAX + 1,
+                datap->cAlternateFileName, PATH_MAX + 1);
+        }
+
+        if (!error) {
+            DWORD attr;
+
+            /* Length of file name excluding zero terminator */
+            entry->d_namlen = n - 1;
+
+            /* File attributes */
+            attr = datap->dwFileAttributes;
+            if ((attr & FILE_ATTRIBUTE_DEVICE) != 0) {
+                entry->d_type = DT_CHR;
+            } else if ((attr & FILE_ATTRIBUTE_DIRECTORY) != 0) {
+                entry->d_type = DT_DIR;
+            } else {
+                entry->d_type = DT_REG;
+            }
+
+            /* Reset dummy fields */
+            entry->d_ino = 0;
+            entry->d_off = 0;
+            entry->d_reclen = sizeof (struct dirent);
+
+        } else {
+
+            /*
+             * Cannot convert file name to multi-byte string so construct
+             * an erroneous directory entry and return that.  Note that
+             * we cannot return NULL as that would stop the processing
+             * of directory entries completely.
+             */
+            entry->d_name[0] = '?';
+            entry->d_name[1] = '\0';
+            entry->d_namlen = 1;
+            entry->d_type = DT_UNKNOWN;
+            entry->d_ino = 0;
+            entry->d_off = -1;
+            entry->d_reclen = 0;
+
+        }
+
+        /* Return pointer to directory entry */
+        *result = entry;
+
+    } else {
+
+        /* No more directory entries */
+        *result = NULL;
+
+    }
+
+    return /*OK*/0;
+}
+
+/*
+ * Close directory stream.
+ */
+static int
+closedir(
+    DIR *dirp)
+{
+    int ok;
+    if (dirp) {
+
+        /* Close wide-character directory stream */
+        ok = _wclosedir (dirp->wdirp);
+        dirp->wdirp = NULL;
+
+        /* Release multi-byte character version */
+        free (dirp);
+
+    } else {
+
+        /* Invalid directory stream */
+        dirent_set_errno (EBADF);
+        ok = /*failure*/-1;
+
+    }
+    return ok;
+}
+
+/*
+ * Rewind directory stream to beginning.
+ */
+static void
+rewinddir(
+    DIR* dirp)
+{
+    /* Rewind wide-character string directory stream */
+    _wrewinddir (dirp->wdirp);
+}
+
+/*
+ * Scan directory for entries.
+ */
+static int
+scandir(
+    const char *dirname,
+    struct dirent ***namelist,
+    int (*filter)(const struct dirent*),
+    int (*compare)(const void*, const void*))
+{
+    struct dirent **files = NULL;
+    size_t size = 0;
+    size_t allocated = 0;
+    const size_t init_size = 1;
+    DIR *dir = NULL;
+    struct dirent *entry;
+    struct dirent *tmp = NULL;
+    size_t i;
+    int result = 0;
+
+    /* Open directory stream */
+    dir = opendir (dirname);
+    if (dir) {
+
+        /* Read directory entries to memory */
+        while (1) {
+
+            /* Enlarge pointer table to make room for another pointer */
+            if (size >= allocated) {
+                void *p;
+                size_t num_entries;
+
+                /* Compute number of entries in the enlarged pointer table */
+                if (size < init_size) {
+                    /* Allocate initial pointer table */
+                    num_entries = init_size;
+                } else {
+                    /* Double the size */
+                    num_entries = size * 2;
+                }
+
+                /* Allocate first pointer table or enlarge existing table */
+                p = realloc (files, sizeof (void*) * num_entries);
+                if (p != NULL) {
+                    /* Got the memory */
+                    files = (dirent**) p;
+                    allocated = num_entries;
+                } else {
+                    /* Out of memory */
+                    result = -1;
+                    break;
+                }
+
+            }
+
+            /* Allocate room for temporary directory entry */
+            if (tmp == NULL) {
+                tmp = (struct dirent*) malloc (sizeof (struct dirent));
+                if (tmp == NULL) {
+                    /* Cannot allocate temporary directory entry */
+                    result = -1;
+                    break;
+                }
+            }
+
+            /* Read directory entry to temporary area */
+            if (readdir_r (dir, tmp, &entry) == /*OK*/0) {
+
+                /* Did we get an entry? */
+                if (entry != NULL) {
+                    int pass;
+
+                    /* Determine whether to include the entry in result */
+                    if (filter) {
+                        /* Let the filter function decide */
+                        pass = filter (tmp);
+                    } else {
+                        /* No filter function, include everything */
+                        pass = 1;
+                    }
+
+                    if (pass) {
+                        /* Store the temporary entry to pointer table */
+                        files[size++] = tmp;
+                        tmp = NULL;
+
+                        /* Keep up with the number of files */
+                        result++;
+                    }
+
+                } else {
+
+                    /*
+                     * End of directory stream reached => sort entries and
+                     * exit.
+                     */
+                    qsort (files, size, sizeof (void*), compare);
+                    break;
+
+                }
+
+            } else {
+                /* Error reading directory entry */
+                result = /*Error*/ -1;
+                break;
+            }
+
+        }
+
+    } else {
+        /* Cannot open directory */
+        result = /*Error*/ -1;
+    }
+
+    /* Release temporary directory entry */
+    if (tmp) {
+        free (tmp);
+    }
+
+    /* Release allocated memory on error */
+    if (result < 0) {
+        for (i = 0; i < size; i++) {
+            free (files[i]);
+        }
+        free (files);
+        files = NULL;
+    }
+
+    /* Close directory stream */
+    if (dir) {
+        closedir (dir);
+    }
+
+    /* Pass pointer table to caller */
+    if (namelist) {
+        *namelist = files;
+    }
+    return result;
+}
+
+/* Alphabetical sorting */
+static int
+alphasort(
+    const struct dirent **a, const struct dirent **b)
+{
+    return strcoll ((*a)->d_name, (*b)->d_name);
+}
+
+/* Sort versions */
+static int
+versionsort(
+    const struct dirent **a, const struct dirent **b)
+{
+    /* FIXME: implement strverscmp and use that */
+    return alphasort (a, b);
+}
+
+
+/* Convert multi-byte string to wide character string */
+static int
+dirent_mbstowcs_s(
+    size_t *pReturnValue,
+    wchar_t *wcstr,
+    size_t sizeInWords,
+    const char *mbstr,
+    size_t count)
+{
+    int error;
+
+#if defined(_MSC_VER)  &&  _MSC_VER >= 1400
+
+    /* Microsoft Visual Studio 2005 or later */
+    error = mbstowcs_s (pReturnValue, wcstr, sizeInWords, mbstr, count);
+
+#else
+
+    /* Older Visual Studio or non-Microsoft compiler */
+    size_t n;
+
+    /* Convert to wide-character string (or count characters) */
+    n = mbstowcs (wcstr, mbstr, sizeInWords);
+    if (!wcstr  ||  n < count) {
+
+        /* Zero-terminate output buffer */
+        if (wcstr  &&  sizeInWords) {
+            if (n >= sizeInWords) {
+                n = sizeInWords - 1;
+            }
+            wcstr[n] = 0;
+        }
+
+        /* Length of resulting multi-byte string WITH zero terminator */
+        if (pReturnValue) {
+            *pReturnValue = n + 1;
+        }
+
+        /* Success */
+        error = 0;
+
+    } else {
+
+        /* Could not convert string */
+        error = 1;
+
+    }
+
+#endif
+
+    return error;
+}
+
+/* Convert wide-character string to multi-byte string */
+static int
+dirent_wcstombs_s(
+    size_t *pReturnValue,
+    char *mbstr,
+    size_t sizeInBytes, /* max size of mbstr */
+    const wchar_t *wcstr,
+    size_t count)
+{
+    int error;
+
+#if defined(_MSC_VER)  &&  _MSC_VER >= 1400
+
+    /* Microsoft Visual Studio 2005 or later */
+    error = wcstombs_s (pReturnValue, mbstr, sizeInBytes, wcstr, count);
+
+#else
+
+    /* Older Visual Studio or non-Microsoft compiler */
+    size_t n;
+
+    /* Convert to multi-byte string (or count the number of bytes needed) */
+    n = wcstombs (mbstr, wcstr, sizeInBytes);
+    if (!mbstr  ||  n < count) {
+
+        /* Zero-terminate output buffer */
+        if (mbstr  &&  sizeInBytes) {
+            if (n >= sizeInBytes) {
+                n = sizeInBytes - 1;
+            }
+            mbstr[n] = '\0';
+        }
+
+        /* Length of resulting multi-bytes string WITH zero-terminator */
+        if (pReturnValue) {
+            *pReturnValue = n + 1;
+        }
+
+        /* Success */
+        error = 0;
+
+    } else {
+
+        /* Cannot convert string */
+        error = 1;
+
+    }
+
+#endif
+
+    return error;
+}
+
+/* Set errno variable */
+static void
+dirent_set_errno(
+    int error)
+{
+#if defined(_MSC_VER)  &&  _MSC_VER >= 1400
+
+    /* Microsoft Visual Studio 2005 and later */
+    _set_errno (error);
+
+#else
+
+    /* Non-Microsoft compiler or older Microsoft compiler */
+    errno = error;
+
+#endif
+}
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif /*DIRENT_H*/
\ No newline at end of file
diff --git a/Analysis.cpp b/Analysis.cpp
index 2acc89786525b60e22cb2dbd1402069fa9046a8b..748931aeea06199c2f7f1bb7cb9055551e6f99f9 100644
--- a/Analysis.cpp
+++ b/Analysis.cpp
@@ -71,10 +71,10 @@ Analysis::Analysis()
      _building = NULL;
      _projectRootDir="";
      _deltaF=5;   // half of the time interval that used to calculate instantaneous velocity of ped i. Here v_i = (X(t+deltaF) - X(t+deltaF))/(2*deltaF).   X is location.
-     _DoesUseMethodA = false; 						// Method A (Zhang2011a)
-     _DoesUseMethodB = false; 			// Method B (Zhang2011a)
-     _DoesUseMethodC = false; 					// Method C //calculate and save results of classic in separate file
-     _DoesUseMethodD = false;  					// Method D--Voronoi method
+     _DoesUseMethodA = false;                                           // Method A (Zhang2011a)
+     _DoesUseMethodB = false;                   // Method B (Zhang2011a)
+     _DoesUseMethodC = false;                                   // Method C //calculate and save results of classic in separate file
+     _DoesUseMethodD = false;                                   // Method D--Voronoi method
      _cutByCircle = false;  //Adjust whether cut each original voronoi cell by a circle
      _getProfile = false;   // Whether make field analysis or not
      _outputGraph = false;   // Whether output the data for plot the fundamental diagram each frame
@@ -119,27 +119,6 @@ void Analysis::InitArgs(ArgumentParser* args)
 {
      string s = "Parameter:\n";
 
-     switch (args->GetLog()) {
-     case 0:
-          // no log file
-          //Log = new OutputHandler();
-          break;
-     case 1:
-          if(Log) delete Log;
-          Log = new STDIOHandler();
-          break;
-     case 2: {
-          char name[CLENGTH]="";
-          sprintf(name,"%s.P0.dat",args->GetErrorLogFile().c_str());
-          if(Log) delete Log;
-          Log = new FileHandler(name);
-     }
-          break;
-     default:
-          Log->Write("Wrong option for Log file!");
-          exit(0);
-     }
-
      if(args->GetIsMethodA()) {
           _DoesUseMethodA = true;
           vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodA();
@@ -188,6 +167,7 @@ void Analysis::InitArgs(ArgumentParser* args)
      _getProfile = args->GetIsGetProfile();
      _outputGraph = args->GetIsOutputGraph();
      _plotGraph = args->GetIsPlotGraph();
+     _plotIndex = args->GetIsPlotIndex();
      _isOneDimensional=args->GetIsOneDimensional();
      _vComponent = args->GetVComponent();
      _IgnoreBackwardMovement =args->GetIgnoreBackwardMovement();
@@ -273,7 +253,9 @@ std::map<int, polygon_2d> Analysis::ReadGeometry(const std::string& geometryFile
      _highVertexY = geo_maxY;
      _lowVertexX = geo_minX;
      _lowVertexY = geo_minY;
-     //cout<<"INFO: \tGeometry polygon is:\t"<<dsv(geoPoly)<<endl;
+     // using boost::geometry::dsv;
+     // cout<<"INFO: \tGeometry polygon is:\t" << dsv(geoPoly[1])<<endl;
+
      return geoPoly;
 }
 
@@ -375,7 +357,7 @@ int Analysis::RunAnalysis(const string& filename, const string& path)
                     Log->Write("INFO:\tSuccess with Method C using measurement area id %d!\n",_areaForMethod_C[i]->_id);
                     if(_plotTimeseriesC[i])
                     {
-                         string parameters_Timeseries="python3 \""+_scriptsLocation+"/_Plot_timeseries_rho_v.py\" -p \""+ _projectRootDir+VORO_LOCATION + "\" -n "+filename+
+                         string parameters_Timeseries="python \""+_scriptsLocation+"/_Plot_timeseries_rho_v.py\" -p \""+ _projectRootDir+VORO_LOCATION + "\" -n "+filename+
                               " -f "+boost::lexical_cast<std::string>(data.GetFps());
                          int res=system(parameters_Timeseries.c_str());
                          Log->Write("INFO:\t time series result: %d ",res);
@@ -403,6 +385,7 @@ int Analysis::RunAnalysis(const string& filename, const string& path)
                method_D.SetGridSize(_grid_size_X, _grid_size_Y);
                method_D.SetOutputVoronoiCellData(_outputGraph);
                method_D.SetPlotVoronoiGraph(_plotGraph);
+               method_D.SetPlotVoronoiIndex(_plotIndex);
                method_D.SetDimensional(_isOneDimensional);
                method_D.SetCalculateProfiles(_getProfile);
                method_D.SetTrajectoriesLocation(path);
@@ -417,7 +400,7 @@ int Analysis::RunAnalysis(const string& filename, const string& path)
                     Log->Write("INFO:\tSuccess with Method D using measurement area id %d!\n",_areaForMethod_D[i]->_id);
                     if(_plotTimeseriesD[i])
                     {
-                         string parameters_Timeseries="python3 \""+_scriptsLocation+"/_Plot_timeseries_rho_v.py\" -p \""+ _projectRootDir+VORO_LOCATION + "\" -n "+filename+
+                         string parameters_Timeseries="python \""+_scriptsLocation+"/_Plot_timeseries_rho_v.py\" -p \""+ _projectRootDir+VORO_LOCATION + "\" -n "+filename+
                               " -f "+boost::lexical_cast<std::string>(data.GetFps());
                          int res=system(parameters_Timeseries.c_str());
                          Log->Write("INFO:\t time series result: %d ",res);
@@ -493,9 +476,5 @@ int Analysis::mkpath(char* file_path, mode_t mode)
      }
      return 0;
 }
-
+// delete
 #endif
-
-
-
-
diff --git a/Analysis.h b/Analysis.h
index 64f127b6645d7c9134b6e660b11a34e40c24d9c9..91397a763cbb72163618a540263cc0b9c5ff8936 100644
--- a/Analysis.h
+++ b/Analysis.h
@@ -2,7 +2,7 @@
  * \file        Analysis.h
  * \date        Oct 10, 2014
  * \version     v0.7
- * \copyright   <2009-2015> Forschungszentrum J�lich GmbH. All rights reserved.
+ * \copyright   <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
  *
  * \section License
  * This file is part of JuPedSim.
@@ -133,7 +133,11 @@ private:
      int _circleEdges;
      bool _getProfile;        // Whether make field analysis or not
      bool _outputGraph;       // Whether output the data for plot the voronoi diagram each frame
-     bool _plotGraph;       // Whether plot the voronoi diagram each frame
+     bool _plotGraph;       // Whether plot the voronoi diagram each
+                            // frame. if (outputGraph==true)
+     bool _plotIndex;       // Whether plot the voronoi diagram each
+                            // frame with index of pedesrians
+                            // if (outputGraph==true and _polotGraph==true)
      std::vector<bool> _plotTimeseriesA;
      std::vector<bool> _plotTimeseriesC;
      std::vector<bool> _plotTimeseriesD;
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f1002596c4e8122fbfc450c874fcfd1dd3dd84c1..7d2f200eedbdfea3658d4ff36e17249b16125b10 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,21 +1,45 @@
 # Change Log
 All notable changes to `jpsreport` will be documented in this file.
-## v0.8.2 [unreleased]
 
+## v0.8.3 [unreleased]
 ### Added
+- Option to output log in a file instead of the screen fe66fa49
+  ```
+  <logfile>log.txt</logfile>
+  ```
+- Output useful debug information like date, git and compiler versions. !6 and discussion in #79 
+- Option to plot Voronoi diagrams with index instead of little blue circles `plot_index`. Use as:
+  ```xml
+  <output_voronoi_cells enabled="true" plot_graphs="true" plot_index="true"/>
+  ```
+- new format of returned polygons `"index | polygon"` 6fa459ad9ffe5a07699c05b655bcf90f114ed635
+- Exit if `framerate` is not found. c1308ef8
 
-- Use the effective distance in method B in case `length_in_movement_direction` is not defined.  (2c321cef) 
+### Changed
+- Pass Matrix by reference bff89e48
+- Better fonts in plot scripts. 56d6a8f7
+
+
+### Fixed
+- Ignore empty line in traj file. 3a3ae04e
+- Fixes for profile plots. a8a1414c
+- Fix parsing of framerate. 2ad0b01d
+
+## v0.8.2 [06.11.2017]
+
+### Added
+
+- Use the effective distance in method B in case `length_in_movement_direction` is not defined.  (2c321cef)
 - Added an error warning when the number of agents in the trajectory is not corresponding to total ids or the ped ids are not continuous in the first frame.
 
 ### Changed
-- Code does not come with Boost anymore. User should install Boost before using jpsreport. (2c0c69f3) 
+- Code does not come with Boost anymore. User should install Boost before using jpsreport. (2c0c69f3)
 - use boost matrix instead of double pointers (9ff5c978)
 - Use own index numbers (9a0d8ec8)
 - Use Python3 in scripts.
 
-### Fixed 
-- Fixed SegFault due to reading files from different OS. (9a42c9dd)
-
+### Fixed
+- Fix SegFault due to reading files from different OS. (9a42c9dd)
 
 ## v0.8.1 [11.10.2016]
 
@@ -23,7 +47,7 @@ All notable changes to `jpsreport` will be documented in this file.
 
 - Two options `startframe` and `stopframe` are added for each measurement area for method D to assign the time periods for analysis.
 
-- Individual density based on Voronoi method is added for one dimensional case in the output file (Individual headway is moved to the 5th column). 
+- Individual density based on Voronoi method is added for one dimensional case in the output file (Individual headway is moved to the 5th column).
 
 - z-position of each measurement area can be assigned in inifile so that the trajectories in geometries with several floors can be analyzed.
 
@@ -66,7 +90,7 @@ All notable changes to `jpsreport` will be documented in this file.
 
 - Issue a warning when the voronoi cell cannot be calculated.
 
-- A warning will will be given and the program stops if trajectory for a given pedestrian ID is not continuous. 
+- A warning will will be given and the program stops if trajectory for a given pedestrian ID is not continuous.
 
 
 ### Changed
@@ -93,7 +117,7 @@ All notable changes to `jpsreport` will be documented in this file.
 
 - when path of trajectory is not given absolutely, the default location is the same folder with the inifile
 
-	
+
 ## v0.7
 
 ### Added
@@ -108,38 +132,38 @@ All notable changes to `jpsreport` will be documented in this file.
 
 - Changed name of some variables in configuration file:
 
-	**measurementAreas**                --->  **measurement_areas**
-    
-	**Length_in_movement_direction**	---> **length_in_movement_direction**
-	
-	**useXComponent**		            ---> **use_x_component**
-	
-	**useYComponent**		            ---> **use_y_component**
-	
-	**halfFrameNumberToUse**            ---> **frame_step**
-	
-	**timeInterval**	                ---> **frame_interval**
-	
-	**measurementArea**	                ---> **measurement_area**
-	
-	**outputGraph**	                    ---> **output_graph**
-	
-	**individualFDdata**	            ---> **individual_FD**
-	
-	**cutByCircle** 	                ---> **cut_by_circle**
-	
-	**getProfile** 		                ---> **profiles**
-	
-	**scale_x**			                ---> **grid_size_x**
-	
-	**scale_y**			                ---> **grid_size_y**
+    **measurementAreas**                --->  **measurement_areas**
+
+    **Length_in_movement_direction**	---> **length_in_movement_direction**
+
+    **useXComponent**                   ---> **use_x_component**
+
+    **useYComponent**                   ---> **use_y_component**
+
+    **halfFrameNumberToUse**            ---> **frame_step**
+
+    **timeInterval**                    ---> **frame_interval**
+
+    **measurementArea**                 ---> **measurement_area**
+
+    **outputGraph**                     ---> **output_graph**
+
+    **individualFDdata**                ---> **individual_FD**
+
+    **cutByCircle**                     ---> **cut_by_circle**
+
+    **getProfile**                      ---> **profiles**
+
+    **scale_x**                         ---> **grid_size_x**
+
+    **scale_y**                         ---> **grid_size_y**
 - Changed the data type of frame rate (fps) from integer to float
 
-- Changed the way for dealing with pedestrian outside geometry. In old version JPSreport stops when some pedestrians are outside geometry but now it continue working by 
+- Changed the way for dealing with pedestrian outside geometry. In old version JPSreport stops when some pedestrians are outside geometry but now it continue working by
 removing these pedestrians from the list.
 
 - More than one sub rooms in one geometry can be analysed independently.
-	
+
 ### Fixed
-	
-- Fixed bug for dealing with obstacles inside geometry.
\ No newline at end of file
+
+- Fixed bug for dealing with obstacles inside geometry.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 02398ae060665e55f5c1ae9746456ca7273e88a6..6fcb0358aea9c3819ca04d91691fa71b871593ca 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,6 +20,56 @@ if (Boost_NO_SYSTEM_PATHS)
   set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH} ${BOOST_LIBRARY_DIRS})
 endif (Boost_NO_SYSTEM_PATHS)
 
+
+find_package(Git REQUIRED) # no need for this msg. It comes from cmake.findgit()
+
+find_program(GIT_SCM git DOC "Git version control")
+mark_as_advanced(GIT_SCM)
+find_file(GITDIR NAMES .git PATHS ${CMAKE_SOURCE_DIR} NO_DEFAULT_PATH)
+if (GIT_SCM AND GITDIR)
+
+# the commit's SHA1, and whether the building workspace was dirty or not
+# describe --match=NeVeRmAtCh --always --tags --abbrev=40 --dirty
+execute_process(COMMAND
+  "${GIT_EXECUTABLE}" --no-pager describe --tags --always --dirty
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+  OUTPUT_VARIABLE GIT_SHA1
+  ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+# branch
+execute_process(
+  COMMAND "${GIT_EXECUTABLE}" rev-parse --abbrev-ref HEAD
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+  OUTPUT_VARIABLE GIT_BRANCH
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+# the date of the commit
+execute_process(COMMAND
+  "${GIT_EXECUTABLE}" log -1 --format=%ad --date=local
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+  OUTPUT_VARIABLE GIT_DATE
+  ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+# the subject of the commit
+execute_process(COMMAND
+  "${GIT_EXECUTABLE}" log -1 --format=%s
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+  OUTPUT_VARIABLE GIT_COMMIT_SUBJECT
+  ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+
+
+add_definitions("-DGIT_COMMIT_HASH=\"${GIT_SHA1}\"")
+add_definitions("-DGIT_COMMIT_DATE=\"${GIT_DATE}\"")
+add_definitions("-DGIT_COMMIT_SUBJECT=\"${GIT_COMMIT_SUBJECT}\"")
+add_definitions("-DGIT_BRANCH=\"${GIT_BRANCH}\"")
+
+else()
+    message(STATUS "Not in a git repo")
+endif()
+
+
 #set(EXECUTABLE_OUTPUT_PATH "../")
 #INCLUDE_DIRECTORIES("./")
 # message( STATUS "CMAKE_BINARY_DIR: " ${CMAKE_BINARY_DIR} )
@@ -29,12 +79,13 @@ set(CMAKE_TEST_DIR ${CMAKE_SOURCE_DIR}/Utest)
 
 set(JPSREPORT_MAJOR_VERSION 0)
 set(JPSREPORT_MINOR_VERSION 8)
-set(JPSREPORT_PATCH_VERSION 2)
+set(JPSREPORT_PATCH_VERSION 3)
 set(JPSREPORT_VERSION
   ${JPSREPORT_MAJOR_VERSION}.${JPSREPORT_MINOR_VERSION}.${JPSREPORT_PATCH_VERSION})
 message( STATUS "JPSREPORT_VERSION: " ${JPSREPORT_VERSION} )
+add_definitions("-DJPSREPORT_VERSION=\"${JPSREPORT_VERSION}\"")
 if(NOT CMAKE_BUILD_TYPE)
-  set (CMAKE_BUILD_TYPE Release) 
+  set (CMAKE_BUILD_TYPE Release)
 endif(NOT CMAKE_BUILD_TYPE)
 message( STATUS "CMAKE_BUILD_TYPE: " ${CMAKE_BUILD_TYPE} )
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/bin")
@@ -74,7 +125,7 @@ endif(NOT DEFINED PROCESSOR_COUNT)
 
 if(PROCESSOR_COUNT)
   # add 1 should be magic! http://www.kitware.com/blog/home/post/63
-  #math(EXPR PROCESSOR_COUNT "${PROCESSOR_COUNT} + 1") 
+  #math(EXPR PROCESSOR_COUNT "${PROCESSOR_COUNT} + 1")
   message( STATUS "PROCESSOR_COUNT " ${PROCESSOR_COUNT})
   set(CTEST_BUILD_FLAGS "-j${PROCESSOR_COUNT}")
 endif(PROCESSOR_COUNT)
@@ -102,7 +153,6 @@ if(BUILD_TESTING)
 endif(BUILD_TESTING)
 
 set(source_files
-  getRSS.c
   Analysis.cpp
   IO/OutputHandler.cpp
   general/ArgumentParser.cpp
@@ -116,45 +166,45 @@ set(source_files
   methods/Method_B.cpp
   methods/Method_C.cpp
   methods/Method_D.cpp
-  geometry/Building.cpp  
-  geometry/Line.cpp      
-  geometry/Point.cpp    
+  geometry/Building.cpp
+  geometry/Line.cpp
+  geometry/Point.cpp
   geometry/Transition.cpp
-  geometry/Crossing.cpp  
-  geometry/NavLine.cpp   
-  geometry/Room.cpp     
+  geometry/Crossing.cpp
+  geometry/NavLine.cpp
+  geometry/Room.cpp
   geometry/Wall.cpp
-  geometry/Hline.cpp     
-  geometry/Obstacle.cpp  
+  geometry/Hline.cpp
+  geometry/Obstacle.cpp
   geometry/SubRoom.cpp
   geometry/Goal.cpp
-)  
+)
 set (  header_files
   Analysis.h
-  methods/MeasurementArea.h   
+  methods/MeasurementArea.h
   methods/VoronoiDiagram.h
   methods/PedData.h
   methods/Method_A.h
   methods/Method_B.h
   methods/Method_C.h
-  methods/Method_D.h  
+  methods/Method_D.h
   IO/OutputHandler.h
   general/ArgumentParser.h
   general/Macros.h
   tinyxml/tinyxml.h
   tinyxml/tinystr.h
-  geometry/Building.h  
-  geometry/Line.h      
-  geometry/Point.h    
+  geometry/Building.h
+  geometry/Line.h
+  geometry/Point.h
   geometry/Transition.h
-  geometry/Crossing.h  
-  geometry/NavLine.h   
-  geometry/Room.h     
+  geometry/Crossing.h
+  geometry/NavLine.h
+  geometry/Room.h
   geometry/Wall.h
-  geometry/Hline.h     
-  geometry/Obstacle.h  
+  geometry/Hline.h
+  geometry/Obstacle.h
   geometry/SubRoom.h
-  geometry/Goal.h  
+  geometry/Goal.h
   )
 
 
@@ -205,19 +255,25 @@ message( STATUS "Boost_INCLUDE_DIR: " ${Boost_INCLUDE_DIR} )
 message( STATUS "Boost_LIBRARY_DIR: " ${Boost_LIBRARY_DIR} )
 
 include_directories(${Boost_INCLUDE_DIR})
-link_directories(${Boost_LIBRARY_DIR})
+#link_directories(${Boost_LIBRARY_DIR})
+
+if (WIN32)
+find_library (PSAPI Psapi PATH_SUFFIXES "x64")
+message (STATUS "PSAPI: ${PSAPI}")
+endif()
 
 add_library ( geometrycore STATIC ${source_files} )
 add_executable(
-	jpsreport main.cpp
+        jpsreport main.cpp
 )
-target_link_libraries( jpsreport ${Boost_LIBRARIES} )
-target_link_libraries( jpsreport geometrycore )
+
 
 if(WIN32)
 target_link_libraries (jpsreport wsock32)
+target_link_libraries( jpsreport ${PSAPI} )
 endif()
 
+target_link_libraries( jpsreport geometrycore )
 
 # ----------------------------- cTest ------------------------------------------
 if(BUILD_TESTING)
@@ -228,11 +284,11 @@ if(BUILD_TESTING)
   include(CTest) #adding Dart support
 
   #test if code compiles and runs default setting. Takes about 30 seconds
-  add_test (jpsreport_compile ${CMAKE_CTEST_COMMAND} 
+  add_test (jpsreport_compile ${CMAKE_CTEST_COMMAND}
   --build-and-test "${CMAKE_SOURCE_DIR}" "${EXECUTABLE_OUTPUT_PATH}" #"${CMAKE_BINARY_DIR}"
   --build-generator ${CMAKE_GENERATOR}
   --build-makeprogram ${CMAKE_MAKE_PROGRAM} -j${PROCESSOR_COUNT}
-  --build-two-config  
+  --build-two-config
   --build-exe-dir ${EXECUTABLE_OUTPUT_PATH}  # todo wo soll der exe hin?: ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
   --build-project JPSreport
   --test-command jpsreport --ini=${CMAKE_TEST_DIR}/files/input_UO_180.xml
diff --git a/demos/HowTo.ipynb b/demos/HowTo.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..fc1e0c99e3a205f9b0c9fa08f0cd122654d27eb7
--- /dev/null
+++ b/demos/HowTo.ipynb
@@ -0,0 +1,339 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Getting started with `JPSreport`\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`JPSreport` needs three different files as input: \n",
+    "- A trajectory file (txt or xml) (see [documentation](http://www.jupedsim.org/jpsreport/2016-11-03-trajectory))\n",
+    "- A geometry file (see [documentation](http://www.jupedsim.org/jpsreport/2016-11-02-geometry)) \n",
+    "- and a project file, called inifile (see [documentation](http://www.jupedsim.org/jpsreport/2016-11-01-inifile))\n",
+    "\n",
+    "For example in **demos/bottleneck** we can find the following files:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Calling `JPSreport` without an inifile yields:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "INFO: \tTrying to load the default configuration from the file <ini.xml>\n",
+      "INFO: \tParsing the ini file <ini.xml>\n",
+      "Usage: \n",
+      "\n",
+      "\t../bin/jpsreport input.xml\n",
+      "\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "ERROR: \tFailed to open file\n",
+      "ERROR: \tCould not parse the ini file\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%bash\n",
+    "../bin/jpsreport"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`JPSreport` rightly complains about a missing inifile. Now, let's have a look in the directory `/demos/bottleneck`: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "geo_AO_300.xml   geometry.png     ini_AO_300.xml   traj_AO_300.txt\r\n"
+     ]
+    }
+   ],
+   "source": [
+    " ls ./bottleneck"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here we can find the three necessary files to get us started. \n",
+    "Let's visualize the geometry with [JPSvis](https://gitlab.version.fz-juelich.de/jupedsim/jpsvis)\n",
+    "\n",
+    "![geometry](./bottleneck/geometry.png)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Preparing the inifile\n",
+    "\n",
+    "here do whatever you want"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Call JPSreport\n",
+    "\n",
+    "To run the analysis run from a terminal the following:\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "INFO: \tParsing the ini file <./bottleneck/ini_AO_300.xml>\n",
+      "INFO: \tGeometry File is: <./bottleneck/geo_AO_300.xml>\n",
+      "INFO: \tFormat of the trajectory file is: <.txt>\n",
+      "INFO: \tInput trajectory file is\t<traj_AO_300.txt>\n",
+      "./\n",
+      "INFO: \tInput directory for loading trajectory is:\t<./bottleneck/./>\n",
+      "INFO: \tInput directory for loading scripts is:\t<./bottleneck/../../scripts/>\n",
+      "INFO: \tMeasure area id  <1> with type <BoundingBox>\n",
+      "\t\tMeasure area points  < 2.400, 0.530>\n",
+      "\t\tMeasure area points  < 2.400, -0.530>\n",
+      "\t\tMeasure area points  < -0.600, -0.530>\n",
+      "\t\tMeasure area points  < -0.600, 0.530>\n",
+      "\t\tLength in movement direction 1.000\n",
+      "INFO: \tMeasure area id  <2> with type <Line>\n",
+      "\t\tMeasurement line starts from  <-2.250, 0.000> to <4.000, 0.000>\n",
+      "INFO: \tMeasure area id  <4> with type <Line>\n",
+      "\t\tMeasurement line starts from  <-2.250, 0.500> to <4.000, 0.500>\n",
+      "INFO: \tBoth x and y-component of coordinates will be used to calculate instantaneous velocity over <10 frames>\n",
+      "INFO: \tMethod A is selected\n",
+      "INFO: \tMeasurement area id <2> will be used for analysis\n",
+      "\tFrame interval used for calculating flow is <100> frame\n",
+      "\tThe Time series N-t measured will be plotted!! \n",
+      "INFO: \tMeasurement area id <4> will be used for analysis\n",
+      "\tFrame interval used for calculating flow is <150> frame\n",
+      "\tThe Time series N-t measured will be plotted!! \n",
+      "INFO: \tFinish parsing inifile\n",
+      "\n",
+      "INFO: \tStart Analysis for the file: traj_AO_300.txt\n",
+      "**********************************************************************\n",
+      "INFO: \tLoading building file successful!!!\n",
+      "\n",
+      "INFO: \tInit Geometry\n",
+      "INFO: \tInit Geometry successful!!!\n",
+      "\n",
+      "INFO:\tthe name of the trajectory is: <traj_AO_300.txt>\n",
+      "INFO:\tfull name of the trajectory is: <./bottleneck/.//traj_AO_300.txt>\n",
+      "INFO:\tFrame rate fps: <16.00>\n",
+      "INFO: pos_id: 0\n",
+      "INFO: pos_fr: 1\n",
+      "INFO: pos_x: 2\n",
+      "INFO: pos_y: 3\n",
+      "INFO: pos_z: 4\n",
+      "INFO: pos_vd: 5\n",
+      "lineNr 100000\n",
+      "INFO:\t Finished reading the data\n",
+      "INFO: Got 128783 lines\n",
+      "INFO: minID: 1\n",
+      "INFO: minFrame: 0\n",
+      "INFO: numFrames: 959\n",
+      "INFO: Total number of Agents: 348\n",
+      "INFO: Enter CreateGlobalVariables with numPeds=348 and numFrames=959\n",
+      "INFO: allocate memory for xCor\n",
+      "INFO: allocate memory for yCor\n",
+      "INFO: allocate memory for zCor\n",
+      "INFO: allocate memory for vComp\n",
+      " Finished memory allocation\n",
+      "INFO: Leave CreateGlobalVariables()\n",
+      "INFO: Create Global Variables done\n",
+      "convert x and y\n",
+      "Save the data for each frame\n",
+      "------------------------Analyzing with Method A-----------------------------\n",
+      "frame ID = 0\n",
+      "------------------------Analyzing with Method A-----------------------------\n",
+      "frame ID = 0\n",
+      "frame ID = 100\n",
+      "frame ID = 100\n",
+      "frame ID = 200\n",
+      "frame ID = 200\n",
+      "frame ID = 300\n",
+      "frame ID = 300\n",
+      "frame ID = 400\n",
+      "frame ID = 400\n",
+      "frame ID = 500\n",
+      "frame ID = 500\n",
+      "frame ID = 600\n",
+      "frame ID = 600\n",
+      "frame ID = 700\n",
+      "frame ID = 800\n",
+      "frame ID = 900\n",
+      "frame ID = 700\n",
+      "frame ID = 800\n",
+      "frame ID = 900\n",
+      "INFO:\tPlotting N-t diagram! Status: 0\n",
+      "INFO:\tSuccess with Method A using measurement area id 4!\n",
+      "\n",
+      "INFO:\tPlotting N-t diagram! Status: 0\n",
+      "INFO:\tSuccess with Method A using measurement area id 2!\n",
+      "\n",
+      "**********************************************************************\n",
+      "INFO: \tEnd Analysis for the file: traj_AO_300.txt\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%bash\n",
+    "../bin/jpsreport ./bottleneck/ini_AO_300.xml"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The output of `JPSreport` gives some logging information on its running process.  It may inform on any errors during the compilation or misconception of the input files (inifile, geometry or trajectory file).\n",
+    "\n",
+    "Uppon a succesful run, the directory **output** is created with the following content (depends on the configuration in inifile):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "bottleneck/Output\n",
+      "└── Fundamental_Diagram\n",
+      "    └── FlowVelocity\n",
+      "        ├── FDFlowVelocity_traj_AO_300.txt_id_2.dat\n",
+      "        ├── FDFlowVelocity_traj_AO_300.txt_id_4.dat\n",
+      "        ├── Flow_NT_traj_AO_300.txt_id_2.dat\n",
+      "        ├── Flow_NT_traj_AO_300.txt_id_2.png\n",
+      "        ├── Flow_NT_traj_AO_300.txt_id_4.dat\n",
+      "        └── Flow_NT_traj_AO_300.txt_id_4.png\n",
+      "\n",
+      "2 directories, 6 files\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%bash\n",
+    "tree bottleneck/Output"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "#Flow rate(1/s)\t\t Mean velocity(m/s)\n",
+      "8.404\t1.148\n",
+      "7.673\t0.869\n",
+      "7.347\t0.774\n",
+      "7.216\t0.723\n",
+      "7.129\t0.736\n",
+      "6.880\t0.726\n",
+      "6.275\t0.740\n",
+      "3.918\t0.748\n",
+      "2.500\t0.806\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%bash\n",
+    "tail bottleneck/Output/Fundamental_Diagram/FlowVelocity/FDFlowVelocity_traj_AO_300.txt_id_2.dat"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![](bottleneck/Output/Fundamental_Diagram/FlowVelocity/Flow_NT_traj_AO_300.txt_id_2.png)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "![](Output/Fundamental_Diagram/FlowVelocity/Flow_NT_traj_AO_300.txt_id_2.png)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Result: \n",
+    "![](bottleneck/Output/Fundamental_Diagram/FlowVelocity/Flow_NT_traj_AO_300.txt_id_2.png)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/demos/T-Junction/ini_KO_240_050_240.xml b/demos/T-Junction/ini_KO_240_050_240.xml
index 88fb87f3415297102098c0a0b0afc95fd275462d..67f5ae75d4c0e7944224ee0a6f9c31ca76525b22 100644
--- a/demos/T-Junction/ini_KO_240_050_240.xml
+++ b/demos/T-Junction/ini_KO_240_050_240.xml
@@ -5,7 +5,7 @@
 	<!-- trajectories file and format -->
 	<!-- either a file name or a path location. In the latter case all files in the directory will be used-->
 	<trajectories format="txt" unit="m">
-		<file name="traj_KO_240_050_240_x.txt" />
+		<file name="traj_KO_240_050_240.txt" />
 		<path location="./" />  
 	</trajectories>
 	<!-- give relative path based on the location inifile or give the absolute path- -->
@@ -37,14 +37,14 @@
 		</area_L>
 	</measurement_areas>
 
-<!--	<velocity>
+	<velocity>
 		<use_x_component>true</use_x_component>
 		<use_y_component>true</use_y_component>
         		<!-- The time interval that used to calculate instantaneous velocity 
-			of ped i [fr] here v_i = (X(t+frame_step/2) - X(t+frame_step/2))/frame_step. X is location. 
+			of ped i [fr] here v_i = (X(t+frame_step/2) - X(t+frame_step/2))/frame_step. X is location. -->
 		<frame_step>10</frame_step>
 	</velocity>
--->   
+
     <velocity frame_step="10" set_movement_direction="None" ignore_backward_movement="false"/>
     <!-- frame_step is the time interval that used to calculate instantaneous velocity 
 	of ped i [fr] here v_i = (X(t+frame_step/2) - X(t+frame_step/2))/frame_step. X is location. -->
@@ -62,19 +62,18 @@
 	</method_B>
 
 	<!-- Method C (Zhang2011a) Classical density and Vel -->
-	<method_C enabled="true">
+	<method_C enabled="false">
 		<measurement_area id="1" plot_time_series="true"/>
         <measurement_area id="2" plot_time_series="true"/>
 	</method_C>
 	
 	<!-- Method D (Zhang2011a) Voronoi density and Vel -->
 	<method_D enabled="true"> 
-        <measurement_area id="1" start_frame="None" stop_frame="None" get_individual_FD="false" plot_time_series="true"/> 
-        <measurement_area id="2" start_frame="500" stop_frame="800" get_individual_FD="false" plot_time_series="true"/> 
+        <measurement_area id="1" start_frame="500" stop_frame="800" get_individual_FD="false" plot_time_series="true"/> 
 		<one_dimensional enabled="false"/>
         <cut_by_circle enabled="false" radius="1.0" edges="10"/>
         <output_voronoi_cells enabled="false" plot_graphs="false"/>
-        <profiles enabled="false" grid_size_x="0.20" grid_size_y="0.20"/> 
+        <profiles enabled="true" grid_size_x="0.20" grid_size_y="0.20"/> 
     </method_D> 
 </JPSreport>
 
diff --git a/demos/bottleneck/geometry.png b/demos/bottleneck/geometry.png
new file mode 100644
index 0000000000000000000000000000000000000000..4e3b5a95671d64fe4b72d5b66e618f8a80c5345d
Binary files /dev/null and b/demos/bottleneck/geometry.png differ
diff --git a/demos/bottleneck/ini_AO_300.xml b/demos/bottleneck/ini_AO_300.xml
index 5410aebfe8ca1f43715f00c3bf986619e100766c..7784b57ec9ebc029bc55a169e9604d2f49e990da 100644
--- a/demos/bottleneck/ini_AO_300.xml
+++ b/demos/bottleneck/ini_AO_300.xml
@@ -1,62 +1,43 @@
-<?xml version="1.0" encoding="UTF-8"?> 
-<JPSreport project="JPS-Project" version="0.8" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="../../xsd/jps_report.xsd"> 
-    <!-- geometry file --> 
-    <geometry file = "geo_AO_300.xml" /> 
-    <!-- trajectories file and format --> 
-    <!-- either a file name or a path location. In the latter case all files in the directory will be used--> 
-    <trajectories format="txt" unit="m"> 
-        <file name="traj_AO_300.txt" /> 
-        <path location="./" /> 
-    </trajectories> 
-	<!----give relative path based on the location inifile or give the absolute path--->
-	<scripts location="../../scripts/"/> 
-	
-    <measurement_areas unit="m"> 
-        <area_B id="1" type="BoundingBox" zPos="None"> 
-            <vertex x="2.40" y="0.53" /> 
-            <vertex x="2.40" y="-0.53" /> 
-            <vertex x="-0.60" y="-0.53" /> 
-            <vertex x="-0.60" y="0.53" /> 
-            <length_in_movement_direction distance="1.0"/> 
-        </area_B> 
-        <area_L id="2" type="Line" zPos="None"> 
-            <start x="-2.25" y="0.00" /> 
-            <end x="4.00" y="0.00" /> 
-        </area_L> 
-	<area_L id="4" type="Line" zPos="None">
+<?xml version="1.0" encoding="UTF-8"?>
+<JPSreport project="JPS-Project" version="0.8" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="../../xsd/jps_report.xsd">
+    <!-- geometry file -->
+    <geometry file = "geo_AO_300.xml" />
+    <!-- trajectories file and format -->
+    <!-- either a file name or a path location. In the latter case all files in the directory will be used-->
+    <trajectories format="txt" unit="m">
+        <file name="traj_AO_300.txt" />
+        <path location="./" />
+    </trajectories>
+        <!-- give relative path based on the location inifile or give the absolute path -->
+        <scripts location="../../scripts/"/>
+
+    <measurement_areas unit="m">
+        <area_B id="1" type="BoundingBox" zPos="None">
+            <vertex x="2.40" y="0.53" />
+            <vertex x="2.40" y="-0.53" />
+            <vertex x="-0.60" y="-0.53" />
+            <vertex x="-0.60" y="0.53" />
+            <length_in_movement_direction distance="1.0"/>
+        </area_B>
+        <area_L id="2" type="Line" zPos="None">
+            <start x="-2.25" y="0.00" />
+            <end x="4.00" y="0.00" />
+        </area_L>
+        <area_L id="4" type="Line" zPos="None">
             <start x="-2.25" y="0.50" />
             <end x="4.00" y="0.50" />
         </area_L>
-    </measurement_areas> 
+    </measurement_areas>
 
     <velocity frame_step="10" set_movement_direction="None" ignore_backward_movement="false"/>
-    <!-- frame_step is the time interval that used to calculate instantaneous velocity 
-	of ped i [fr] here v_i = (X(t+frame_step/2) - X(t+frame_step/2))/frame_step. X is location. -->
-
-    <!-- Method A (Zhang2011a) Flow and Vel --> 
-    <method_A enabled="true"> 
-        <!-- Time interval used to count the flow [fr] --> 
-        <!-- The coordinate of the line used to calculate the flow and velocity --> 
-        <measurement_area id="2" frame_interval="100" plot_time_series="true"/> 
-		<measurement_area id="4" frame_interval="150" plot_time_series="true"/>
-    </method_A> 
-
-    <!-- Method B (Zhang2011a) Vel and Dens based on Tin and Tout --> 
-    <method_B enabled="false"> 
-        <measurement_area id="1" /> 
-    </method_B> 
-
-    <!-- Method C (Zhang2011a) Classical density and Vel --> 
-    <method_C enabled="true"> 
-        <measurement_area id="1" plot_time_series="true"/> 
-    </method_C> 
+    <!-- frame_step is the time interval that used to calculate instantaneous velocity  -->
+        <!-- of ped i [fr] here v_i = (X(t+frame_step/2) - X(t+frame_step/2))/frame_step. X is location. -->
 
-    <!-- Method D (Zhang2011a) Voronoi density and Vel --> 
-	<method_D enabled="true"> 
-        <measurement_area id="1" start_frame="None" stop_frame="None" plot_time_series="true" get_individual_FD="false"/>
-		<one_dimensional enabled="false"/>		
-        <cut_by_circle enabled="true" radius="1.0" edges="10"/>
-        <output_voronoi_cells enabled="true" plot_graphs="true"/>
-        <profiles enabled="false" grid_size_x="0.20" grid_size_y="0.20"/> 
-    </method_D> 
-</JPSreport> 
+    <!-- Method A (Zhang2011a) Flow and Vel -->
+    <method_A enabled="true">
+        <!-- Time interval used to count the flow [fr] -->
+        <!-- The coordinate of the line used to calculate the flow and velocity -->
+        <measurement_area id="2" frame_interval="100" plot_time_series="true"/>
+        <measurement_area id="4" frame_interval="150" plot_time_series="true"/>
+    </method_A>
+</JPSreport>
diff --git a/general/ArgumentParser.cpp b/general/ArgumentParser.cpp
index f88102af32b81522dca8c7ec61e3213aaaa645cb..202878589f8e97c65fb53f179a0821a31925b2c2 100644
--- a/general/ArgumentParser.cpp
+++ b/general/ArgumentParser.cpp
@@ -2,7 +2,7 @@
  * \file        ArgumentParser.cpp
  * \date        Oct 10, 2014
  * \version     v0.7
- * \copyright   <2009-2015> Forschungszentrum J��lich GmbH. All rights reserved.
+ * \copyright   <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
  *
  * \section License
  * This file is part of JuPedSim.
@@ -34,8 +34,11 @@
 #include <string>
 #include <sstream>
 #include <math.h>
+#ifdef _MSC_VER
+#include "../.vs/dirent.h"
+#else
 #include <dirent.h>
-
+#endif
 #ifdef _OPENMP
 #include <omp.h>
 #else
@@ -50,11 +53,39 @@
 
 using namespace std;
 
+/* https://stackoverflow.com/questions/38530981/output-compiler-version-in-a-c-program#38531037 */
+std::string ver_string(int a, int b, int c) {
+      std::ostringstream ss;
+      ss << a << '.' << b << '.' << c;
+      return ss.str();
+}
+
+  std::string true_cxx =
+#ifdef __clang__
+   "clang++";
+#else
+   "g++";
+#endif
+
+  std::string true_cxx_ver =
+#ifdef __clang__
+    ver_string(__clang_major__, __clang_minor__, __clang_patchlevel__);
+#else
+    ver_string(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__);
+#endif
+
+// todo: handle Visual Studio
+/* #ifdef _MSC_VER */
+/*     std::to_string(_MSC_VER) */
+/* #endif */
+
+
+
 void ArgumentParser::Usage(const std::string file)
 {
 
      Log->Write("Usage: \n");
-     Log->Write("\t%s input.xml\n",file.c_str());
+     Log->Write("%s inifile.xml\n",file.c_str());
      exit(EXIT_SUCCESS);
 }
 
@@ -73,14 +104,15 @@ ArgumentParser::ArgumentParser()
      _isCutByCircle = false;
      _isOutputGraph= false;
      _isPlotGraph= false;
+     _isPlotIndex = false;
      _isOneDimensional=false;
      _isGetProfile =false;
      _steadyStart =100;
      _steadyEnd = 1000;
      _grid_size_X = 10;
      _grid_size_Y = 10;
-     _errorLogFile="./Logfile.dat";
-     _log=1; //no output wanted
+     _errorLogFile="log.txt";
+     _log=2; //no output wanted
      _trajectoriesLocation="./";
      _trajectoriesFilename="";
      _projectRootDir="./";
@@ -112,11 +144,6 @@ bool ArgumentParser::ParseArgs(int argc, char **argv)
           Usage(argv[0]);
           return false;
      }
-     else if(argument == "-v" || argument == "--version")
-     {
-          fprintf(stderr,"You are actually using JuPedsim (jpsreport) version %s  \n\n",JPS_VERSION);
-          return false;
-     }
 
      // other special case where a single configuration file is submitted
      //check if inifile options are given
@@ -151,7 +178,14 @@ const string& ArgumentParser::GetProjectRootDir() const
 
 bool ArgumentParser::ParseIniFile(const string& inifile)
 {
-
+     // first logs will go to stdout
+     Log->Write("----\nJuPedSim - JPSreport\n");
+     Log->Write("Current date   : %s %s", __DATE__, __TIME__);
+     Log->Write("Version        : %s", JPSREPORT_VERSION);
+     Log->Write("Compiler       : %s (%s)", true_cxx.c_str(), true_cxx_ver.c_str());
+     Log->Write("Commit hash    : %s", GIT_COMMIT_HASH);
+     Log->Write("Commit date    : %s", GIT_COMMIT_DATE);
+     Log->Write("Branch         : %s\n----\n", GIT_BRANCH);
      Log->Write("INFO: \tParsing the ini file <%s>",inifile.c_str());
 
      //extract and set the project root dir
@@ -180,6 +214,45 @@ bool ArgumentParser::ParseIniFile(const string& inifile)
           return false;
      }
 
+     if (xMainNode->FirstChild("logfile")) {
+          this->SetErrorLogFile(
+               this->GetProjectRootDir()+xMainNode->FirstChild("logfile")->FirstChild()->Value());
+          this->SetLog(2);
+          Log->Write("INFO:\tlogfile <%s>", this->GetErrorLogFile().c_str());
+     }
+     switch (this->GetLog()) {
+     case 0:
+          // no log file
+          //Log = new OutputHandler();
+          break;
+     case 1:
+          if(Log) delete Log;
+          Log = new STDIOHandler();
+          break;
+     case 2: {
+          char name[CLENGTH]="";
+          sprintf(name,"%s", this->GetErrorLogFile().c_str());
+          if(Log) delete Log;
+          Log = new FileHandler(name);
+     }
+          break;
+     default:
+          Log->Write("ERROR: \tWrong option for Log file!");
+          exit(0);
+     }
+     // from this point if case 2, the logs will go to a logfile
+     if(this->GetLog() == 2)
+     {
+          Log->Write("----\nJuPedSim - JPSreport\n");
+          Log->Write("Current date   : %s %s", __DATE__, __TIME__);
+          Log->Write("Version        : %s", JPSREPORT_VERSION);
+          Log->Write("Compiler       : %s (%s)", true_cxx.c_str(), true_cxx_ver.c_str());
+          Log->Write("Commit hash    : %s", GIT_COMMIT_HASH);
+          Log->Write("Commit date    : %s", GIT_COMMIT_DATE);
+          Log->Write("Branch         : %s\n----\n", GIT_BRANCH);
+     }
+
+
      //geometry
      if(xMainNode->FirstChild("geometry"))
      {
@@ -373,7 +446,7 @@ bool ArgumentParser::ParseIniFile(const string& inifile)
                }
                correct(poly); // in the case the Polygone is not closed
                areaB->_poly=poly;
-               
+
                TiXmlElement* xLength=xMeasurementArea_B->FirstChildElement("length_in_movement_direction");
                if(xLength)
                {
@@ -413,7 +486,7 @@ bool ArgumentParser::ParseIniFile(const string& inifile)
                Log->Write("\t\tMeasurement line starts from  <%.3f, %.3f> to <%.3f, %.3f>",areaL->_lineStartX*CMtoM,areaL->_lineStartY*CMtoM,areaL->_lineEndX*CMtoM,areaL->_lineEndY*CMtoM);
           }
      }
-     
+
      //instantaneous velocity
      /*    TiXmlNode* xVelocity=xMainNode->FirstChild("velocity");
            if(xVelocity)
@@ -733,19 +806,31 @@ bool ArgumentParser::ParseIniFile(const string& inifile)
 
                if ( xMethod_D->FirstChildElement("output_voronoi_cells"))
                {
-                    if ( string(xMethod_D->FirstChildElement("output_voronoi_cells")->Attribute("enabled"))=="true")
-                    {
-                         _isOutputGraph=true;
-                    	 Log->Write("INFO: \tData of voronoi diagram is asked to output" );
-                    	 if(string(xMethod_D->FirstChildElement("output_voronoi_cells")->Attribute("plot_graphs"))=="true")
-                    	 {
-                              _isPlotGraph=true;
-                              if(_isPlotGraph)
+                    auto enabled = xMethod_D->FirstChildElement("output_voronoi_cells")->Attribute("enabled");
+                    if(enabled)
+                         if ( string(enabled)=="true")
+                         {
+                              _isOutputGraph=true;
+                              Log->Write("INFO: \tData of voronoi diagram is asked to output" );
+                              auto plot_graphs = xMethod_D->FirstChildElement("output_voronoi_cells")->Attribute("plot_graphs");
+                              if(plot_graphs)
                               {
-                                   Log->Write("INFO: \tGraph of voronoi diagram will be plotted" );
-                              }
-                    	 }
-                    }
+                                   if (string(plot_graphs)=="true")
+                                   {
+                                        _isPlotGraph = true;
+                                        Log->Write("INFO: \tGraph of voronoi diagram will be plotted");
+                                   }
+                                   auto plot_index = xMethod_D->FirstChildElement("output_voronoi_cells")->Attribute(
+                                           "plot_index");
+                                   if (plot_index)
+                                        if (string(plot_index)=="true")
+                                        {
+                                             _isPlotIndex = true;
+                                             Log->Write(
+                                                     "INFO: \tVoronoi diagram will be plotted with index of pedestrians");
+                                        } // plot_index
+                              } // plot_graphs
+                         }// enabled
                }
 
                if ( xMethod_D->FirstChildElement("steadyState"))
@@ -873,6 +958,11 @@ bool ArgumentParser::GetIsPlotGraph() const
      return _isPlotGraph;
 }
 
+bool ArgumentParser::GetIsPlotIndex() const
+{
+     return _isPlotIndex;
+}
+
 vector<bool> ArgumentParser::GetIsPlotTimeSeriesA() const
 {
      return _isPlotTimeSeriesA;
@@ -966,3 +1056,12 @@ MeasurementArea* ArgumentParser::GetMeasurementArea(int id)
      return _measurementAreas[id];
 
 }
+
+void ArgumentParser::SetErrorLogFile(std::string errorLogFile)
+{
+     _errorLogFile = errorLogFile;
+};
+
+void ArgumentParser::SetLog(int log) {
+     _log = log;
+};
diff --git a/general/ArgumentParser.h b/general/ArgumentParser.h
index d85f1c3b822fb1edefdf18a621b6c4ace08025c2..3a5df0f53dbafe1be6052a1dfc1b99916a9d46d1 100644
--- a/general/ArgumentParser.h
+++ b/general/ArgumentParser.h
@@ -2,7 +2,7 @@
  * \file        ArgumentParser.cpp
  * \date        Oct 10, 2014
  * \version     v0.7
- * \copyright   <2009-2015> Forschungszentrum J��lich GmbH. All rights reserved.
+ * \copyright   <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
  *
  * \section License
  * This file is part of JuPedSim.
@@ -68,6 +68,7 @@ private:
      int _circleEdges;
      bool _isOutputGraph;
      bool _isPlotGraph;
+     bool _isPlotIndex;
      /*bool _isPlotTimeSeriesA;
      bool _isPlotTimeSeriesC;
      bool _isPlotTimeSeriesD;*/
@@ -136,6 +137,7 @@ public:
      int GetCircleEdges() const;
      bool GetIsOutputGraph() const;
      bool GetIsPlotGraph() const;
+     bool GetIsPlotIndex() const;
      std::vector<bool> GetIsPlotTimeSeriesA() const;
      std::vector<bool> GetIsPlotTimeSeriesC() const;
      std::vector<bool> GetIsPlotTimeSeriesD() const;
@@ -149,7 +151,8 @@ public:
      float GetGridSizeY() const;
      int GetLog() const;
      bool ParseArgs(int argc, char **argv);
-
+     void SetErrorLogFile(std::string errorLogFile);
+     void SetLog(int log);
      MeasurementArea* GetMeasurementArea(int id);
 
      /**
diff --git a/getRSS.c b/getRSS.c
deleted file mode 100644
index 574d373f97caf5b2228e792106f71a17158216b5..0000000000000000000000000000000000000000
--- a/getRSS.c
+++ /dev/null
@@ -1,122 +0,0 @@
-/*
- * Author:  David Robert Nadeau
- * Site:    http://NadeauSoftware.com/
- * License: Creative Commons Attribution 3.0 Unported License
- *          http://creativecommons.org/licenses/by/3.0/deed.en_US
- */
-
-#if defined(_WIN32)
-#include <windows.h>
-#include <psapi.h>
-
-#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
-#include <unistd.h>
-#include <sys/resource.h>
-
-#if defined(__APPLE__) && defined(__MACH__)
-#include <mach/mach.h>
-
-#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
-#include <fcntl.h>
-#include <procfs.h>
-
-#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
-#include <stdio.h>
-
-#endif
-
-#else
-#error "Cannot define getPeakRSS( ) or getCurrentRSS( ) for an unknown OS."
-#endif
-
-
-
-
-
-/**
- * Returns the peak (maximum so far) resident set size (physical
- * memory use) measured in bytes, or zero if the value cannot be
- * determined on this OS.
- */
-size_t getPeakRSS( )
-{
-#if defined(_WIN32)
-	/* Windows -------------------------------------------------- */
-	PROCESS_MEMORY_COUNTERS info;
-	GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
-	return (size_t)info.PeakWorkingSetSize;
-
-#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
-	/* AIX and Solaris ------------------------------------------ */
-	struct psinfo psinfo;
-	int fd = -1;
-	if ( (fd = open( "/proc/self/psinfo", O_RDONLY )) == -1 )
-		return (size_t)0L;		/* Can't open? */
-	if ( read( fd, &psinfo, sizeof(psinfo) ) != sizeof(psinfo) )
-	{
-		close( fd );
-		return (size_t)0L;		/* Can't read? */
-	}
-	close( fd );
-	return (size_t)(psinfo.pr_rssize * 1024L);
-
-#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
-	/* BSD, Linux, and OSX -------------------------------------- */
-	struct rusage rusage;
-	getrusage( RUSAGE_SELF, &rusage );
-#if defined(__APPLE__) && defined(__MACH__)
-	return (size_t)rusage.ru_maxrss;
-#else
-	return (size_t)(rusage.ru_maxrss * 1024L);
-#endif
-
-#else
-	/* Unknown OS ----------------------------------------------- */
-	return (size_t)0L;			/* Unsupported. */
-#endif
-}
-
-
-
-
-
-/**
- * Returns the current resident set size (physical memory use) measured
- * in bytes, or zero if the value cannot be determined on this OS.
- */
-size_t getCurrentRSS( )
-{
-#if defined(_WIN32)
-	/* Windows -------------------------------------------------- */
-	PROCESS_MEMORY_COUNTERS info;
-	GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
-	return (size_t)info.WorkingSetSize;
-
-#elif defined(__APPLE__) && defined(__MACH__)
-	/* OSX ------------------------------------------------------ */
-	struct mach_task_basic_info info;
-	mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
-	if ( task_info( mach_task_self( ), MACH_TASK_BASIC_INFO,
-		(task_info_t)&info, &infoCount ) != KERN_SUCCESS )
-		return (size_t)0L;		/* Can't access? */
-	return (size_t)info.resident_size;
-
-#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
-	/* Linux ---------------------------------------------------- */
-	long rss = 0L;
-	FILE* fp = NULL;
-	if ( (fp = fopen( "/proc/self/statm", "r" )) == NULL )
-		return (size_t)0L;		/* Can't open? */
-	if ( fscanf( fp, "%*s%ld", &rss ) != 1 )
-	{
-		fclose( fp );
-		return (size_t)0L;		/* Can't read? */
-	}
-	fclose( fp );
-	return (size_t)rss * (size_t)sysconf( _SC_PAGESIZE);
-
-#else
-	/* AIX, BSD, Solaris, and Unknown OS ------------------------ */
-	return (size_t)0L;			/* Unsupported. */
-#endif
-}
diff --git a/main.cpp b/main.cpp
index d02f63090e6ecd12ed343b5b1a6c33089abdf3c7..06a0855ce8f991fcb99abf1e6d8dcfc910016b1b 100644
--- a/main.cpp
+++ b/main.cpp
@@ -25,7 +25,7 @@
  *
  * Some useful links:
  *
- * 	1: <a href="http://jupedsim.org">jupedsim.org</a> <br>
+ *      1: <a href="http://jupedsim.org">jupedsim.org</a> <br>
  *
  **/
 
@@ -37,10 +37,9 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-     Log = new STDIOHandler();
-
+      Log = new STDIOHandler();
      // Parsing the arguments
-     ArgumentParser* args = new ArgumentParser();
+      ArgumentParser* args = new ArgumentParser();
 
 
      if(args->ParseArgs(argc, argv))
@@ -48,8 +47,6 @@ int main(int argc, char **argv)
           // get the number of file to analyse
           const vector<string>& files = args->GetTrajectoriesFiles();
           const string& path = args->GetTrajectoriesLocation();
-          //path="";
-
           // create and initialize the analysis engine
           for (unsigned int i = 0; i < files.size(); i++)
           {
diff --git a/methods/Method_A.cpp b/methods/Method_A.cpp
index f4c40cf18562050f13472dcc7c7dfce646b65121..575b8ac89321e23251fe36870f2da68636b49b2f 100644
--- a/methods/Method_A.cpp
+++ b/methods/Method_A.cpp
@@ -78,8 +78,8 @@ bool Method_A::Process (const PedData& peddata,const string& scriptsLocation, co
      bool PedInGeometry =false;
      for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();ite++)
      {
-    	  int frameNr = ite->first;
-    	  int frid =  frameNr + peddata.GetMinFrame();
+          int frameNr = ite->first;
+          int frid =  frameNr + peddata.GetMinFrame();
           if(!(frid%100))
           {
                Log->Write("frame ID = %d",frid);
@@ -91,21 +91,21 @@ bool Method_A::Process (const PedData& peddata,const string& scriptsLocation, co
           const vector<double> VInFrame = peddata.GetVInFrame(frameNr, ids, zPos_measureArea);
           if(IdInFrame.size()>0)
           {
-        	  GetAccumFlowVelocity(frameNr, IdInFrame, VInFrame);
-			  char tmp[30];
-			  sprintf(tmp, "%.2f\t%d\n", frid/_fps, _classicFlow);
-			  outputRhoV.append(tmp);
-			  PedInGeometry=true;
+               GetAccumFlowVelocity(frameNr, IdInFrame, VInFrame);
+               char tmp[30];
+               sprintf(tmp, "%.2f\t%d\n", frid/_fps, _classicFlow);
+               outputRhoV.append(tmp);
+               PedInGeometry=true;
           }
      }
      if(PedInGeometry)
      {
-		 FlowRate_Velocity(peddata.GetFps(), _accumPedsPassLine,_accumVPassLine);
-		 WriteFile_N_t(outputRhoV);
+          FlowRate_Velocity(peddata.GetFps(), _accumPedsPassLine,_accumVPassLine);
+          WriteFile_N_t(outputRhoV);
      }
      else
      {
-    	 Log->Write("Warning: No pedestrian exists on the plane of the selected Measurement area!!");
+          Log->Write("Warning: No pedestrian exists on the plane of the selected Measurement area!!");
      }
      delete []_passLine;
      return true;
@@ -123,9 +123,9 @@ void Method_A::WriteFile_N_t(string data)
           string file_N_t ="Flow_NT_"+_trajName+"_id_"+_measureAreaId+".dat";
           if(_plotTimeSeries)
           {
-			  string parameters_N_t="python3 \""+_scriptsLocation+"/_Plot_N_t.py\" -p \""+ METHOD_A_LOCATION + "\" -n "+file_N_t;
-			  int res = system(parameters_N_t.c_str());
-			  Log->Write("INFO:\tPlotting N-t diagram! Status: %d", res);
+               string parameters_N_t="python \""+_scriptsLocation+"/_Plot_N_t.py\" -p \""+ METHOD_A_LOCATION + "\" -n "+file_N_t;
+               int res = system(parameters_N_t.c_str());
+               Log->Write("INFO:\tPlotting N-t diagram! Status: %d", res);
           }
      }
      else
@@ -136,18 +136,18 @@ void Method_A::WriteFile_N_t(string data)
 
 void Method_A::GetAccumFlowVelocity(int frame, const vector<int>& ids, const vector<double>& VInFrame)
 {
-	for(unsigned int i=0; i<ids.size();i++)
+     for(unsigned int i=0; i<ids.size();i++)
      {
-		  int id = ids[i];
+          int id = ids[i];
           bool IspassLine=false;
           if(frame >_firstFrame[id-_minId]&&!_passLine[id-_minId])
           {
-        	  IspassLine = IsPassLine(_areaForMethod_A->_lineStartX,
-                                          _areaForMethod_A->_lineStartY,
-                                          _areaForMethod_A->_lineEndX,
-                                          _areaForMethod_A->_lineEndY, _xCor(id-_minId,frame - 1),
-                                          _yCor(id-_minId,frame - 1), _xCor(id-_minId,frame), 
-                                          _yCor(id-_minId,frame));
+               IspassLine = IsPassLine(_areaForMethod_A->_lineStartX,
+                                       _areaForMethod_A->_lineStartY,
+                                       _areaForMethod_A->_lineEndX,
+                                       _areaForMethod_A->_lineEndY, _xCor(id-_minId,frame - 1),
+                                       _yCor(id-_minId,frame - 1), _xCor(id-_minId,frame),
+                                       _yCor(id-_minId,frame));
           }
           if(IspassLine==true)
           {
diff --git a/methods/Method_D.cpp b/methods/Method_D.cpp
index 2139806828567e3d365c2395036242b721e40036..39ccc0ba305abcd9a808562c30f3dfcf44c944b9 100644
--- a/methods/Method_D.cpp
+++ b/methods/Method_D.cpp
@@ -2,7 +2,7 @@
  * \file        Method_D.cpp
  * \date        Oct 10, 2014
  * \version     v0.7
- * \copyright   <2009-2015> Forschungszentrum J��lich GmbH. All rights reserved.
+ * \copyright   <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
  *
  * \section License
  * This file is part of JuPedSim.
@@ -28,7 +28,9 @@
 
 #include "Method_D.h"
 #include <cmath>
-
+#include<map>
+#include<iostream>
+#include<vector>
 //using std::string;
 //using std::vector;
 //using std::ofstream;
@@ -64,35 +66,45 @@ Method_D::~Method_D()
 {
 
 }
-
+// auto outOfRange = [&](int number, int start, int end) -> bool
+// {
+//      return (number < start || number > end);
+// };
 bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocation, const double& zPos_measureArea)
 {
-	 _scriptsLocation = scriptsLocation;
+     _scriptsLocation = scriptsLocation;
      _peds_t = peddata.GetPedsFrame();
      _trajName = peddata.GetTrajName();
      _projectRootDir = peddata.GetProjectRootDir();
      _measureAreaId = boost::lexical_cast<string>(_areaForMethod_D->_id);
      _fps =peddata.GetFps();
-     Log->Write("INFO:\tFrame rate fps: <%.2f>", _fps);
-
+     int mycounter = 0;
      int minFrame = peddata.GetMinFrame();
+     Log->Write("INFO:\tMethod D: frame rate fps: <%.2f>, start: <%d>, stop: <%d> (minFrame = %d)", _fps, _startFrame, _stopFrame, minFrame);
      if(_startFrame!=_stopFrame)
      {
-         	 if(_startFrame==-1)
-         	 {
-         		_startFrame = minFrame;
-         	 }
-         	 if(_stopFrame==-1)
-         	 {
-         		_stopFrame = peddata.GetNumFrames()+minFrame;
-         	 }
-    	 	 for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();ite++)
-         	 {
-         		 if((ite->first + minFrame)<_startFrame || (ite->first + minFrame) >_stopFrame)
-         		 {
-         			_peds_t.erase(ite);
-         		 }
-         	 }
+          if(_startFrame==-1)
+          {
+               _startFrame = minFrame;
+          }
+          if(_stopFrame==-1)
+          {
+               _stopFrame = peddata.GetNumFrames()+minFrame;
+          }
+          for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();)
+          {
+               if((ite->first + minFrame)<_startFrame || (ite->first + minFrame) >_stopFrame)
+               {
+                    mycounter++;
+                    ite = _peds_t.erase(ite);
+               }
+               else
+               {
+                    ++ite;
+               }
+
+
+          }
      }
      OpenFileMethodD();
      if(_calcIndividualFD)
@@ -101,16 +113,17 @@ bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocatio
      }
      Log->Write("------------------------Analyzing with Method D-----------------------------");
      //for(int frameNr = 0; frameNr < peddata.GetNumFrames(); frameNr++ )
-     for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();ite++)
+     //for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();ite++)
+     for(auto ite: _peds_t)
      {
-    	  int frameNr = ite->first;
-    	  int frid =  frameNr + minFrame;
+          int frameNr = ite.first;
+          int frid =  frameNr + minFrame;
           //padd the frameid with 0
           std::ostringstream ss;
           ss << std::setw(5) << std::setfill('0') << frid;
           const std::string str_frid = ss.str();
 
-          if(!(frid%100))
+          if(!(frid%10))
           {
                Log->Write("frame ID = %d",frid);
           }
@@ -123,67 +136,81 @@ bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocatio
           //------------------------------Remove peds outside geometry------------------------------------------
           for( int i=0;i<(int)IdInFrame.size();i++)
           {
-        	  if(false==within(point_2d(round(XInFrame[i]), round(YInFrame[i])), _geoPoly))
-        	  {
-        		  //Log->Write("Warning:\tPedestrian at <x=%.4f, y=%.4f> is not in geometry and not considered in analysis!", XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
-        		  IdInFrame.erase(IdInFrame.begin() + i);
-        		  XInFrame.erase(XInFrame.begin() + i);
-        		  YInFrame.erase(YInFrame.begin() + i);
-        		  VInFrame.erase(VInFrame.begin() + i);
-        		  i--;
-        	  }
+               if(false==within(point_2d(round(XInFrame[i]), round(YInFrame[i])), _geoPoly))
+               {
+                    Log->Write("Warning:\tPedestrian at <x=%.4f, y=%.4f> is not in geometry and not considered in analysis!", XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
+                    IdInFrame.erase(IdInFrame.begin() + i);
+                    XInFrame.erase(XInFrame.begin() + i);
+                    YInFrame.erase(YInFrame.begin() + i);
+                    VInFrame.erase(VInFrame.begin() + i);
+                    i--;
+               }
           }
           int NumPeds = IdInFrame.size();
           //---------------------------------------------------------------------------------------------------------------
+          // std::cout << "NumPeds: " << NumPeds << "  Frame: "<< frid << "\n";
           if(NumPeds>3)
           {
-        	  if(_isOneDimensional)
-        	  {
-                  CalcVoronoiResults1D(XInFrame, VInFrame, IdInFrame, _areaForMethod_D->_poly,str_frid);
-        	  }
-        	  else
-        	  {
-        		  if(IsPointsOnOneLine(XInFrame, YInFrame))
-        		  {
-        			  if(fabs(XInFrame[1]-XInFrame[0])<dmin)
-        			  {
-        				  XInFrame[1]+= offset;
-        			  }
-        			  else
-        			  {
-        				  YInFrame[1]+= offset;
-        			  }
-        		  }
-        		   vector<polygon_2d> polygons = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
-				   if(!polygons.empty())
-				   {
-					   OutputVoronoiResults(polygons, str_frid, VInFrame);
-					   if(_calcIndividualFD)
-					   {
-							if(false==_isOneDimensional)
-							{
-								 GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_D->_poly, str_frid);
-							}
-					   }
-					   if(_getProfile)
-					   { //	field analysis
-							GetProfiles(str_frid, polygons, VInFrame);
-					   }
-					   if(_outputVoronoiCellData)
-					   { // output the Voronoi polygons of a frame
-							OutputVoroGraph(str_frid, polygons, NumPeds, XInFrame, YInFrame,VInFrame);
-					   }
-				   }
-				   else
-				   {
-					   Log->Write("Warning:\tVoronoi Diagrams are not obtained!\n");
-					   return false;
-				   }
-        	  }
+               if(_isOneDimensional)
+               {
+                    CalcVoronoiResults1D(XInFrame, VInFrame, IdInFrame, _areaForMethod_D->_poly,str_frid);
+               }
+               else
+               {
+                    if(IsPointsOnOneLine(XInFrame, YInFrame))
+                    {
+                         if(fabs(XInFrame[1]-XInFrame[0])<dmin)
+                         {
+                              XInFrame[1]+= offset;
+                         }
+                         else
+                         {
+                              YInFrame[1]+= offset;
+                         }
+                    }
+                    std::vector<std::pair<polygon_2d, int> > polygons_id = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
+                    // std::cout << ">> polygons_id " << polygons_id.size() << "\n";
+                    vector<polygon_2d> polygons;
+                    for (auto p: polygons_id)
+                         polygons.push_back(p.first);
+
+                    if(!polygons.empty())
+                    {
+                         OutputVoronoiResults(polygons, str_frid, VInFrame); // TODO polygons_id
+                         if(_calcIndividualFD)
+                         {
+                              if(false==_isOneDimensional)
+                              {
+                                   GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_D->_poly, str_frid); // TODO polygons_id
+                              }
+                         }
+                         if(_getProfile)
+                         { //	field analysis
+                              GetProfiles(str_frid, polygons, VInFrame); // TODO polygons_id
+                         }
+                         if(_outputVoronoiCellData)
+                         { // output the Voronoi polygons of a frame
+                              OutputVoroGraph(str_frid, polygons_id, NumPeds, XInFrame, YInFrame,VInFrame);
+                         }
+                    }
+                    else
+                    {
+
+                         for(int i=0;i<(int)IdInFrame.size();i++)
+                         {
+                              std::cout << XInFrame[i]*CMtoM << "   " << YInFrame[i]*CMtoM <<  "   "  << IdInFrame[i] << "\n";
+                         }
+                         Log->Write("WARNING: \tVoronoi Diagrams are not obtained!. Frame: %d\n", frid);
+                              return true; // todo: return false and
+                                           // check when does
+                                           // intersetion with circle fail
+                    }
+               }
           }
           else
           {
-               //Log->Write("INFO: \tThe number of the pedestrians is less than 2 !!");
+               Log->Write("DEBUG: \tThe number of the pedestrians is less than 2. Frame = %d\n", frid);
+               return true;
           }
      }
      fclose(_fVoronoiRhoV);
@@ -199,19 +226,19 @@ bool Method_D::OpenFileMethodD()
      string results_V=  _projectRootDir+VORO_LOCATION+"rho_v_Voronoi_"+_trajName+"_id_"+_measureAreaId+".dat";
      if((_fVoronoiRhoV=Analysis::CreateFile(results_V))==NULL)
      {
-          Log->Write("cannot open the file to write Voronoi density and velocity\n");
+          Log->Write("ERROR: \tcannot open the file to write Voronoi density and velocity\n");
           return false;
      }
      else
      {
-         if(_isOneDimensional)
-         {
-        	 fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-1))\t	Voronoi velocity(m/s)\n",_fps);
-         }
-         else
-         {
-        	 fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-2))\t	Voronoi velocity(m/s)\n",_fps);
-         }
+          if(_isOneDimensional)
+          {
+               fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-1))\t	Voronoi velocity(m/s)\n",_fps);
+          }
+          else
+          {
+               fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-2))\t	Voronoi velocity(m/s)\n",_fps);
+          }
           return true;
      }
 }
@@ -221,40 +248,48 @@ bool Method_D::OpenFileIndividualFD()
      string Individualfundment=_projectRootDir+"./Output/Fundamental_Diagram/Individual_FD/IndividualFD"+_trajName+"_id_"+_measureAreaId+".dat";
      if((_fIndividualFD=Analysis::CreateFile(Individualfundment))==NULL)
      {
-          Log->Write("cannot open the file individual\n");
+          Log->Write("ERROR:\tcannot open the file individual\n");
           return false;
      }
      else
      {
-    	 if(_isOneDimensional)
-    	 {
-    		 fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PedId	\t	Individual density(m^(-1)) \t 	Individual velocity(m/s)	\t	Headway(m)\n",_fps);
-    	 }
-    	 else
-    	 {
-    		 fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PedId	\t	Individual density(m^(-2)) \t 	Individual velocity(m/s)\n",_fps);
-    	 }
+          if(_isOneDimensional)
+          {
+               fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PedId	\t	Individual density(m^(-1)) \t   Individual velocity(m/s)	\t	Headway(m)\n",_fps);
+          }
+          else
+          {
+               fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PedId	\t	Individual density(m^(-2)) \t   Individual velocity(m/s)\n",_fps);
+          }
           return true;
      }
 }
 
-vector<polygon_2d> Method_D::GetPolygons(vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& VInFrame, vector<int>& IdInFrame)
+std::vector<std::pair<polygon_2d, int> > Method_D::GetPolygons(vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& VInFrame, vector<int>& IdInFrame)
 {
      VoronoiDiagram vd;
      //int NrInFrm = ids.size();
      double boundpoint =10*max(max(fabs(_geoMinX),fabs(_geoMinY)), max(fabs(_geoMaxX), fabs(_geoMaxY)));
-     vector<polygon_2d>  polygons = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
+     std::vector<std::pair<polygon_2d, int> > polygons_id;
+     polygons_id = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
+     // std:: cout << " GetPolygons " << polygons_id.size() << "\n";
+
+     polygon_2d poly ;
      if(_cutByCircle)
      {
-          polygons = vd.cutPolygonsWithCircle(polygons, XInFrame, YInFrame, _cutRadius,_circleEdges);
+          polygons_id = vd.cutPolygonsWithCircle(polygons_id, XInFrame, YInFrame, _cutRadius,_circleEdges);
      }
-     polygons = vd.cutPolygonsWithGeometry(polygons, _geoPoly, XInFrame, YInFrame);
-
-     for(auto && p:polygons)
+     // std:: cout << " GetPolygons cirlces " << polygons_id.size() << "\n";
+     polygons_id = vd.cutPolygonsWithGeometry(polygons_id, _geoPoly, XInFrame, YInFrame);
+     // std:: cout << " GetPolygons geometry " << polygons_id.size() << "\n";
+     for(auto && p:polygons_id)
      {
-          ReducePrecision(p);
+          poly = p.first;
+          ReducePrecision(poly);
+          // TODO update polygon_id?
      }
-     return polygons;
+     // std:: cout << " GetPolygons leave " << polygons_id.size() << "\n";
+     return polygons_id;
 }
 /**
  * Output the Voronoi density and velocity in the corresponding file
@@ -283,7 +318,7 @@ double Method_D::GetVoronoiDensity(const vector<polygon_2d>& polygon, const poly
                     std::cout<<"measure area: \t"<<std::setprecision(16)<<dsv(measureArea)<<"\n";
                     std::cout<<"Original polygon:\t"<<std::setprecision(16)<<dsv(polygon_iterator)<<"\n";
                     std::cout<<"intersected polygon: \t"<<std::setprecision(16)<<dsv(v[0])<<"\n";
-                    std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<<"\n";
+                    std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<<  "  (diff=" << (area(v[0]) - area(polygon_iterator)) << ")" << "\n";
                     //exit(EXIT_FAILURE);
                }
           }
@@ -333,7 +368,7 @@ double Method_D::GetVoronoiVelocity(const vector<polygon_2d>& polygon, const vec
                meanV+=(Velocity[temp]*area(v[0])/area(measureArea));
                if((area(v[0]) - area(polygon_iterator))>J_EPS)
                {
-                    std::cout<<"this is a wrong result in calculating velocity\t"<<area(v[0])<<'\t'<<area(polygon_iterator)<<std::endl;
+                    std::cout<<"this is a wrong result in calculating velocity\t"<<area(v[0])<<'\t'<<area(polygon_iterator)<< "  (diff=" << area(v[0]) - area(polygon_iterator) << ")"<< std::endl;
                }
           }
           temp++;
@@ -366,9 +401,9 @@ void Method_D::GetProfiles(const string& frameId, const vector<polygon_2d>& poly
                polygon_2d measurezoneXY;
                {
                     const double coor[][2] = {
-                              {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
-                              {_geoMinX+colum_j*_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
-                              {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y} // closing point is opening point
+                         {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
+                         {_geoMinX+colum_j*_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
+                         {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y} // closing point is opening point
                     };
                     assign_points(measurezoneXY, coor);
                }
@@ -385,10 +420,11 @@ void Method_D::GetProfiles(const string& frameId, const vector<polygon_2d>& poly
      fclose(Prf_density);
 }
 
-void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polygons, int numPedsInFrame, vector<double>& XInFrame, vector<double>& YInFrame,const vector<double>& VInFrame)
+void Method_D::OutputVoroGraph(const string & frameId,  std::vector<std::pair<polygon_2d, int> >& polygons_id, int numPedsInFrame, vector<double>& XInFrame, vector<double>& YInFrame,const vector<double>& VInFrame)
 {
      //string voronoiLocation=_projectRootDir+"./Output/Fundamental_Diagram/Classical_Voronoi/VoronoiCell/id_"+_measureAreaId;
      string voronoiLocation=_projectRootDir+VORO_LOCATION+"VoronoiCell/";
+     polygon_2d poly;
 
 
 #if defined(_WIN32)
@@ -397,28 +433,29 @@ void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polyg
      mkdir(voronoiLocation.c_str(), 0777);
 #endif
 
-
      string polygon=voronoiLocation+"/polygon"+_trajName+"_id_"+_measureAreaId+"_"+frameId+".dat";
      ofstream polys (polygon.c_str());
      if(polys.is_open())
      {
           //for(vector<polygon_2d> polygon_iterator=polygons.begin(); polygon_iterator!=polygons.end(); polygon_iterator++)
-    	  for(auto && poly:polygons)
+          for(auto && p:polygons_id)
           {
-        	  for(auto&& point:poly.outer())
-        	  {
-        		  point.x(point.x()*CMtoM);
-        		  point.y(point.y()*CMtoM);
-        	  }
-        	  for(auto&& innerpoly:poly.inners())
-			  {
-        		  for(auto&& point:innerpoly)
-				  {
-					  point.x(point.x()*CMtoM);
-					  point.y(point.y()*CMtoM);
-				  }
-			  }
-        	  polys << dsv(poly) << endl;
+               poly = p.first;
+               for(auto&& point:poly.outer())
+               {
+                    point.x(point.x()*CMtoM);
+                    point.y(point.y()*CMtoM);
+               }
+               for(auto&& innerpoly:poly.inners())
+               {
+                    for(auto&& point:innerpoly)
+                    {
+                         point.x(point.x()*CMtoM);
+                         point.y(point.y()*CMtoM);
+                    }
+               }
+               polys << p.second << " | " << dsv(poly) << endl;
+               //polys  <<dsv(poly)<< endl;
           }
      }
      else
@@ -443,26 +480,30 @@ void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polyg
      }
 
      /*string point=voronoiLocation+"/points"+_trajName+"_id_"+_measureAreaId+"_"+frameId+".dat";
-     ofstream points (point.c_str());
-     if( points.is_open())
-     {
-          for(int pts=0; pts<numPedsInFrame; pts++)
-          {
-               points << XInFrame[pts]*CMtoM << "\t" << YInFrame[pts]*CMtoM << endl;
-          }
-     }
-     else
-     {
-          Log->Write("ERROR:\tcannot create the file <%s>",point.c_str());
-          exit(EXIT_FAILURE);
-     }*/
+       ofstream points (point.c_str());
+       if( points.is_open())
+       {
+       for(int pts=0; pts<numPedsInFrame; pts++)
+       {
+       points << XInFrame[pts]*CMtoM << "\t" << YInFrame[pts]*CMtoM << endl;
+       }
+       }
+       else
+       {
+       Log->Write("ERROR:\tcannot create the file <%s>",point.c_str());
+       exit(EXIT_FAILURE);
+       }*/
 
      if(_plotVoronoiCellData)
      {
-          string parameters_rho="python3 "+_scriptsLocation+"/_Plot_cell_rho.py -f \""+ voronoiLocation + "\" -n "+ _trajName+"_id_"+_measureAreaId+"_"+frameId+
+          string parameters_rho="python "+_scriptsLocation+"/_Plot_cell_rho.py -f \""+ voronoiLocation + "\" -n "+ _trajName+"_id_"+_measureAreaId+"_"+frameId+
                " -g "+_geometryFileName+" -p "+_trajectoryPath;
-          string parameters_v="python3 "+_scriptsLocation+"/_Plot_cell_v.py -f \""+ voronoiLocation + "\" -n "+ _trajName+"_id_"+_measureAreaId+"_"+frameId+
+          string parameters_v="python "+_scriptsLocation+"/_Plot_cell_v.py -f \""+ voronoiLocation + "\" -n "+ _trajName+"_id_"+_measureAreaId+"_"+frameId+
                " -g "+_geometryFileName+" -p "+_trajectoryPath;
+
+          if(_plotVoronoiIndex)
+               parameters_rho += " -i";
+
           Log->Write("INFO:\t%s",parameters_rho.c_str());
           Log->Write("INFO:\tPlotting Voronoi Cell at the frame <%s>",frameId.c_str());
           system(parameters_rho.c_str());
@@ -501,19 +542,19 @@ void Method_D::SetCalculateIndividualFD(bool individualFD)
 
 void Method_D::SetStartFrame(int startFrame)
 {
-	_startFrame=startFrame;
+     _startFrame=startFrame;
 }
 
 void Method_D::SetStopFrame(int stopFrame)
 {
-	_stopFrame=stopFrame;
+     _stopFrame=stopFrame;
 }
 
 void Method_D::Setcutbycircle(double radius,int edges)
 {
-	_cutByCircle=true;
-	_cutRadius = radius;
-	_circleEdges = edges;
+     _cutByCircle=true;
+     _cutRadius = radius;
+     _circleEdges = edges;
 }
 
 void Method_D::SetGeometryPolygon(polygon_2d geometryPolygon)
@@ -523,15 +564,15 @@ void Method_D::SetGeometryPolygon(polygon_2d geometryPolygon)
 
 void Method_D::SetGeometryBoundaries(double minX, double minY, double maxX, double maxY)
 {
-	_geoMinX = minX;
-	_geoMinY = minY;
-	_geoMaxX = maxX;
-	_geoMaxY = maxY;
+     _geoMinX = minX;
+     _geoMinY = minY;
+     _geoMaxX = maxX;
+     _geoMaxY = maxY;
 }
 
 void Method_D::SetGeometryFileName(const std::string& geometryFile)
 {
-	_geometryFileName=geometryFile;
+     _geometryFileName=geometryFile;
 }
 
 void Method_D::SetTrajectoriesLocation(const std::string& trajectoryPath)
@@ -559,6 +600,10 @@ void Method_D::SetPlotVoronoiGraph(bool plotVoronoiGraph)
 {
      _plotVoronoiCellData = plotVoronoiGraph;
 }
+void Method_D::SetPlotVoronoiIndex(bool plotVoronoiIndex)
+{
+     _plotVoronoiIndex = plotVoronoiIndex;
+}
 
 void Method_D::SetMeasurementArea (MeasurementArea_B* area)
 {
@@ -572,7 +617,7 @@ void Method_D::SetDimensional (bool dimension)
 
 void Method_D::ReducePrecision(polygon_2d& polygon)
 {
-	for(auto&& point:polygon.outer())
+     for(auto&& point:polygon.outer())
      {
           point.x(round(point.x() * 100000000000.0) / 100000000000.0);
           point.y(round(point.y() * 100000000000.0) / 100000000000.0);
@@ -581,133 +626,133 @@ void Method_D::ReducePrecision(polygon_2d& polygon)
 
 bool Method_D::IsPedInGeometry(int frames, int peds, double **Xcor, double **Ycor, int  *firstFrame, int *lastFrame)
 {
-	for(int i=0; i<peds; i++)
-		for(int j =0; j<frames; j++)
-		{
-			if (j>firstFrame[i] && j< lastFrame[i] && (false==within(point_2d(round(Xcor[i][j]), round(Ycor[i][j])), _geoPoly)))
-			{
-				Log->Write("Error:\tPedestrian at the position <x=%.4f, y=%.4f> is outside geometry. Please check the geometry or trajectory file!", Xcor[i][j]*CMtoM, Ycor[i][j]*CMtoM );
-				return false;
-			}
-		}
-	return true;
+     for(int i=0; i<peds; i++)
+          for(int j =0; j<frames; j++)
+          {
+               if (j>firstFrame[i] && j< lastFrame[i] && (false==within(point_2d(round(Xcor[i][j]), round(Ycor[i][j])), _geoPoly)))
+               {
+                    Log->Write("Error:\tPedestrian at the position <x=%.4f, y=%.4f> is outside geometry. Please check the geometry or trajectory file!", Xcor[i][j]*CMtoM, Ycor[i][j]*CMtoM );
+                    return false;
+               }
+          }
+     return true;
 }
 
 void Method_D::CalcVoronoiResults1D(vector<double>& XInFrame, vector<double>& VInFrame, vector<int>& IdInFrame, const polygon_2d & measureArea,const string& frid)
 {
-	vector<double> measurearea_x;
-	for(unsigned int i=0;i<measureArea.outer().size();i++)
-	{
-		measurearea_x.push_back(measureArea.outer()[i].x());
-	}
-	double left_boundary=*min_element(measurearea_x.begin(),measurearea_x.end());
-	double right_boundary=*max_element(measurearea_x.begin(),measurearea_x.end());
-
-	vector<double> voronoi_distance;
-	vector<double> Xtemp=XInFrame;
-	vector<double> dist;
-	vector<double> XLeftNeighbor;
-	vector<double> XLeftTemp;
-	vector<double> XRightNeighbor;
-	vector<double> XRightTemp;
-	sort(Xtemp.begin(),Xtemp.end());
-	dist.push_back(Xtemp[1]-Xtemp[0]);
-	XLeftTemp.push_back(2*Xtemp[0]-Xtemp[1]);
-	XRightTemp.push_back(Xtemp[1]);
-	for(unsigned int i=1;i<Xtemp.size()-1;i++)
-	{
-		dist.push_back((Xtemp[i+1]-Xtemp[i-1])/2.0);
-		XLeftTemp.push_back(Xtemp[i-1]);
-		XRightTemp.push_back(Xtemp[i+1]);
-	}
-	dist.push_back(Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
-	XLeftTemp.push_back(Xtemp[Xtemp.size()-2]);
-	XRightTemp.push_back(2*Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
-	for(unsigned int i=0;i<XInFrame.size();i++)
-	{
-		for(unsigned int j=0;j<Xtemp.size();j++)
-		{
-			if(fabs(XInFrame[i]-Xtemp[j])<1.0e-5)
-			{
-				voronoi_distance.push_back(dist[j]);
-				XLeftNeighbor.push_back(XLeftTemp[j]);
-				XRightNeighbor.push_back(XRightTemp[j]);
-				break;
-			}
-		}
-	}
-
-	double VoronoiDensity=0;
-	double VoronoiVelocity=0;
-	for(unsigned int i=0; i<XInFrame.size(); i++)
-	{
-		double ratio=getOverlapRatio((XInFrame[i]+XLeftNeighbor[i])/2.0, (XRightNeighbor[i]+XInFrame[i])/2.0,left_boundary,right_boundary);
-		VoronoiDensity+=ratio;
-		VoronoiVelocity+=(VInFrame[i]*voronoi_distance[i]*ratio*CMtoM);
-		if(_calcIndividualFD)
-		{
-			double headway=(XRightNeighbor[i] - XInFrame[i])*CMtoM;
-			double individualDensity = 2.0/((XRightNeighbor[i] - XLeftNeighbor[i])*CMtoM);
-			fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%.3f\n",frid.c_str(), IdInFrame[i], individualDensity,VInFrame[i], headway);
-		}
-	}
-	VoronoiDensity/=((right_boundary-left_boundary)*CMtoM);
-	VoronoiVelocity/=((right_boundary-left_boundary)*CMtoM);
-	fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
+     vector<double> measurearea_x;
+     for(unsigned int i=0;i<measureArea.outer().size();i++)
+     {
+          measurearea_x.push_back(measureArea.outer()[i].x());
+     }
+     double left_boundary=*min_element(measurearea_x.begin(),measurearea_x.end());
+     double right_boundary=*max_element(measurearea_x.begin(),measurearea_x.end());
+
+     vector<double> voronoi_distance;
+     vector<double> Xtemp=XInFrame;
+     vector<double> dist;
+     vector<double> XLeftNeighbor;
+     vector<double> XLeftTemp;
+     vector<double> XRightNeighbor;
+     vector<double> XRightTemp;
+     sort(Xtemp.begin(),Xtemp.end());
+     dist.push_back(Xtemp[1]-Xtemp[0]);
+     XLeftTemp.push_back(2*Xtemp[0]-Xtemp[1]);
+     XRightTemp.push_back(Xtemp[1]);
+     for(unsigned int i=1;i<Xtemp.size()-1;i++)
+     {
+          dist.push_back((Xtemp[i+1]-Xtemp[i-1])/2.0);
+          XLeftTemp.push_back(Xtemp[i-1]);
+          XRightTemp.push_back(Xtemp[i+1]);
+     }
+     dist.push_back(Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
+     XLeftTemp.push_back(Xtemp[Xtemp.size()-2]);
+     XRightTemp.push_back(2*Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
+     for(unsigned int i=0;i<XInFrame.size();i++)
+     {
+          for(unsigned int j=0;j<Xtemp.size();j++)
+          {
+               if(fabs(XInFrame[i]-Xtemp[j])<1.0e-5)
+               {
+                    voronoi_distance.push_back(dist[j]);
+                    XLeftNeighbor.push_back(XLeftTemp[j]);
+                    XRightNeighbor.push_back(XRightTemp[j]);
+                    break;
+               }
+          }
+     }
+
+     double VoronoiDensity=0;
+     double VoronoiVelocity=0;
+     for(unsigned int i=0; i<XInFrame.size(); i++)
+     {
+          double ratio=getOverlapRatio((XInFrame[i]+XLeftNeighbor[i])/2.0, (XRightNeighbor[i]+XInFrame[i])/2.0,left_boundary,right_boundary);
+          VoronoiDensity+=ratio;
+          VoronoiVelocity+=(VInFrame[i]*voronoi_distance[i]*ratio*CMtoM);
+          if(_calcIndividualFD)
+          {
+               double headway=(XRightNeighbor[i] - XInFrame[i])*CMtoM;
+               double individualDensity = 2.0/((XRightNeighbor[i] - XLeftNeighbor[i])*CMtoM);
+               fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%.3f\n",frid.c_str(), IdInFrame[i], individualDensity,VInFrame[i], headway);
+          }
+     }
+     VoronoiDensity/=((right_boundary-left_boundary)*CMtoM);
+     VoronoiVelocity/=((right_boundary-left_boundary)*CMtoM);
+     fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
 
 }
 
 double Method_D::getOverlapRatio(const double& left, const double& right, const double& measurearea_left, const double& measurearea_right)
 {
-	double OverlapRatio=0;
-	double PersonalSpace=right-left;
-	if(left > measurearea_left && right < measurearea_right) //case1
-	{
-		OverlapRatio=1;
-	}
-	else if(right > measurearea_left && right < measurearea_right && left < measurearea_left)
-	{
-		OverlapRatio=(right-measurearea_left)/PersonalSpace;
-	}
-	else if(left < measurearea_left && right > measurearea_right)
-	{
-		OverlapRatio=(measurearea_right - measurearea_left)/PersonalSpace;
-	}
-	else if(left > measurearea_left && left < measurearea_right && right > measurearea_right)
-	{
-		OverlapRatio=(measurearea_right-left)/PersonalSpace;
-	}
-	return OverlapRatio;
+     double OverlapRatio=0;
+     double PersonalSpace=right-left;
+     if(left > measurearea_left && right < measurearea_right) //case1
+     {
+          OverlapRatio=1;
+     }
+     else if(right > measurearea_left && right < measurearea_right && left < measurearea_left)
+     {
+          OverlapRatio=(right-measurearea_left)/PersonalSpace;
+     }
+     else if(left < measurearea_left && right > measurearea_right)
+     {
+          OverlapRatio=(measurearea_right - measurearea_left)/PersonalSpace;
+     }
+     else if(left > measurearea_left && left < measurearea_right && right > measurearea_right)
+     {
+          OverlapRatio=(measurearea_right-left)/PersonalSpace;
+     }
+     return OverlapRatio;
 }
 
 bool Method_D::IsPointsOnOneLine(vector<double>& XInFrame, vector<double>& YInFrame)
 {
-	double deltaX=XInFrame[1] - XInFrame[0];
-	bool isOnOneLine=true;
-	if(fabs(deltaX)<dmin)
-	{
-		for(unsigned int i=2; i<XInFrame.size();i++)
-		{
-			if(fabs(XInFrame[i] - XInFrame[0])> dmin)
-			{
-				isOnOneLine=false;
-				break;
-			}
-		}
-	}
-	else
-	{
-		double slope=(YInFrame[1] - YInFrame[0])/deltaX;
-		double intercept=YInFrame[0] - slope*XInFrame[0];
-		for(unsigned int i=2; i<XInFrame.size();i++)
-		{
-			double dist=fabs(slope*XInFrame[i] - YInFrame[i] + intercept)/sqrt(slope*slope +1);
-			if(dist > dmin)
-			{
-				isOnOneLine=false;
-				break;
-			}
-		}
-	}
-	return isOnOneLine;
+     double deltaX=XInFrame[1] - XInFrame[0];
+     bool isOnOneLine=true;
+     if(fabs(deltaX)<dmin)
+     {
+          for(unsigned int i=2; i<XInFrame.size();i++)
+          {
+               if(fabs(XInFrame[i] - XInFrame[0])> dmin)
+               {
+                    isOnOneLine=false;
+                    break;
+               }
+          }
+     }
+     else
+     {
+          double slope=(YInFrame[1] - YInFrame[0])/deltaX;
+          double intercept=YInFrame[0] - slope*XInFrame[0];
+          for(unsigned int i=2; i<XInFrame.size();i++)
+          {
+               double dist=fabs(slope*XInFrame[i] - YInFrame[i] + intercept)/sqrt(slope*slope +1);
+               if(dist > dmin)
+               {
+                    isOnOneLine=false;
+                    break;
+               }
+          }
+     }
+     return isOnOneLine;
 }
diff --git a/methods/Method_D.h b/methods/Method_D.h
index a115b780b80c2ee01d780c3fae5ede9493ac0a2f..ffe2f55369cf3b1a38519b969492d9e33bc8381c 100644
--- a/methods/Method_D.h
+++ b/methods/Method_D.h
@@ -63,6 +63,7 @@ public:
      void SetCalculateProfiles(bool calcProfile);
      void SetOutputVoronoiCellData(bool outputCellData);
      void SetPlotVoronoiGraph(bool plotVoronoiGraph);
+     void SetPlotVoronoiIndex(bool plotVoronoiIndex);
      void SetMeasurementArea (MeasurementArea_B* area);
      void SetDimensional (bool dimension);
      void SetTrajectoriesLocation(const std::string& trajectoryPath);
@@ -81,6 +82,7 @@ private:
      bool _getProfile;
      bool _outputVoronoiCellData;
      bool _plotVoronoiCellData;
+     bool _plotVoronoiIndex;
      bool _isOneDimensional;
      bool _cutByCircle;       //Adjust whether cut each original voronoi cell by a circle
      double _cutRadius;
@@ -103,14 +105,14 @@ private:
      int _stopFrame;
 
 
-     std::vector<polygon_2d> GetPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame,
-               std::vector<double>& VInFrame, std::vector<int>& IdInFrame);
+     std::vector<std::pair<polygon_2d, int> >  GetPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame,
+                                                           std::vector<double>& VInFrame, std::vector<int>& IdInFrame);
      void OutputVoronoiResults(const std::vector<polygon_2d>&  polygons, const std::string& frid, const std::vector<double>& VInFrame);
      double GetVoronoiDensity(const std::vector<polygon_2d>& polygon, const polygon_2d & measureArea);
      double GetVoronoiDensity2(const std::vector<polygon_2d>& polygon, double* XInFrame, double* YInFrame, const polygon_2d& measureArea);
      double GetVoronoiVelocity(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const polygon_2d & measureArea);
      void GetProfiles(const std::string& frameId, const std::vector<polygon_2d>& polygons, const std::vector<double>& velocity);
-     void OutputVoroGraph(const std::string & frameId, std::vector<polygon_2d>& polygons, int numPedsInFrame,std::vector<double>& XInFrame,
+     void OutputVoroGraph(const std::string & frameId,  std::vector<std::pair<polygon_2d, int> >& polygons, int numPedsInFrame,std::vector<double>& XInFrame,
                std::vector<double>& YInFrame,const std::vector<double>& VInFrame);
      void GetIndividualFD(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const std::vector<int>& Id, const polygon_2d& measureArea, const std::string& frid);
      /**
diff --git a/methods/PedData.cpp b/methods/PedData.cpp
index 0e27ec17fed46eb093260a4772401dde748cb570..731b668efc43c1a6cd3be36784dbf3e3cbfa958a 100644
--- a/methods/PedData.cpp
+++ b/methods/PedData.cpp
@@ -35,8 +35,6 @@ using std::map;
 using std::vector;
 using std::ifstream;
 
-#include "../getRSS.c"
-
 PedData::PedData()
 {
 }
@@ -68,7 +66,7 @@ bool PedData::ReadData(const string& projectRootDir, const string& path, const s
                Log->Write("ERROR: \t could not parse the trajectories file <%s>",fullTrajectoriesPathName.c_str());
                return false;
           }
-           TiXmlElement* xRootNode = docGeo.RootElement();
+          TiXmlElement* xRootNode = docGeo.RootElement();
           result=InitializeVariables(xRootNode);	//initialize some global variables
      }
 
@@ -106,8 +104,10 @@ bool PedData::InitializeVariables(const string& filename)
           int pos_y=3;
           int pos_z=4;
           int pos_vd=5; //velocity direction
+          int fps_found = 0;
           while ( getline(fdata,line) )
           {
+               boost::algorithm::trim(line);
                //looking for the framerate which is suppposed to be at the second position
                if(line[0] == '#')
                {
@@ -118,11 +118,21 @@ bool PedData::InitializeVariables(const string& filename)
                          line.erase(0,1); // remove #
                          boost::split(strs, line , boost::is_any_of(":"),boost::token_compress_on);
                          if(strs.size()>1)
+                         {
                               _fps= atof(strs[1].c_str());
-
+                              if(_fps == 0.0) // in case not valid fps found
+                              {
+                                   Log->Write("ERROR:\t Could not convert fps <%s>", strs[1].c_str());
+                                   exit(1);
+                              }
+                         }
+                         else{
+                              Log->Write("ERROR:\tFrame rate fps not defined");
+                              exit(1);
+                         }
                          Log->Write("INFO:\tFrame rate fps: <%.2f>", _fps);
+                         fps_found = 1;
                     }
-
                     if(line.find("ID") != std::string::npos &&
                        line.find("FR") != std::string::npos &&
                        line.find("X") != std::string::npos &&
@@ -130,34 +140,40 @@ bool PedData::InitializeVariables(const string& filename)
                        line.find("Z") != std::string::npos)
                     {
                          // looking for this line
-                         // #ID	 FR  X Y Z
-                    	std::vector<std::string> strs1;
-                    	line.erase(0,1);
-                    	boost::split(strs1, line , boost::is_any_of("\t\r "),boost::token_compress_on);
-                    	vector<string>::iterator it_id;
-                    	it_id=find(strs1.begin(),strs1.end(),"ID");
-                    	pos_id = std::distance(strs1.begin(), it_id);
-                    	it_id=find(strs1.begin(),strs1.end(),"FR");
-                    	pos_fr = std::distance(strs1.begin(), it_id);
-                    	it_id=find(strs1.begin(),strs1.end(),"X");
-                    	pos_x = std::distance(strs1.begin(), it_id);
-                    	it_id=find(strs1.begin(),strs1.end(),"Y");
-                    	pos_y = std::distance(strs1.begin(), it_id);
-                    	it_id=find(strs1.begin(),strs1.end(),"Z");
-                        pos_z = std::distance(strs1.begin(), it_id);
-                        it_id=find(strs1.begin(),strs1.end(),"VD");
-                        pos_vd = std::distance(strs1.begin(), it_id);
+                         // #ID  FR  X Y Z
+                         //std::cout << "0. line <" << line << ">\n";
+                         std::vector<std::string> strs1;
+                         line.erase(0,1);
+                         //std::cout << "1. line <" << line << ">\n";
+                         boost::split(strs1, line , boost::is_any_of("\t\r "),boost::token_compress_on);
+                         // std::cout << "str size = " << strs1.size() << "\n";
+                         // for(auto s: strs1)
+                         //      std::cout << "s <" << s  << ">\n";
+
+                         vector<string>::iterator it_id;
+                         it_id=find(strs1.begin(),strs1.end(),"ID");
+                         pos_id = std::distance(strs1.begin(), it_id);
+                         it_id=find(strs1.begin(),strs1.end(),"FR");
+                         pos_fr = std::distance(strs1.begin(), it_id);
+                         it_id=find(strs1.begin(),strs1.end(),"X");
+                         pos_x = std::distance(strs1.begin(), it_id);
+                         it_id=find(strs1.begin(),strs1.end(),"Y");
+                         pos_y = std::distance(strs1.begin(), it_id);
+                         it_id=find(strs1.begin(),strs1.end(),"Z");
+                         pos_z = std::distance(strs1.begin(), it_id);
+                         it_id=find(strs1.begin(),strs1.end(),"VD");
+                         pos_vd = std::distance(strs1.begin(), it_id);
                     }
 
                }
                else if ( line[0] != '#' && !(line.empty()) )
                {
                     static int once = 1;
-                    if (lineNr % 100000 == 0)
-                         std::cout << "lineNr " << lineNr<< std::endl;
                     std::vector<std::string> strs;
-                    boost::algorithm::trim_right(line);
                     boost::split(strs, line , boost::is_any_of("\t "),boost::token_compress_on);
+                    if (lineNr % 100000 == 0)
+                         std::cout << "lineNr " << lineNr<< std::endl;
+
                     if(once) // && strs.size() < 5
                     {
                          once = 0;
@@ -166,14 +182,13 @@ bool PedData::InitializeVariables(const string& filename)
                          Log->Write("INFO: pos_x: %d", pos_x);
                          Log->Write("INFO: pos_y: %d", pos_y);
                          Log->Write("INFO: pos_z: %d", pos_z);
-                         // Log->Write("WARNING:\t assuming z=0 for all data");
+                         Log->Write("INFO: pos_vd: %d", pos_vd);
                     }
-
                     _IdsTXT.push_back(atoi(strs[pos_id].c_str()));
                     _FramesTXT.push_back(atoi(strs[pos_fr].c_str()));
                     xs.push_back(atof(strs[pos_x].c_str()));
                     ys.push_back(atof(strs[pos_y].c_str()));
-                    
+
                     if(strs.size() >= 5)
                          zs.push_back(atof(strs[pos_z].c_str()));
                     else
@@ -193,9 +208,14 @@ bool PedData::InitializeVariables(const string& filename)
                     }
                }
                lineNr++;
+          }// while
+          if(fps_found == 0)
+          {
+               Log->Write("ERROR:\tFrame rate fps not defined ");
+               exit(1);
           }
           Log->Write("INFO:\t Finished reading the data");
-          
+
      }
      fdata.close();
      Log->Write("INFO: Got %d lines", _IdsTXT.size());
@@ -210,7 +230,7 @@ bool PedData::InitializeVariables(const string& filename)
      //Total number of agents
      std::vector<int> unique_ids = _IdsTXT;
      // no need to
-     //sort. Assume that ids are ascendant 
+     //sort. Assume that ids are ascendant
      sort(unique_ids.begin(), unique_ids.end());
      std::vector<int>::iterator it;
      it = unique(unique_ids.begin(), unique_ids.end());
@@ -220,43 +240,40 @@ bool PedData::InitializeVariables(const string& filename)
      CreateGlobalVariables(_numPeds, _numFrames);
      Log->Write("INFO: Create Global Variables done");
 
-     
-     for(int i=_minID;i<_minID+_numPeds; i++)
-     {
-    	 int firstFrameIndex=INT_MAX;   //The first frame index of a pedestrian
-    	 int lastFrameIndex=-1;    //The last frame index of a pedestrian
-    	 int actual_totalframe=0;  //The total data points of a pedestrian in the trajectory
-    	 for (auto j = _IdsTXT.begin(); j != _IdsTXT.end(); ++j)
-    	     {
-    	         if (*j ==i)
-    	         {
-    	             int pos = std::distance(_IdsTXT.begin(), j);
-    	             if(pos<firstFrameIndex)
-    	             {
-    	            	 firstFrameIndex=pos;
-    	             }
-    	             if(pos>lastFrameIndex)
-    	             {
-    	            	 lastFrameIndex=pos;
-    	             }
-    	             actual_totalframe++;
-    	         }
-    	     }
-         
-    	 if(lastFrameIndex==0)
-    	 {
-    		 Log->Write("Warning:\tThere is no trajectory for ped with ID <%d>!",i);
-    	 }
-    	 _firstFrame[i-_minID] = _FramesTXT[firstFrameIndex] - _minFrame;
-    	 _lastFrame[i-_minID] = _FramesTXT[lastFrameIndex] - _minFrame;
-
-         int expect_totalframe=_lastFrame[i-_minID]-_firstFrame[i-_minID]+1;
-         if(actual_totalframe != expect_totalframe)
-         {
-              Log->Write("Error:\tThe trajectory of ped with ID <%d> is not continuous. Please modify the trajectory file!",i);
-              Log->Write("Error:\t actual_totalfame = <%d>, expected_totalframe = <%d> ", actual_totalframe, expect_totalframe);
-              return false;
-         }
+
+     for(int i=_minID;i<_minID+_numPeds; i++){
+          int firstFrameIndex=INT_MAX;   //The first frame index of a pedestrian
+          int lastFrameIndex=-1;    //The last frame index of a pedestrian
+          int actual_totalframe=0;  //The total data points of a pedestrian in the trajectory
+          for (auto j = _IdsTXT.begin(); j != _IdsTXT.end(); ++j){
+               if (*j ==i){
+                    int pos = std::distance(_IdsTXT.begin(), j);
+                    if(pos<firstFrameIndex)
+                    {
+                         firstFrameIndex=pos;
+                    }
+                    if(pos>lastFrameIndex)
+                    {
+                         lastFrameIndex=pos;
+                    }
+                    actual_totalframe++;
+               }
+          }
+          if(lastFrameIndex <=0 || lastFrameIndex == INT_MAX)
+          {
+               Log->Write("Warning:\tThere is no trajectory for ped with ID <%d>!",i);
+               continue;
+          }
+          _firstFrame[i-_minID] = _FramesTXT[firstFrameIndex] - _minFrame;
+          _lastFrame[i-_minID] = _FramesTXT[lastFrameIndex] - _minFrame;
+
+          int expect_totalframe=_lastFrame[i-_minID]-_firstFrame[i-_minID]+1;
+          if(actual_totalframe != expect_totalframe)
+          {
+               Log->Write("Error:\tThe trajectory of ped with ID <%d> is not continuous. Please modify the trajectory file!",i);
+               Log->Write("Error:\t actual_totalfame = <%d>, expected_totalframe = <%d> ", actual_totalframe, expect_totalframe);
+               return false;
+          }
      }
      Log->Write("convert x and y");
      for(unsigned int i = 0; i < _IdsTXT.size(); i++)
@@ -271,7 +288,7 @@ bool PedData::InitializeVariables(const string& filename)
           {
                Log->Write("Error:\t Id %d does not exist in file", _IdsTXT[i]);
                return false;
-          } 
+          }
           else
           {
                ID  = std::distance(unique_ids.begin(), it);
@@ -281,7 +298,7 @@ bool PedData::InitializeVariables(const string& filename)
           double x = xs[i]*M2CM;
           double y = ys[i]*M2CM;
           double z = zs[i]*M2CM;
-          
+
           _xCor(ID,frm) = x;
           _yCor(ID,frm) = y;
           _zCor(ID,frm) = z;
@@ -295,7 +312,7 @@ bool PedData::InitializeVariables(const string& filename)
           }
      }
      Log->Write("Save the data for each frame");
-     
+
      //save the data for each frame
      for (unsigned int i = 0; i < _FramesTXT.size(); i++ )
      {
@@ -308,12 +325,12 @@ bool PedData::InitializeVariables(const string& filename)
           {
                Log->Write("Error2:\t Id %d does not exist in file", _IdsTXT[i]);
                return false;
-          } 
+          }
           else
           {
                id  = std::distance(unique_ids.begin(), it);
           }
-          
+
           int t =_FramesTXT[i]- _minFrame;
           _peds_t[t].push_back(id);
 
@@ -359,7 +376,7 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
      // obtaining the minimum id and minimum frame
      int maxFrame=0;
      for(TiXmlElement* xFrame = xRootNode->FirstChildElement("frame"); xFrame;
-               xFrame = xFrame->NextSiblingElement("frame"))
+         xFrame = xFrame->NextSiblingElement("frame"))
      {
           int frm = atoi(xFrame->Attribute("ID"));
           if(frm < _minFrame)
@@ -368,10 +385,10 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
           }
           if(frm>maxFrame)
           {
-        	  maxFrame=frm;
+               maxFrame=frm;
           }
           for(TiXmlElement* xAgent = xFrame->FirstChildElement("agent"); xAgent;
-                    xAgent = xAgent->NextSiblingElement("agent"))
+              xAgent = xAgent->NextSiblingElement("agent"))
           {
                int id= atoi(xAgent->Attribute("ID"));
                if(id < _minID)
@@ -388,16 +405,16 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
      vector<int> totalframes;
      for (int i = 0; i <_numPeds; i++)
      {
-    	 totalframes.push_back(0);
+          totalframes.push_back(0);
      }
      //int frameNr=0;
      for(TiXmlElement* xFrame = xRootNode->FirstChildElement("frame"); xFrame;
-               xFrame = xFrame->NextSiblingElement("frame"))
+         xFrame = xFrame->NextSiblingElement("frame"))
      {
-    	  int frameNr = atoi(xFrame->Attribute("ID")) - _minFrame;
+          int frameNr = atoi(xFrame->Attribute("ID")) - _minFrame;
           //todo: can be parallelized with OpenMP
           for(TiXmlElement* xAgent = xFrame->FirstChildElement("agent"); xAgent;
-                    xAgent = xAgent->NextSiblingElement("agent"))
+              xAgent = xAgent->NextSiblingElement("agent"))
           {
                //get agent id, x, y
                double x= atof(xAgent->Attribute("x"));
@@ -406,8 +423,8 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
                int ID= atoi(xAgent->Attribute("ID"))-_minID;
                if(ID>=_numPeds)
                {
-            	   Log->Write("ERROR:\t The number of agents are not corresponding to IDs. Maybe Ped IDs are not continuous in the first frame, please check!!");
-            	   return false;
+                    Log->Write("ERROR:\t The number of agents are not corresponding to IDs. Maybe Ped IDs are not continuous in the first frame, please check!!");
+                    return false;
                }
 
                _peds_t[frameNr].push_back(ID);
@@ -417,15 +434,15 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
                _zCor(ID,frameNr) =  z*M2CM;
                if(_vComponent == "F")
                {
-            	   if(xAgent->Attribute("VD"))
-            	   {
-                        _vComp(ID,frameNr) = *string(xAgent->Attribute("VD")).c_str();
-            	   }
-            	   else
-            	   {
-                        Log->Write("ERROR:\t There is no indicator for velocity component in trajectory file or ini file!!");
-                        return false;
-            	   }
+                    if(xAgent->Attribute("VD"))
+                    {
+                         _vComp(ID,frameNr) = *string(xAgent->Attribute("VD")).c_str();
+                    }
+                    else
+                    {
+                         Log->Write("ERROR:\t There is no indicator for velocity component in trajectory file or ini file!!");
+                         return false;
+                    }
                }
                else
                {
@@ -446,14 +463,14 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
      }
      for(int id = 0; id<_numPeds; id++)
      {
-         int actual_totalframe= totalframes[id];
-         int expect_totalframe=_lastFrame[id]-_firstFrame[id]+1;
-         if(actual_totalframe != expect_totalframe)
-         {
-             Log->Write("Error:\tThe trajectory of ped with ID <%d> is not continuous. Please modify the trajectory file!",id+_minID);
-             Log->Write("Error:\t actual_totalfame = <%d>, expected_totalframe = <%d> ", actual_totalframe, expect_totalframe);
-             return false;
-         }
+          int actual_totalframe= totalframes[id];
+          int expect_totalframe=_lastFrame[id]-_firstFrame[id]+1;
+          if(actual_totalframe != expect_totalframe)
+          {
+               Log->Write("Error:\tThe trajectory of ped with ID <%d> is not continuous. Please modify the trajectory file!",id+_minID);
+               Log->Write("Error:\t actual_totalfame = <%d>, expected_totalframe = <%d> ", actual_totalframe, expect_totalframe);
+               return false;
+          }
      }
 
      return true;
@@ -466,8 +483,8 @@ vector<double> PedData::GetVInFrame(int frame, const vector<int>& ids, double zP
      {
           int id = ids[i];
           int Tpast = frame - _deltaF;
-		  int Tfuture = frame + _deltaF;
-		  double v = GetInstantaneousVelocity1(frame, Tpast, Tfuture, id, _firstFrame, _lastFrame, _xCor, _yCor);
+          int Tfuture = frame + _deltaF;
+          double v = GetInstantaneousVelocity1(frame, Tpast, Tfuture, id, _firstFrame, _lastFrame, _xCor, _yCor);
           if(zPos<1000000.0)
           {
                if(fabs(_zCor(id,frame)-zPos*M2CM)<J_EPS_EVENT)
@@ -488,17 +505,17 @@ vector<double> PedData::GetXInFrame(int frame, const vector<int>& ids, double zP
      vector<double> XInFrame;
      for(int id:ids)
      {
-    	 if(zPos<1000000.0)
-		  {
-                       if(fabs(_zCor(id,frame)-zPos*M2CM)<J_EPS_EVENT)
-			  {
-                               XInFrame.push_back(_xCor(id,frame));
-			  }
-		  }
-		  else
-		  {
-                       XInFrame.push_back(_xCor(id,frame));
-		  }
+          if(zPos<1000000.0)
+          {
+               if(fabs(_zCor(id,frame)-zPos*M2CM)<J_EPS_EVENT)
+               {
+                    XInFrame.push_back(_xCor(id,frame));
+               }
+          }
+          else
+          {
+               XInFrame.push_back(_xCor(id,frame));
+          }
      }
      return XInFrame;
 }
@@ -591,7 +608,7 @@ double PedData::GetInstantaneousVelocity(int Tnow,int Tpast, int Tfuture, int ID
 {
      std::string vcmp = _vComp(ID,Tnow);
      double v=0.0;
-    //check the component used in the calculation of velocity
+     //check the component used in the calculation of velocity
      if(vcmp == "X" || vcmp == "X+"|| vcmp == "X-")
      {
           if((Tpast >=Tfirst[ID])&&(Tfuture <= Tlast[ID]))
@@ -653,41 +670,41 @@ double PedData::GetInstantaneousVelocity1(int Tnow,int Tpast, int Tfuture, int I
      std::string vcmp = _vComp(ID,Tnow);  // the vcmp is the angle from 0 to 360
      if(vcmp=="X+")
      {
-    	 vcmp="0";
+          vcmp="0";
      }
      else if(vcmp=="Y+")
-	 {
-		 vcmp="90";
-	 }
+     {
+          vcmp="90";
+     }
      if(vcmp=="X-")
-	 {
-		 vcmp="180";
-	 }
+     {
+          vcmp="180";
+     }
      if(vcmp=="Y-")
-	 {
-		 vcmp="270";
-	 }
+     {
+          vcmp="270";
+     }
      double v=0.0;
      if(vcmp != "B")  //check the component used in the calculation of velocity
      {
-    	 float alpha=atof(vcmp.c_str())*2*M_PI/360.0;
-
-    	 if((Tpast >=Tfirst[ID])&&(Tfuture <= Tlast[ID]))
-		  {
-                       v=  _fps*CMtoM*((Xcor(ID,Tfuture) - Xcor(ID,Tpast))*cos(alpha)+(Ycor(ID,Tfuture) - Ycor(ID,Tpast))*sin(alpha))/(2.0 * _deltaF);
-		  }
-		  else if((Tpast <Tfirst[ID])&&(Tfuture <= Tlast[ID]))
-		  {
-                       v = _fps*CMtoM*((Xcor(ID,Tfuture) - Xcor(ID,Tnow))*cos(alpha)+(Ycor(ID,Tfuture) - Ycor(ID,Tnow))*sin(alpha))/(_deltaF);  //one dimensional velocity
-		  }
-		  else if((Tpast >=Tfirst[ID])&&(Tfuture > Tlast[ID]))
-		  {
-                       v = _fps*CMtoM*((Xcor(ID,Tnow) - Xcor(ID,Tpast))*cos(alpha)+(Ycor(ID,Tnow) - Ycor(ID,Tpast))*sin(alpha))/( _deltaF);  //one dimensional velocity
-		  }
-		  if(_IgnoreBackwardMovement && v<0)           //if no move back and pedestrian moves back, his velocity is set as 0;
-		  {
-			   v=0;
-		  }
+          float alpha=atof(vcmp.c_str())*2*M_PI/360.0;
+
+          if((Tpast >=Tfirst[ID])&&(Tfuture <= Tlast[ID]))
+          {
+               v=  _fps*CMtoM*((Xcor(ID,Tfuture) - Xcor(ID,Tpast))*cos(alpha)+(Ycor(ID,Tfuture) - Ycor(ID,Tpast))*sin(alpha))/(2.0 * _deltaF);
+          }
+          else if((Tpast <Tfirst[ID])&&(Tfuture <= Tlast[ID]))
+          {
+               v = _fps*CMtoM*((Xcor(ID,Tfuture) - Xcor(ID,Tnow))*cos(alpha)+(Ycor(ID,Tfuture) - Ycor(ID,Tnow))*sin(alpha))/(_deltaF);  //one dimensional velocity
+          }
+          else if((Tpast >=Tfirst[ID])&&(Tfuture > Tlast[ID]))
+          {
+               v = _fps*CMtoM*((Xcor(ID,Tnow) - Xcor(ID,Tpast))*cos(alpha)+(Ycor(ID,Tnow) - Ycor(ID,Tpast))*sin(alpha))/( _deltaF);  //one dimensional velocity
+          }
+          if(_IgnoreBackwardMovement && v<0)           //if no move back and pedestrian moves back, his velocity is set as 0;
+          {
+               v=0;
+          }
 
      }
      else if(vcmp == "B")
@@ -714,23 +731,13 @@ void PedData::CreateGlobalVariables(int numPeds, int numFrames)
      Log->Write("INFO: Enter CreateGlobalVariables with numPeds=%d and numFrames=%d", numPeds, numFrames);
      Log->Write("INFO: allocate memory for xCor");
      _xCor = ub::matrix<double>(numPeds, numFrames);
-     size_t currentSize = getCurrentRSS( )/1000000;
-     size_t peakSize    = getPeakRSS( )/1000000;
-     std::cout << "\tcurrentSize: " << currentSize  <<" (MB)" << std::endl;
-     std::cout << "\tpeakSize: " << peakSize << " (MB)" << std::endl;
      Log->Write("INFO: allocate memory for yCor");
      _yCor = ub::matrix<double>(numPeds, numFrames);
-     std::cout << "\tcurrentSize: " << currentSize  <<" (MB)" << std::endl;
-     std::cout << "\tpeakSize: " << peakSize << " (MB)" << std::endl;
      Log->Write("INFO: allocate memory for zCor");
      _zCor = ub::matrix<double>(numPeds, numFrames);
-     std::cout << "\tcurrentSize: " << currentSize  <<" (MB)" << std::endl;
-     std::cout << "\tpeakSize: " << peakSize << " (MB)" << std::endl;
      Log->Write("INFO: allocate memory for vComp");
      _vComp = ub::matrix<std::string>(numPeds, numFrames);
-     std::cout << "\tcurrentSize: " << currentSize  <<" (MB)" << std::endl;
-     std::cout << "\tpeakSize: " << peakSize << " (MB)" << std::endl;
-     Log->Write(" Finished memory allocation");   
+     Log->Write(" Finished memory allocation");
      _firstFrame = new int[numPeds];  // Record the first frame of each pedestrian
      _lastFrame = new int[numPeds];  // Record the last frame of each pedestrian
      for(int i = 0; i <numPeds; i++) {
diff --git a/methods/VoronoiDiagram.cpp b/methods/VoronoiDiagram.cpp
index 89ca23fdd94ef1b0b752cfea46f84328a0e1d457..a31985b6710a86d2311dd0dc432e6b8ca9ef1834 100644
--- a/methods/VoronoiDiagram.cpp
+++ b/methods/VoronoiDiagram.cpp
@@ -44,9 +44,9 @@ VoronoiDiagram::VoronoiDiagram()
 VoronoiDiagram::~VoronoiDiagram()
 {
 }
-
+//typedef std::vector<std::pair<polygon_2d>, int> >  poly_id;
 // Traversing Voronoi edges using cell iterator.
-vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame, vector<double>& YInFrame,
+std::vector<std::pair<polygon_2d, int> > VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame, vector<double>& YInFrame,
           vector<double>& VInFrame, vector<int>& IdInFrame, const double Bound_Max)
 {
      const int numPedsInFrame = IdInFrame.size();
@@ -59,16 +59,16 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
      {
           points.push_back(point_type2(round(XInFrame[i]), round(YInFrame[i])));
           XInFrame_temp.push_back(round(XInFrame[i]));
-		  YInFrame_temp.push_back(round(YInFrame[i]));
-		  VInFrame_temp.push_back(VInFrame[i]);
-		  IdInFrame_temp.push_back(IdInFrame[i]);
+          YInFrame_temp.push_back(round(YInFrame[i]));
+          VInFrame_temp.push_back(VInFrame[i]);
+          IdInFrame_temp.push_back(IdInFrame[i]);
      }
 
      VD voronoidiagram;
      construct_voronoi(points.begin(), points.end(), &voronoidiagram);
      int Ncell = 0;
-     std::vector<polygon_2d> polygons;
-
+     //std::vector<polygon_2d> polygons;
+     std::vector<std::pair<polygon_2d, int> >    polygons_id;
      double Bd_Box_minX = -Bound_Max;
      double Bd_Box_minY = -Bound_Max;
      double Bd_Box_maxX = Bound_Max;
@@ -83,7 +83,7 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
      for (voronoi_diagram<double>::const_cell_iterator it = voronoidiagram.cells().begin();
                it != voronoidiagram.cells().end(); ++it)
      {
-    	  polygon_2d poly;
+          polygon_2d poly;
           vector<point_type2> polypts;
           point_type2 pt_s;
           point_type2 pt_e;
@@ -170,7 +170,6 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
                }
                edge = edge->next();
           } while (edge != cell.incident_edge());
-          Ncell++;
 
           if (infinite_s && infinite_e)
           {
@@ -190,9 +189,14 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
           }
 
           boost::geometry::correct(poly);
-          polygons.push_back(poly);
+          //polygons.push_back(poly);
+          //cout << "poly is: " << typeid(poly).name() << '\n'
+          int id_ped = IdInFrame[Ncell];
+          std::pair<polygon_2d, int>  poly_id = std::make_pair(poly, id_ped);
+          polygons_id.push_back(poly_id);
+          Ncell++;
      }
-     return polygons;
+     return polygons_id;
 }
 
 point_type2 VoronoiDiagram::retrieve_point(const cell_type& cell)
@@ -359,7 +363,7 @@ vector<point_type2> VoronoiDiagram::add_bounding_points(const point_type2& pt1,
                else
                {
                     double tempx = pt1.x()
-                            		                      + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y());
+                                                              + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y());
                     if (tempx < pt.x())
                     {
                          bounding_vertex.push_back(point_type2(minX, maxY));
@@ -391,7 +395,7 @@ vector<point_type2> VoronoiDiagram::add_bounding_points(const point_type2& pt1,
                else
                {
                     double tempy = pt1.y()
-                            		                      + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x());
+                                                              + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x());
                     if (tempy < pt.y())
                     {
                          bounding_vertex.push_back(point_type2(minX, minY));
@@ -433,26 +437,33 @@ point_type2 VoronoiDiagram::getIntersectionPoint(const point_2d& pt0, const poin
      return interpts;
 }
 
-std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon,
+std::vector<std::pair<polygon_2d, int> > VoronoiDiagram::cutPolygonsWithGeometry(const std::vector<std::pair<polygon_2d, int> > & polygon,
           const polygon_2d& Geometry, const vector<double>& xs, const vector<double>& ys)
 {
-     vector<polygon_2d> intersetionpolygons;
+     // vector<polygon_2d> intersetionpolygons;
+     std::vector<std::pair<polygon_2d, int> > intersetionpolygons;
      int temp = 0;
-     for (const auto & polygon_iterator:polygon)
+     // std::cout << "geometry = " << dsv(Geometry) << "\n";
+     for (const auto & polygon_iterator : polygon)
      {
           polygon_list v;
-          intersection(Geometry, polygon_iterator, v);
+          polygon_2d p = polygon_iterator.first;
+
+          intersection(Geometry, p, v);
+          // std::cout << "p" << polygon_iterator.second << " = " << dsv(p) << "\n";
           BOOST_FOREACH(auto & it, v)
           {
+               // std::cout << "v" << temp << " = " << dsv(it) << "\n";
                if (within(point_2d(xs[temp], ys[temp]), it))
                {
+                    // std::cout << "within " << temp << ": " << xs[temp] << ", " << ys[temp] << "\n";
                     //check and remove duplicates
                     //dispatch::unique (it);
                     polygon_2d simplified;
                     simplify(it,simplified,J_EPS);
 
                     correct(simplified);
-                    intersetionpolygons.push_back(simplified);
+                    intersetionpolygons.push_back(std::make_pair(simplified, polygon_iterator.second));
                }
           }
           temp++;
@@ -461,14 +472,16 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vecto
 }
 
 
-std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon,
+std::vector<std::pair<polygon_2d, int> >  VoronoiDiagram::cutPolygonsWithCircle(const std::vector<std::pair<polygon_2d, int> > & polygon,
           const vector<double>& xs, const vector<double>& ys, double radius, int edges)
 {
-     std::vector<polygon_2d> intersetionpolygons;
+     // std::cout << "ENTER CIRCLE\n";
+     std::vector<std::pair<polygon_2d, int> >  intersetionpolygons;
      int temp = 0;
      for (const auto & polygon_iterator:polygon)
      {
           polygon_2d circle;
+          polygon_2d p = polygon_iterator.first;
           {
                for (int angle = 0; angle <=edges; angle++)
                {
@@ -479,20 +492,27 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<
           }
           correct(circle);
           polygon_list v;
-          intersection(circle, polygon_iterator, v);
+          intersection(circle, p, v);
+          // std::cout << "p" << polygon_iterator.second << " = " << dsv(p) << "\n";
+          // std::cout << "circle" << polygon_iterator.second << " = " << dsv(circle) << "\n";
+
           BOOST_FOREACH(auto & it, v)
           {
+               // std::cout << "v" << temp << " = " << dsv(it) << "\n";
+               // std::cout << "check " << temp << ": " << xs[temp] << ", " << ys[temp] << "\n";
                if (within(point_2d(xs[temp], ys[temp]), it))
                {
+                    // std::cout << "within " << temp << ": " << xs[temp] << ", " << ys[temp] << "\n";
                     //check and remove duplicates
                     //dispatch::unique (it);
                     polygon_2d simplified;
                     simplify(it,simplified,J_EPS);
                     correct(simplified);
-                    intersetionpolygons.push_back(simplified);
+                    intersetionpolygons.push_back(std::make_pair(simplified, polygon_iterator.second));
                }
           }
           temp++;
      }
+     // std::cout << "LEAVE CIRCLE\n";
      return intersetionpolygons;
 }
diff --git a/methods/VoronoiDiagram.h b/methods/VoronoiDiagram.h
index d8f0c58d1f8a172b6d3391c42ab3e98313b4b6f4..88ec2091e796b1ac6d9f5eb07fe6a90ce2521869 100644
--- a/methods/VoronoiDiagram.h
+++ b/methods/VoronoiDiagram.h
@@ -79,10 +79,10 @@ public:
      VoronoiDiagram();
      virtual ~VoronoiDiagram();
 
-     std::vector<polygon_2d> getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame,
+    std::vector<std::pair<polygon_2d, int> > getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame,
                std::vector<int>& IdInFrame, const double Bound_Max);
-     std::vector<polygon_2d> cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys);
-     std::vector<polygon_2d> cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges);
+     std::vector<std::pair<polygon_2d, int> >  cutPolygonsWithGeometry(const std::vector<std::pair<polygon_2d, int> > & polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys);
+     std::vector<std::pair<polygon_2d, int> >  cutPolygonsWithCircle(const std::vector<std::pair<polygon_2d, int> > & polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges);
 };
 
-#endif /* VORONOIDIAGRAM_H_ */
\ No newline at end of file
+#endif /* VORONOIDIAGRAM_H_ */
diff --git a/scripts/_Plot_cell_rho.py b/scripts/_Plot_cell_rho.py
index 5e66b11b03ad26ece80bb6ab0be6637559aadbea..1b1fd034e70277def494230b939ac16e535a4336 100644
--- a/scripts/_Plot_cell_rho.py
+++ b/scripts/_Plot_cell_rho.py
@@ -1,32 +1,39 @@
-import numpy as np
-import matplotlib
-#matplotlib.use('Agg') 
-import matplotlib.pyplot as plt
-from matplotlib.patches import Polygon as pgon
-import Polygon as pol
-import matplotlib.cm as cm
-import pylab
-import pandas as pd
-from mpl_toolkits.axes_grid1 import make_axes_locatable
-import argparse
+# todo https://stackoverflow.com/questions/12881848/draw-polygons-more-efficiently-with-matplotlib
+import os
 import sys
 import logging
+import argparse
 from xml.dom import minidom
-import os
 try:
     import xml.etree.cElementTree as ET
 except ImportError:
     import xml.etree.ElementTree as ET
-    
+
+import numpy as np
+#matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib.patches import Polygon as pgon
+import matplotlib.cm as cm
+
+import Polygon as pol
+#import pylab
+import pandas as pd
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+logfile = "rho_log.txt"
+logging.basicConfig(level=logging.INFO,
+                    format='%(asctime)s - %(levelname)s - %(message)s')
 
 def getParserArgs():
-	parser = argparse.ArgumentParser(description='Combine French data to one file')
-	parser.add_argument("-f", "--filepath", default="./", help='give the path of source file')
-	parser.add_argument("-n", "--namefile", help='give the name of the source file')
-	parser.add_argument("-g", "--geoname", help='give the name of the geometry file')
-	parser.add_argument("-p", "--trajpath", help='give the path of the trajectory file')
-	args = parser.parse_args()
-	return args
+    parser = argparse.ArgumentParser(description='Combine French data to one file')
+    parser.add_argument("-f", "--filepath", default="./", help='give the path of source file')
+    parser.add_argument("-n", "--namefile", help='give the name of the source file')
+    parser.add_argument("-g", "--geoname", help='give the name of the geometry file')
+    parser.add_argument("-p", "--trajpath", help='give the path of the trajectory file')
+    parser.add_argument("-i", "--index", action='store_const', const=True,
+                        default=False, help='plot index of pedestrians along with the Voronoi polygons')
+    args = parser.parse_args()
+    return args
 
 def get_polygon(poly):
     X = []
@@ -52,7 +59,7 @@ def get_geometry_boundary(geometry):
     tree = ET.parse(geometry)
     root = tree.getroot()
     geominX = []
-    geomaxX = [] 
+    geomaxX = []
     geominY = []
     geomaxY = []
     for node in root.iter():
@@ -101,7 +108,7 @@ def plot_geometry(geometry):
 
 def parse_xml_traj_file(filename):
     """
-    parse trajectories in Travisto-format and output results
+    parse trajectories in jpsvis-format and output results
     in the following  format: id    frame    x    y
     (no sorting of the data is performed)
     returns:
@@ -109,16 +116,16 @@ def parse_xml_traj_file(filename):
     N: number of pedestrians
     data: trajectories (numpy.array) [id fr x y]
     """
-    logging.info("parsing <%s>"%filename)
+    logging.info("parsing <%s>", filename)
     try:
         xmldoc = minidom.parse(filename)
     except:
         logging.critical('could not parse file. exit')
-        exit(FAILURE)
+        exit(0)
     N = int(xmldoc.getElementsByTagName('agents')[0].childNodes[0].data)
     fps = xmldoc.getElementsByTagName('frameRate')[0].childNodes[0].data #type unicode
     fps = round(float(fps))
-    
+
     #print ("fps=", fps)
     #fps = int(xmldoc.getElementsByTagName('frameRate')[0].childNodes[0].data)
     logging.info("Npeds = %d, fps = %d"%(N, fps))
@@ -141,89 +148,113 @@ def main():
     filepath = args.filepath
     sys.path.append(filepath)
     namefile = args.namefile
-    geoFile=args.geoname
-    trajpath=args.trajpath
+    geoFile = args.geoname
+    trajpath = args.trajpath
     #geoLocation = filepath.split("Output")[0]
     trajName = namefile.split(".")[0]
     trajType = namefile.split(".")[1].split("_")[0]
     trajFile = os.path.join(trajpath, trajName+"."+trajType) #trajpath+trajName+"."+trajType
-    print(">> plot_cell_rho trajpath %s" %trajpath)
-    print(">> plot_cell_rho  trajName %s"%trajName+"."+trajType)
-    print(">> plot_cell_rho trajFile %s" % trajFile)
+    plotIndex = args.index
+    logging.info("plot_cell_rho trajpath %s", trajpath)
+    logging.info("plot_cell_rho  trajName %s", trajName+"."+trajType)
+    logging.info("plot_cell_rho trajFile %s", trajFile)
     try:
-        frameNr=int(namefile.split("_")[-1])
+        frameNr = int(namefile.split("_")[-1])
     except ValueError:
-        exit("Could not parse fps")
-        
+        exit("ERROR: Could not parse fps")
+
     geominX, geomaxX, geominY, geomaxY = get_geometry_boundary(geoFile)
-    
     fig = plt.figure(figsize=(16*(geomaxX-geominX)/(geomaxY-geominY)+2, 16), dpi=100)
-    ax1 = fig.add_subplot(111,aspect='equal')
+    ax1 = fig.add_subplot(111, aspect='equal')
     plt.rc("font", size=30)
-    plt.rc('pdf',fonttype = 42)
-    ax1.set_yticks([int(1*j) for j in range(-2,5)])
+    plt.rc('pdf', fonttype=42)
+    ax1.set_yticks([int(1*j) for j in range(-2, 5)])
     for label in ax1.get_xticklabels() + ax1.get_yticklabels():
         label.set_fontsize(30)
-        
+
     for tick in ax1.get_xticklines() + ax1.get_yticklines():
         tick.set_markeredgewidth(2)
         tick.set_markersize(6)
         ax1.set_aspect("equal")
 
     plot_geometry(geoFile)
-        
+
     density = np.array([])
     density_orig = np.array([])
-    polygons = [] # polygons converted from string 
-    # TODO use os.path.join
-    polys = open("%s/polygon%s.dat"%(filepath,namefile)).readlines()
+    polygons = [] # polygons converted from string
+    polygon_filename = os.path.join(filepath, "polygon"+namefile+".dat")
+    File = open(polygon_filename)
+    polys = File.readlines()
+    File.close()
+    # plot the blue dots (head) first, to avoid intervening with "text"
+    # ---------------------------------------------------
+    # if  not plotIndex:
+    #     if(trajType == "xml"):
+    #         fps, N, Trajdata = parse_xml_traj_file(trajFile)
+    #     elif(trajType == "txt"):
+    #         try:
+    #             Trajdata = np.array(pd.read_csv(trajFile, comment="#", delimiter=r"\s+"))
+    #         except:
+    #             Trajdata = np.loadtxt(trajFile)
+
+    #     Trajdata = Trajdata[Trajdata[:, 1] == frameNr]
+    #     ax1.plot(Trajdata[:, 2], Trajdata[:, 3], "bo", markersize=20, markeredgewidth=2)
+    # ---------------------------------------------------
+
+    #polys = open("%s/polygon%s.dat"%(filepath,namefile)).readlines()
+    poly_index = []
+    areas = []
     for poly in polys:
-        exec("p = %s"%poly)
+        poly = poly.split("|")
+        poly_index.append(poly[0].strip())
+        Poly = poly[1].strip()
+        exec("p = %s"%Poly)
         pp = locals()['p']
         polygons.append(pp)
-        xx=1.0/pol.Polygon(pp).area()
-        if xx>rho_max:
-            xx=rho_max
-            
-        density=np.append(density,xx)
-        
+        area = pol.Polygon(pp).area()
+        xx = 1.0/area
+        if xx > rho_max:
+            xx = rho_max
+
+        density = np.append(density, xx)
+        areas.append(area)
+
     density_orig = np.copy(density)
-    Maxd=density.max()
-    Mind=density.min()
-    erro=np.ones(np.shape(density))*Mind
-    density=rho_max*(density-erro)/(Maxd-Mind)
+    Maxd = density.max()
+    Mind = density.min()
+    erro = np.ones(np.shape(density))*Mind
+    density = rho_max*(density-erro)/(Maxd-Mind)
     sm = cm.ScalarMappable(cmap=cm.jet)
     sm.set_array(density)
     sm.autoscale_None()
-    sm.set_clim(vmin=0,vmax=5)
+    sm.set_clim(vmin=0, vmax=5)
+    maxArea = np.max(areas)
+    meanArea = np.mean(areas)
     for j, poly in enumerate(polys):
-        ax1.add_patch(pgon(polygons[j],facecolor=sm.to_rgba(density_orig[j]), edgecolor='white',linewidth=2))
-
-    if(trajType=="xml"):
-       fps, N, Trajdata = parse_xml_traj_file(trajFile)
-    elif(trajType=="txt"):
-        try:            
-            Trajdata = np.array(pd.read_csv(trajFile, comment="#", delimiter="\s+"))
-        except :
-            Trajdata = np.loadtxt(trajFile)   
-
-
-    Trajdata = Trajdata[ Trajdata[:, 1] == frameNr ]
-    ax1.plot(Trajdata[:,2], Trajdata[:,3], "bo", markersize = 20, markeredgewidth=2)
-    
-    ax1.set_xlim(geominX-0.2,geomaxX+0.2)
-    ax1.set_ylim(geominY-0.2,geomaxY+0.2)
+        ax1.add_patch(pgon(polygons[j], fc=sm.to_rgba(density_orig[j]), ec='white', lw=2))
+        bcolor = sm.to_rgba(density_orig[j]) #inverse background color
+        icolor = [1 - c for c in bcolor]
+        icolor[-1] = bcolor[-1] # alpha
+        if plotIndex:
+            ax1.text(pol.Polygon(polygons[j]).center()[0],
+                     pol.Polygon(polygons[j]).center()[1],
+                     poly_index[j],
+                     fontsize=25*areas[j]/maxArea, color=icolor)
+
+    ax1.set_xlim(geominX-0.2, geomaxX+0.2)
+    ax1.set_ylim(geominY-0.2, geomaxY+0.2)
     plt.xlabel("x [m]")
     plt.ylabel("y [m]")
     plt.title("%s"%namefile)
     divider = make_axes_locatable(ax1)
     cax = divider.append_axes("right", size="2.5%", pad=0.2)
-    cb=fig.colorbar(sm,ax=ax1,cax=cax,format='%.1f')
-    cb.set_label('Density [$m^{-2}$]') 
-    plt.savefig("%s/rho_%s.png"%(filepath,namefile))
+    cb = fig.colorbar(sm, ax=ax1, cax=cax, format='%.1f')
+    cb.set_label('Density [$m^{-2}$]')
+    figname = os.path.join(filepath, "rho_"+namefile+".png")
+    logging.info("SAVE: %s", figname)
+    plt.savefig(figname)
     plt.close()
 
 
 if __name__ == '__main__':
     main()
-
diff --git a/scripts/_Plot_cell_v.py b/scripts/_Plot_cell_v.py
index cb11ed1d4c076b6f2c15de455c34c7fe85223df2..8657f6555f2a7af753058abd494135688ff0e128 100644
--- a/scripts/_Plot_cell_v.py
+++ b/scripts/_Plot_cell_v.py
@@ -1,90 +1,118 @@
 #!/usr/bin/env python3
-from numpy import *
+import argparse
+import sys
+import os
+import logging
+import numpy as np
+
+import Polygon as pol
+
+import matplotlib.cm as cm
 import matplotlib.pyplot as plt
 from matplotlib.patches import Polygon as pgon
-from Polygon import *         
-import matplotlib.cm as cm
-import pylab
+
+#import pylab
 from mpl_toolkits.axes_grid1 import make_axes_locatable
-import matplotlib
-import argparse
-import sys
-from _Plot_cell_rho import *
 
+from _Plot_cell_rho import get_geometry_boundary
+from _Plot_cell_rho import parse_xml_traj_file
+from _Plot_cell_rho import plot_geometry
 
-def getParserArgs():
-	parser = argparse.ArgumentParser(description='Combine French data to one file')
-	parser.add_argument("-f", "--filepath", default="./", help='give the path of source file')
-	parser.add_argument("-n", "--namefile", help='give the name of the source file')
-	parser.add_argument("-g", "--geoname", help='give the name of the geometry file')
-	parser.add_argument("-p", "--trajpath", help='give the path of the trajectory file')
-	args = parser.parse_args()
-	return args
+logfile = "v_log.txt"
+logging.basicConfig(filename=logfile,
+                    level=logging.INFO,
+                    format='%(asctime)s - %(levelname)s - %(message)s')
 
 
-if __name__ == '__main__':
-   args = getParserArgs()
-   filepath = args.filepath
-   sys.path.append(filepath)
-   namefile = args.namefile
-   geoFile=args.geoname
-   trajpath=args.trajpath
-   geoLocation = filepath.split("Output")[0]
-   trajName = namefile.split(".")[0]
-   trajType = namefile.split(".")[1].split("_")[0]
-   trajFile = os.path.join(trajpath, trajName+"."+trajType)
-   print(">> plot_cell_rho trajpath %s" %trajpath)
-   print(">> plot_cell_rho  trajName %s"%trajName+"."+trajType)
-   print(">> plot_cell_rho trajFile %s" % trajFile)
+def getParserArgs():
+    parser = argparse.ArgumentParser(description='Combine French data to one file')
+    parser.add_argument("-f", "--filepath", default="./", help='give the path of source file')
+    parser.add_argument("-n", "--namefile", help='give the name of the source file')
+    parser.add_argument("-g", "--geoname", help='give the name of the geometry file')
+    parser.add_argument("-p", "--trajpath", help='give the path of the trajectory file')
+    parser.add_argument("-i", "--index",
+                        action='store_const',
+                        const=True,
+                        default=False,
+                        help='plot index of pedestrians along with the Voronoi polygons')
+    arguments = parser.parse_args()
+    return arguments
 
-   frameNr=int(namefile.split("_")[-1])
-   geominX, geomaxX, geominY, geomaxY = get_geometry_boundary(geoFile)
+if __name__ == '__main__':
+    args = getParserArgs()
+    filepath = args.filepath
+    sys.path.append(filepath)
+    namefile = args.namefile
+    geoFile = args.geoname
+    trajpath = args.trajpath
+    geoLocation = filepath.split("Output")[0]
+    trajName = namefile.split(".")[0]
+    trajType = namefile.split(".")[1].split("_")[0]
+    trajFile = os.path.join(trajpath, trajName+"."+trajType)
+    print("plot_cell_v trajpath %s" %trajpath)
+    print("plot_cell_v trajName %s"%trajName+"."+trajType)
+    print("plot_cell_v trajFile %s" % trajFile)
 
-   fig = plt.figure(figsize=(16*(geomaxX-geominX)/(geomaxY-geominY)+2, 16), dpi=100)
-   ax1 = fig.add_subplot(111,aspect='equal')
-   plt.rc("font", size=30)
-   ax1.set_yticks([int(1.00*j) for j in range(-2,5)])
-   for label in ax1.get_xticklabels() + ax1.get_yticklabels():
-      label.set_fontsize(30)
-   for tick in ax1.get_xticklines() + ax1.get_yticklines():
-      tick.set_markeredgewidth(2)
-      tick.set_markersize(6)
-   ax1.set_aspect("equal")
-   plot_geometry(geoFile)
-   
-   velocity=array([])
-   polys = open("%s/polygon%s.dat"%(filepath,namefile)).readlines()
-   velocity=loadtxt("%s/speed%s.dat"%(filepath,namefile))
-   sm = cm.ScalarMappable(cmap=cm.jet)
-   sm.set_array(velocity)
-   sm.autoscale_None()
-   sm.set_clim(vmin=0,vmax=1.5)
-   index=0
-   for poly in polys:
-      exec("p = %s"%poly)
-      xx = velocity[index]
-      index += 1
-      ax1.add_patch(pgon(p,facecolor=sm.to_rgba(xx), edgecolor='white',linewidth=2))
+    frameNr = int(namefile.split("_")[-1])
+    geominX, geomaxX, geominY, geomaxY = get_geometry_boundary(geoFile)
 
-   if(trajType=="xml"):
-       fps, N, Trajdata = parse_xml_traj_file(trajFile)
-   elif(trajType=="txt"):
-        Trajdata = loadtxt(trajFile)   
-   Trajdata = Trajdata[ Trajdata[:, 1] == frameNr ]
-   ax1.plot(Trajdata[:,2], Trajdata[:,3], "bo", markersize = 20, markeredgewidth=2)
-   
-   ax1.set_xlim(geominX-0.2,geomaxX+0.2)
-   ax1.set_ylim(geominY-0.2,geomaxY+0.2)
-   plt.xlabel("x [m]")
-   plt.ylabel("y [m]")
-   plt.title("%s"%namefile)
-   divider = make_axes_locatable(ax1)
-   cax = divider.append_axes("right", size="2.5%", pad=0.2)
-   cb=fig.colorbar(sm,ax=ax1,cax=cax,format='%.1f')
-   cb.set_label('Velocity [$m/s$]') 
-   plt.savefig("%s/v_%s.png"%(filepath,namefile))
-   plt.close()
+    fig = plt.figure(figsize=(16*(geomaxX-geominX)/(geomaxY-geominY)+2, 16), dpi=100)
+    ax1 = fig.add_subplot(111, aspect='equal')
+    plt.rc("font", size=30)
+    ax1.set_yticks([int(1.00*j) for j in range(-2, 5)])
+    for label in ax1.get_xticklabels() + ax1.get_yticklabels():
+        label.set_fontsize(30)
+    for tick in ax1.get_xticklines() + ax1.get_yticklines():
+        tick.set_markeredgewidth(2)
+        tick.set_markersize(6)
+    ax1.set_aspect("equal")
+    plot_geometry(geoFile)
 
+    velocity = np.array([])
+    polygon_filename = os.path.join(filepath, "polygon"+namefile+".dat")
+    File = open(polygon_filename)
+    polys = File.readlines()
+    velocity_filename = os.path.join(filepath, "speed"+namefile+".dat")
+    velocity = np.loadtxt(velocity_filename)
+    sm = cm.ScalarMappable(cmap=cm.jet)
+    sm.set_array(velocity)
+    sm.autoscale_None()
+    sm.set_clim(vmin=0, vmax=1.5)
+    index = 0
+    for poly in polys:
+        poly = poly.split("|")
+        poly_index = poly[0].strip()
+        Poly = poly[1].strip()
+        exec("p = %s"%Poly)
+        xx = velocity[index]
+        index += 1
+        ax1.add_patch(pgon(p,facecolor=sm.to_rgba(xx), edgecolor='white',linewidth=2))
+        # todo
+        #if plotIndex:
+        # ax1.text(pol.Polygon(polygons[j]).center()[0], pol.Polygon(polygons[j]).center()[1], poly_index,
+        #          horizontalalignment='center',
+        #          verticalalignment='center',
+        #          fontsize=20, color='red',
+        #          transform=ax1.transAxes)
 
+    if(trajType == "xml"):
+        fps, N, Trajdata = parse_xml_traj_file(trajFile)
+    elif(trajType == "txt"):
+        Trajdata = np.loadtxt(trajFile)
 
+    Trajdata = Trajdata[Trajdata[:, 1] == frameNr]
+    ax1.plot(Trajdata[:, 2], Trajdata[:, 3], "bo", markersize=20, markeredgewidth=2)
 
+    ax1.set_xlim(geominX-0.2, geomaxX+0.2)
+    ax1.set_ylim(geominY-0.2, geomaxY+0.2)
+    plt.xlabel("x [m]")
+    plt.ylabel("y [m]")
+    plt.title("%s"%namefile)
+    divider = make_axes_locatable(ax1)
+    cax = divider.append_axes("right", size="2.5%", pad=0.2)
+    cb = fig.colorbar(sm, ax=ax1, cax=cax, format='%.1f')
+    cb.set_label('Velocity [$m/s$]')
+    figname = os.path.join(filepath, "v_"+namefile+".png")
+    logging.info("SAVE: %s", figname)
+    plt.savefig(figname)
+    plt.close()
diff --git a/scripts/_Plot_profiles.py b/scripts/_Plot_profiles.py
index 745eb80d0b1cbba91577c3a10e8ce53aa35e0aee..8d5ea6fd4c7837159ac02946ba5a01cf311a6194 100644
--- a/scripts/_Plot_profiles.py
+++ b/scripts/_Plot_profiles.py
@@ -22,10 +22,10 @@ def get_parser_args():
                                                                              'steady state')
     parser.add_argument("-e", "--endsteady", type=int, required=True, help='the frame for the ending of '
                                                                            'steady state')
-    parser.add_argument("-x1", "--geominx", type=float, required=True, help='the minmum x of the geometry')
-    parser.add_argument("-x2", "--geomaxx", type=float, required=True, help='the maxmum x of the geometry')
-    parser.add_argument("-y1", "--geominy", type=float, required=True, help='the minmum y of the geometry')
-    parser.add_argument("-y2", "--geomaxy", type=float, required=True, help='the maxmum y of the geometry')
+    parser.add_argument("-x1", "--geominx", type=float, required=True, help='the minimum x of the geometry')
+    parser.add_argument("-x2", "--geomaxx", type=float, required=True, help='the maximum x of the geometry')
+    parser.add_argument("-y1", "--geominy", type=float, required=True, help='the minimum y of the geometry')
+    parser.add_argument("-y2", "--geomaxy", type=float, required=True, help='the maximum y of the geometry')
     arguments = parser.parse_args()
     return arguments
 
@@ -66,11 +66,14 @@ if __name__ == '__main__':
     ax1 = fig.add_subplot(111, aspect='1')
     plt.rc("font", size=40)
     for j in range(beginsteady, endsteady):
-        density_file = os.path.join(
-            pathfile, "density", "Prf_d_%s
-_id_1_%.5d.dat" % (nametraj, j))
-        print("loading: %s" % density_file)
-        density += np.loadtxt(density_file)
+        density_file = os.path.join(pathfile, "density", "Prf_d_%s_id_1_%.5d.dat" %(nametraj, j))
+        if os.path.exists(density_file):
+            print("loading: %s" % density_file)
+            density += np.loadtxt(density_file)
+            
+        else:
+            print("WARNING: file not found %s"%density_file)
+        
 
     density = density / (endsteady - beginsteady)
 
@@ -97,9 +100,11 @@ _id_1_%.5d.dat" % (nametraj, j))
     for j in range(beginsteady, endsteady):
         velocity_file = os.path.join(
             pathfile, "velocity", "Prf_v_%s_id_1_%.5d.dat" % (nametraj, j))
-        print("loading: %s" % velocity_file)
-        velocity += np.loadtxt(velocity_file)
-
+        if os.path.exists(velocity_file):
+                print("loading: %s" % velocity_file)
+                velocity += np.loadtxt(velocity_file)
+        else:
+            print("WARNING: file not found %s"%velocity_file)
     velocity = velocity / (endsteady - beginsteady)
     im = plt.imshow(velocity,
                     cmap=cm.jet,
diff --git a/scripts/_Plot_timeseries_rho_v.py b/scripts/_Plot_timeseries_rho_v.py
index d89b21bdc1fd0aea91faacbeb6978c6b4460785f..c92d0060cf81a44736ae6991da0f84e6caafebf0 100644
--- a/scripts/_Plot_timeseries_rho_v.py
+++ b/scripts/_Plot_timeseries_rho_v.py
@@ -1,5 +1,5 @@
 #!/usr/bin/env python3
-from numpy import *
+import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
 from Polygon import *         
@@ -24,13 +24,16 @@ def plotRhoT(pathfile,figname,fps,title,data_Classic=None, data_Voronoi=None):
         plt.rc('pdf',fonttype = 42)
         if data_Classic is not None:
                 plt.plot(data_Classic[:,0]/fps,data_Classic[:,1], 'b--', label="Classic method")
+                y_max = np.max(data_Classic[:,1])
         if data_Voronoi is not None:
                 plt.plot(data_Voronoi[:,0]/fps,data_Voronoi[:,1], 'r-', lw=3, label="Voronoi method")
-        plt.xlabel("t [$s$]")
-        plt.ylabel("density [$m^{-2}$]")
+                y_max = np.max(data_Voronoi[:,1])
+
+        plt.xlabel("$t\;\; [s$]", size=18)
+        plt.ylabel(r"$\rho\;\; [m^{-2}$]", size=18)
         plt.gca().set_xlim(left=0)
         #plt.gca().set_ylim(bottom=0)
-        plt.ylim(0,8)
+        plt.ylim(0, y_max + 0.5)
         plt.title("%s"%title)
         plt.legend()
         plt.savefig("%s/%s.png"%(pathfile,figname))
@@ -45,8 +48,8 @@ def plotVT(pathfile,figname,fps,title,data_Classic=None, data_Voronoi=None):
                 plt.plot(data_Classic[:,0]/fps,data_Classic[:,2], 'b--', label="Classic method")
         if data_Voronoi is not None:
                 plt.plot(data_Voronoi[:,0]/fps,data_Voronoi[:,2], 'r-', lw=3, label="Voronoi method")
-        plt.xlabel("t [$s$]")
-        plt.ylabel("velocity [$m/s$]")
+        plt.xlabel("$t\;\; [s]$", size=18)
+        plt.ylabel("$v\;\; [m/s]$", size=18)
         plt.gca().set_xlim(left=0)
         #plt.gca().set_ylim(bottom=0)
         plt.ylim(0,2)
@@ -70,7 +73,7 @@ if __name__ == '__main__':
                    figname_rho="rho_t_%s_id_%s"%(nametraj,ID)
                    figname_v="v_t_%s_id_%s"%(nametraj,ID)
                    title = "%s_id_%s"%(nametraj,ID)
-                   data_Classic = loadtxt(fC)
+                   data_Classic = np.loadtxt(fC)
                    plotRhoT(pathfile,figname_rho,fps,title,data_Classic)
                    plotVT(pathfile,figname_v,fps,title,data_Classic)
            
@@ -80,10 +83,10 @@ if __name__ == '__main__':
                    figname_rho="rho_t_%s_id_%s"%(nametraj,ID)
                    figname_v="v_t_%s_id_%s"%(nametraj,ID)
                    title = "%s_id_%s"%(nametraj,ID)
-                   data_Voronoi = loadtxt(fV)
+                   data_Voronoi = np.loadtxt(fV)
                    fC = "%s/rho_v_Classic_%s_id_%s.dat"%(pathfile,nametraj,ID)
                    if (os.path.isfile(fC)):
-                           data_Classic = loadtxt(fC)
+                           data_Classic = np.loadtxt(fC)
                            plotRhoT(pathfile,figname_rho,fps,title,data_Classic,data_Voronoi)
                            plotVT(pathfile,figname_v,fps,title,data_Classic,data_Voronoi)
                    else: