diff --git a/.travis.yml b/.travis.yml
index 5cc25d9..55e1471 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -6,20 +6,17 @@ language: cpp
 #env:
 
 before_install:
-  - docker pull ubuntu:xenial
-  - docker run -itd -v "${TRAVIS_BUILD_DIR}:/travis/workdir" --name ubuntu-build ubuntu:xenial
+  - docker pull fdiblen/ubuntu1604-test:latest
+  - docker run -itd -v "${TRAVIS_BUILD_DIR}:/travis/workdir" --name ubuntu-build fdiblen/ubuntu1604-test:latest
+  - docker pull fdiblen/sl7-test:latest
+  - docker run -itd -v "${TRAVIS_BUILD_DIR}:/travis/workdir" --name sl-build fdiblen/sl7-test:latest
 #  - docker pull base/archlinux
 #  - docker run -itd -v "${TRAVIS_BUILD_DIR}:/travis/workdir" --name arch-build base/archlinux
-#  - docker pull scientificlinux/sl:7
-#  - docker run -itd -v "${TRAVIS_BUILD_DIR}:/travis/workdir" --name sl-build scientificlinux/sl:7
 
 script:
   - export BRANCH=$(if [ "$TRAVIS_PULL_REQUEST" == "false" ]; then echo $TRAVIS_BRANCH; else echo $TRAVIS_PULL_REQUEST_BRANCH; fi)
   - echo "TRAVIS_BRANCH=$TRAVIS_BRANCH, PR=$PR, BRANCH=$BRANCH"
-  - docker exec --env BRANCH=$BRANCH ubuntu-build /bin/bash -c "/travis/workdir/build-tests/ubuntu/prepare-ubuntu.sh"
-  - docker exec --env BRANCH=$BRANCH ubuntu-build /bin/bash -c "/travis/workdir/build-tests/ubuntu/compile-ubuntu.sh"
+  - docker exec -e BRANCH=$BRANCH -e IMAGE='ubuntu' ubuntu-build /bin/bash -c "/travis/workdir/build-tests/compile_sagecal.sh"
+  - docker exec -e BRANCH=$BRANCH -e IMAGE='sl7' sl-build /bin/bash -c "/travis/workdir/build-tests/compile_sagecal.sh"
 #  - docker exec --env BRANCH=$BRANCH arch-build /bin/bash -c "/travis/workdir/build-tests/arch/prepare-arch.sh"
 #  - docker exec --env BRANCH=$BRANCH arch-build /bin/bash -c "/travis/workdir/build-tests/arch/compile-arch.sh"
-#  - docker exec --env BRANCH=$BRANCH sl-build /bin/bash -c "/travis/workdir/build-tests/sl/prepare-sl.sh"
-#  - docker exec --env BRANCH=$BRANCH sl-build /bin/bash -c "/travis/workdir/build-tests/sl/compile-sl.sh"
-
diff --git a/CMakeLists.txt b/CMakeLists.txt
index aedcefa..34a923f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,14 +50,15 @@ IF("${isSystemDir}" STREQUAL "-1")
    SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
 ENDIF("${isSystemDir}" STREQUAL "-1")
 
-# libsynthesis has mix of C++ and Fortran which is not handled well by
-# versions before 2.8.
-if(${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION} LESS 2.8)
-  if(NOT LIB_EXTRA_SYNTHESIS)
-    set(LIB_EXTRA_SYNTHESIS gfortran)
-  endif(NOT LIB_EXTRA_SYNTHESIS)
-endif(${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION} LESS 2.8)
 
+# # libsynthesis has mix of C++ and Fortran which is not handled well by
+# # versions before 2.8.
+# if(${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION} LESS 2.8)
+#   if(NOT LIB_EXTRA_SYNTHESIS)
+#     set(LIB_EXTRA_SYNTHESIS gfortran)
+#   endif(NOT LIB_EXTRA_SYNTHESIS)
+# endif(${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION} LESS 2.8)
+#
 
 
 #--------------------------------------- SageCal dependencies
@@ -73,7 +74,7 @@ include_directories(${CASACORE_INCLUDE_DIR})
 
 #cfitsio
 find_package(CfitsIO REQUIRED)
-include_directories(${CFITSIO_INCLUDE_DIR})
+include_directories(${CFITSIO_INCLUDE})
 
 #lapack
 find_package(LAPACK REQUIRED)
@@ -91,10 +92,21 @@ include_directories(${OpenBLAS_INCLUDE_DIR})
 #find_package(BLAS REQUIRED)
 #include_directories(${BLAS_INCLUDE_DIR})
 
-# FIXME: is this really needed?
-##gfortran
+# ##gfortran
+enable_language(Fortran)
 #find_package(GFortranLibs REQUIRED)
-#include_directories(${GFORTRAN_INCLUDE_DIR})
+# include_directories(${GFORTRAN_INCLUDE_DIR})
+
+
+
+
+
+# enable_language(Fortran)
+# # get_filename_component (Fortran_COMPILER_NAME ${CMAKE_Fortran_COMPILER} NAME)
+# set(FC "$ENV{FC}" CACHE INTERNAL "Got from environment variable")
+# set(CMAKE_Fortran_COMPILER "$ENV{FC}" CACHE INTERNAL "Got from environment variable")
+#
+
 
 ##hdf5
 #find_package(HDF5 REQUIRED)
@@ -134,6 +146,7 @@ set (CMAKE_CXX_FLAGS "-std=c++0x -g -O3 -Wall")
 #--------------------------------------- summary
 message(STATUS "\n############################\n#   Configuration summary\n############################")
 message (STATUS "CMAKE_SYSTEM .......... = ${CMAKE_SYSTEM}")
+message (STATUS "CMAKE_INSTALL_PREFIX .. = ${CMAKE_INSTALL_PREFIX}")
 message (STATUS "CMAKE_BUILD_TYPE ...... = ${CMAKE_BUILD_TYPE}")
 message (STATUS "BUILD_SHARED_LIBS ..... = ${BUILD_SHARED_LIBS}")
 message (STATUS "CMAKE_INSTALL_NAME_DIR  = ${CMAKE_INSTALL_NAME_DIR}")
@@ -144,30 +157,36 @@ message (STATUS "CMAKE_CXX_FLAGS ....... = ${CMAKE_CXX_FLAGS}")
 message (STATUS "ENABLE_MPI ............ = ${ENABLE_MPI}")
 message (STATUS "ENABLE_CUDA ........... = ${ENABLE_CUDA}")
 
-message (STATUS "CASACORE_INCLUDE_DIR........... = ${CASACORE_INCLUDE_DIR}")
-message (STATUS "CASACORE_LIBRARIES........... = ${CASACORE_LIBRARIES}")
+message (STATUS "CASACORE_INCLUDE_DIR....= ${CASACORE_INCLUDE_DIR}")
+message (STATUS "CASACORE_LIBRARIES......= ${CASACORE_LIBRARIES}")
 
 message (STATUS "OpenBLAS_LIB .......... = ${OpenBLAS_LIB}")
 
-message (STATUS "GLIB_PKG_INCLUDE_DIRS.............. = ${GLIB_PKG_INCLUDE_DIRS}")
-message (STATUS "GLIB_PKG_LIBRARIES............. = ${GLIB_PKG_LIBRARIES}")
+message (STATUS "GLIB_PKG_INCLUDE_DIRS...= ${GLIB_PKG_INCLUDE_DIRS}")
+message (STATUS "GLIB_PKG_LIBRARIES......= ${GLIB_PKG_LIBRARIES}")
 
 #message (STATUS "LAPACK-INC ............ = ${LAPACK_INCLUDE_DIR}")
-message (STATUS "LAPACK_LIBRARIES........... = ${LAPACK_LIBRARIES}")
+message (STATUS "LAPACK_LIBRARIES........= ${LAPACK_LIBRARIES}")
 
-#message (STATUS "GFORTRAN-INC  ......... = ${GFORTRAN_INCLUDE_DIR}")
-#message (STATUS "GFORTRAN-LIBS ......... = ${LIBGFORTRAN_LIBRARIES}")
+# message (STATUS "GFORTRAN-INC  ......... = ${GFORTRAN_INCLUDE_DIR}")
+# message (STATUS "GFORTRAN-LIBS ......... = ${LIBGFORTRAN_LIBRARIES}")
+message (STATUS "CMAKE_Fortran_COMPILER  = ${CMAKE_Fortran_COMPILER}")
 
-message (STATUS "CFITSIO_INCLUDE_DIR.......... = ${CFITSIO_INCLUDE_DIR}")
-message (STATUS "CFITSIO_LIBRARIES.......... = ${CFITSIO_LIBRARIES}")
+message (STATUS "CFITSIO_ROOT_DIR........= ${CFITSIO_ROOT_DIR}")
+message (STATUS "CFITSIO_INCLUDE.........= ${CFITSIO_INCLUDE}")
+message (STATUS "CFITSIO_LIB.............= ${CFITSIO_LIB}")
 
-message (STATUS "WCSLIB_INCLUDE_DIR........... = ${WCSLIB_INCLUDE_DIR}")
-message (STATUS "WCSLIB_LIBRARIES ........... = ${WCSLIB_LIBRARIES}")
+message (STATUS "WCSLIB_INCLUDE_DIRS.....= ${WCSLIB_INCLUDE_DIRS}")
+message (STATUS "WCSLIB_LIBRARIES .......= ${WCSLIB_LIBRARIES}")
 
-message (STATUS "HDF5_INCLUDE_DIR........... = ${HDF5_INCLUDE_DIRS}")
-message (STATUS "HDF5_LIBRARIES........... = ${HDF5_LIBRARIES}")
+message (STATUS "HDF5_INCLUDE_DIR........= ${HDF5_INCLUDE_DIRS}")
+message (STATUS "HDF5_LIBRARIES..........= ${HDF5_LIBRARIES}")
 
 
+message (STATUS "CMAKE_CURRENT_BINARY_DIR..........= ${CMAKE_CURRENT_BINARY_DIR}")
+message (STATUS "CMAKE_CURRENT_LIST_DIR..........= ${CMAKE_CURRENT_LIST_DIR}")
+message (STATUS "CMAKE_INSTALL_PREFIX..........= ${CMAKE_INSTALL_PREFIX}")
+
 
 #--------------------------------------- include directories
 add_subdirectory(src)
@@ -178,3 +197,5 @@ add_subdirectory(src)
 # FIXME: this will be the final step for testing
 #file(COPY ${PROJECT_SOURCE_DIR}/test DESTINATION ${MAINFOLDER}/dist/test)
 
+
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/dist/ DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/CMakeModules/FindCasaCore.cmake b/CMakeModules/FindCasaCore.cmake
index 6c6e1cd..0939e6e 100644
--- a/CMakeModules/FindCasaCore.cmake
+++ b/CMakeModules/FindCasaCore.cmake
@@ -71,7 +71,15 @@
 # directly fed to the linker.
 #
 #   Usage: casacore_resolve_dependencies(result components...)
-#
+
+
+if (NOT "$ENV{CASACORE_ROOT_DIR}" STREQUAL "")
+set(CASACORE_ROOT_DIR "$ENV{CASACORE_ROOT_DIR}" CACHE INTERNAL "Got from environment variable")
+endif()
+
+
+
+
 macro(casacore_resolve_dependencies _result)
   set(${_result} ${ARGN})
   set(_index 0)
diff --git a/CMakeModules/FindCfitsIO.cmake b/CMakeModules/FindCfitsIO.cmake
index a749701..cbfda03 100644
--- a/CMakeModules/FindCfitsIO.cmake
+++ b/CMakeModules/FindCfitsIO.cmake
@@ -3,11 +3,11 @@
 #  CFITSIO_ROOT_DIR     - CFITSIO root directory
 # Variables defined by this module:
 #  CFITSIO_FOUND        - system has CFITSIO
-#  CFITSIO_INCLUDE_DIR  - the CFITSIO include directory (cached)
-#  CFITSIO_INCLUDE_DIRS - the CFITSIO include directories
-#                         (identical to CFITSIO_INCLUDE_DIR)
+#  CFITSIO_INCLUDE  - the CFITSIO include directory (cached)
+#  CFITSIO_INCLUDES - the CFITSIO include directories
+#                         (identical to CFITSIO_INCLUDE)
 #  CFITSIO_LIBRARY      - the CFITSIO library (cached)
-#  CFITSIO_LIBRARIES    - the CFITSIO libraries
+#  CFITSIO_LIB    - the CFITSIO libraries
 #                         (identical to CFITSIO_LIBRARY)
 #  CFITSIO_VERSION_STRING the found version of CFITSIO, padded to 3 digits
 
@@ -33,12 +33,17 @@
 
 if(NOT CFITSIO_FOUND)
 
-  find_path(CFITSIO_INCLUDE_DIR fitsio.h
+  if (NOT "$ENV{CFITSIO_ROOT_DIR}" STREQUAL "")
+  set(CFITSIO_ROOT_DIR "$ENV{CFITSIO_ROOT_DIR}" CACHE INTERNAL "Got from environment variable")
+  endif()
+
+
+  find_path(CFITSIO_INCLUDE fitsio.h
     HINTS ${CFITSIO_ROOT_DIR} PATH_SUFFIXES include include/cfitsio
              include/libcfitsio0)
 
-  if(CFITSIO_INCLUDE_DIR)
-    FILE(READ "${CFITSIO_INCLUDE_DIR}/fitsio.h" CFITSIO_H)
+  if(CFITSIO_INCLUDE)
+    FILE(READ "${CFITSIO_INCLUDE}/fitsio.h" CFITSIO_H)
     set(CFITSIO_VERSION_REGEX ".*#define CFITSIO_VERSION[^0-9]*([0-9]+)\\.([0-9]+).*")
     if ("${CFITSIO_H}" MATCHES ${CFITSIO_VERSION_REGEX})
       # Pad CFITSIO minor version to three digit because 3.181 is older than 3.35
@@ -51,24 +56,24 @@ if(NOT CFITSIO_FOUND)
     else ()
       set(CFITSIO_VERSION_STRING "Unknown")
     endif()
-  endif(CFITSIO_INCLUDE_DIR)
+  endif(CFITSIO_INCLUDE)
 
   find_library(CFITSIO_LIBRARY cfitsio
     HINTS ${CFITSIO_ROOT_DIR} PATH_SUFFIXES lib)
   find_library(M_LIBRARY m)
-  mark_as_advanced(CFITSIO_INCLUDE_DIR CFITSIO_LIBRARY M_LIBRARY)
+  mark_as_advanced(CFITSIO_INCLUDE CFITSIO_LIBRARY M_LIBRARY)
 
   if(CMAKE_VERSION VERSION_LESS "2.8.3")
     find_package_handle_standard_args(CFITSIO DEFAULT_MSG
-      CFITSIO_LIBRARY M_LIBRARY CFITSIO_INCLUDE_DIR)
+      CFITSIO_LIBRARY M_LIBRARY CFITSIO_INCLUDE)
   else ()
     include(FindPackageHandleStandardArgs)
     find_package_handle_standard_args(CFITSIO
-      REQUIRED_VARS CFITSIO_LIBRARY M_LIBRARY CFITSIO_INCLUDE_DIR
+      REQUIRED_VARS CFITSIO_LIBRARY M_LIBRARY CFITSIO_INCLUDE
       VERSION_VAR CFITSIO_VERSION_STRING)
   endif ()
 
-  set(CFITSIO_INCLUDE_DIRS ${CFITSIO_INCLUDE_DIR})
-  set(CFITSIO_LIBRARIES ${CFITSIO_LIBRARY} ${M_LIBRARY})
+  set(CFITSIO_INCLUDES ${CFITSIO_INCLUDE})
+  set(CFITSIO_LIB ${CFITSIO_LIBRARY} ${M_LIBRARY})
 
 endif(NOT CFITSIO_FOUND)
diff --git a/CMakeModules/FindGFortranLibs.cmake b/CMakeModules/FindGFortranLibs.cmake
index 36b42bf..2e7d9e3 100644
--- a/CMakeModules/FindGFortranLibs.cmake
+++ b/CMakeModules/FindGFortranLibs.cmake
@@ -24,120 +24,128 @@ if(NOT CMAKE_REQUIRED_QUIET)
   message(STATUS "Looking for gfortran related libraries...")
 endif()
 
+# enable_language(Fortran)
+# if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
+#
+#   # Basically, call "gfortran -v" to dump compiler info to the string
+#   # GFORTRAN_VERBOSE_STR, which will be used to get necessary paths
+#   message(STATUS "Extracting library and header information by calling 'gfortran -v'...")
+#   execute_process(COMMAND "${CMAKE_Fortran_COMPILER}" "-v" ERROR_VARIABLE
+#     GFORTRAN_VERBOSE_STR RESULT_VARIABLE FLAG)
+#
+#   # For debugging
+#   message(STATUS "'gfortran -v' returned:")
+#   message(STATUS "${GFORTRAN_VERBOSE_STR}")
+#
+#   # Detect gfortran version
+#   string(REGEX MATCH "gcc version [^\t\n ]+" GFORTRAN_VER_STR "${GFORTRAN_VERBOSE_STR}")
+#   string(REGEX REPLACE "gcc version ([^\t\n ]+)" "\\1" GFORTRAN_VERSION_STRING "${GFORTRAN_VER_STR}")
+#   message(STATUS "Detected gfortran version ${GFORTRAN_VERSION_STRING}")
+#   unset(GFORTRAN_VER_STR)
+#
+#   set(MATCH_REGEX "[^\t\n ]+[\t\n ]+")
+#   set(REPLACE_REGEX "([^\t\n ]+)")
+#
+#   # Find architecture for compiler
+#   string(REGEX MATCH "Target: [^\t\n ]+"
+#     GFORTRAN_ARCH_STR "${GFORTRAN_VERBOSE_STR}")
+#   message(STATUS "Architecture string: ${GFORTRAN_ARCH_STR}")
+#   string(REGEX REPLACE "Target: ([^\t\n ]+)" "\\1"
+#     GFORTRAN_ARCH "${GFORTRAN_ARCH_STR}")
+#   message(STATUS "Detected gfortran architecture: ${GFORTRAN_ARCH}")
+#   unset(GFORTRAN_ARCH_STR)
+#
+#   # Find install prefix, if it exists; if not, use default
+#   string(REGEX MATCH  "--prefix=[^\t\n ]+[\t\n ]+"
+#     GFORTRAN_PREFIX_STR "${GFORTRAN_VERBOSE_STR}")
+#   if(NOT GFORTRAN_PREFIX_STR)
+#     message(STATUS "Detected default gfortran prefix")
+#     set(GFORTRAN_PREFIX_DIR "/usr/local") # default prefix for gcc install
+#   else()
+#     string(REGEX REPLACE "--prefix=([^\t\n ]+)" "\\1"
+#       GFORTRAN_PREFIX_DIR "${GFORTRAN_PREFIX_STR}")
+#   endif()
+#   message(STATUS "Detected gfortran prefix: ${GFORTRAN_PREFIX_DIR}")
+#   unset(GFORTRAN_PREFIX_STR)
+#
+#   # Find install exec-prefix, if it exists; if not, use default
+#   string(REGEX MATCH "--exec-prefix=[^\t\n ]+[\t\n ]+" "\\1"
+#     GFORTRAN_EXEC_PREFIX_STR "${GFORTRAN_VERBOSE_STR}")
+#   if(NOT GFORTRAN_EXEC_PREFIX_STR)
+#     message(STATUS "Detected default gfortran exec-prefix")
+#     set(GFORTRAN_EXEC_PREFIX_DIR "${GFORTRAN_PREFIX_DIR}")
+#   else()
+#     string(REGEX REPLACE "--exec-prefix=([^\t\n ]+)" "\\1"
+#       GFORTRAN_EXEC_PREFIX_DIR "${GFORTRAN_EXEC_PREFIX_STR}")
+#   endif()
+#   message(STATUS "Detected gfortran exec-prefix: ${GFORTRAN_EXEC_PREFIX_DIR}")
+#   UNSET(GFORTRAN_EXEC_PREFIX_STR)
+#
+#   # Find library directory and include directory, if library directory specified
+#   string(REGEX MATCH "--libdir=[^\t\n ]+"
+#     GFORTRAN_LIB_DIR_STR "${GFORTRAN_VERBOSE_STR}")
+#   if(NOT GFORTRAN_LIB_DIR_STR)
+#     message(STATUS "Found --libdir flag -- not found")
+#     message(STATUS "Using default gfortran library & include directory paths")
+#     set(GFORTRAN_LIBRARIES_DIR
+#       "${GFORTRAN_EXEC_PREFIX_DIR}/lib/gcc/${GFORTRAN_ARCH}/${GFORTRAN_VERSION_STRING}")
+#     string(CONCAT GFORTRAN_INCLUDE_DIR "${GFORTRAN_LIBRARIES_DIR}" "/include")
+#   else()
+#     message(STATUS "Found --libdir flag -- yes")
+#     string(REGEX REPLACE "--libdir=([^\t\n ]+)" "\\1"
+#       GFORTRAN_LIBRARIES_DIR "${GFORTRAN_LIB_DIR_STR}")
+#     string(CONCAT GFORTRAN_INCLUDE_DIR "${GFORTRAN_LIBRARIES_DIR}" "/gcc/" "${GFORTRAN_ARCH}" "/" "${GFORTRAN_VERSION_STRING}" "/include")
+#   endif()
+#   message(STATUS "gfortran libraries path: ${GFORTRAN_LIBRARIES_DIR}")
+#   message(STATUS "gfortran include path dir: ${GFORTRAN_INCLUDE_DIR}")
+#   unset(GFORTRAN_LIB_DIR_STR)
+#
+#   # There are lots of other build options for gcc & gfortran. For now, the
+#   # options implemented above should cover a lot of common use cases.
+#
+#   # Clean up be deleting the output string from "gfortran -v"
+#   unset(GFORTRAN_VERBOSE_STR)
+#
+#   # Find paths for libgfortran, libquadmath, libgomp
+#   # libgomp needed for OpenMP support without Clang
+#   find_library(LIBGFORTRAN_LIBRARIES NAMES gfortran libgfortran
+#     HINTS ${GFORTRAN_LIBRARIES_DIR})
+#   find_library(LIBQUADMATH_LIBRARIES NAMES quadmath libquadmath
+#     HINTS ${GFORTRAN_LIBRARIES_DIR})
+#   find_library(LIBGOMP_LIBRARIES NAMES gomp libgomp
+#     HINTS ${GFORTRAN_LIBRARIES_DIR})
+#
+#   # Find OpenMP headers
+#   find_path(LIBGOMP_INCLUDE_DIR NAMES omp.h HINTS ${GFORTRAN_INCLUDE_DIR})
+#
+# else()
+#   message(STATUS "CMAKE_Fortran_COMPILER_ID does not match 'GNU'!")
+# endif()
+#
+# include(FindPackageHandleStandardArgs)
+#
+# # Required: libgfortran, libquadmath, path for gfortran libraries
+# # Optional: libgomp, path for OpenMP headers, path for gcc/gfortran headers
+# find_package_handle_standard_args(GFortranLibs
+#   REQUIRED_VARS LIBGFORTRAN_LIBRARIES LIBQUADMATH_LIBRARIES GFORTRAN_LIBRARIES_DIR
+#   VERSION_VAR GFORTRAN_VERSION_STRING)
+#
+# if(GFORTRANLIBS_FOUND)
+#   message(STATUS "Looking for gfortran libraries -- found")
+#   message(STATUS "gfortran version: ${GFORTRAN_VERSION_STRING}")
+# else()
+#   message(STATUS "Looking for gfortran libraries -- not found")
+# endif()
+#
+# mark_as_advanced(LIBGFORTRAN_LIBRARIES LIBQUADMATH_LIBRARIES
+#   LIBGOMP_LIBRARIES LIBGOMP_INCLUDE_DIR
+#   GFORTRAN_LIBRARIES_DIR GFORTRAN_INCLUDE_DIR)
+# # FindGFortranLIBS.cmake ends here
+
+
 enable_language(Fortran)
-if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
-
-  # Basically, call "gfortran -v" to dump compiler info to the string
-  # GFORTRAN_VERBOSE_STR, which will be used to get necessary paths
-  message(STATUS "Extracting library and header information by calling 'gfortran -v'...")
-  execute_process(COMMAND "${CMAKE_Fortran_COMPILER}" "-v" ERROR_VARIABLE
-    GFORTRAN_VERBOSE_STR RESULT_VARIABLE FLAG)
-
-  # For debugging
-  message(STATUS "'gfortran -v' returned:")
-  message(STATUS "${GFORTRAN_VERBOSE_STR}")
-
-  # Detect gfortran version
-  string(REGEX MATCH "gcc version [^\t\n ]+" GFORTRAN_VER_STR "${GFORTRAN_VERBOSE_STR}")
-  string(REGEX REPLACE "gcc version ([^\t\n ]+)" "\\1" GFORTRAN_VERSION_STRING "${GFORTRAN_VER_STR}")
-  message(STATUS "Detected gfortran version ${GFORTRAN_VERSION_STRING}")
-  unset(GFORTRAN_VER_STR)
-
-  set(MATCH_REGEX "[^\t\n ]+[\t\n ]+")
-  set(REPLACE_REGEX "([^\t\n ]+)")
-
-  # Find architecture for compiler
-  string(REGEX MATCH "Target: [^\t\n ]+"
-    GFORTRAN_ARCH_STR "${GFORTRAN_VERBOSE_STR}")
-  message(STATUS "Architecture string: ${GFORTRAN_ARCH_STR}")
-  string(REGEX REPLACE "Target: ([^\t\n ]+)" "\\1"
-    GFORTRAN_ARCH "${GFORTRAN_ARCH_STR}")
-  message(STATUS "Detected gfortran architecture: ${GFORTRAN_ARCH}")
-  unset(GFORTRAN_ARCH_STR)
-
-  # Find install prefix, if it exists; if not, use default
-  string(REGEX MATCH  "--prefix=[^\t\n ]+[\t\n ]+"
-    GFORTRAN_PREFIX_STR "${GFORTRAN_VERBOSE_STR}")
-  if(NOT GFORTRAN_PREFIX_STR)
-    message(STATUS "Detected default gfortran prefix")
-    set(GFORTRAN_PREFIX_DIR "/usr/local") # default prefix for gcc install
-  else()
-    string(REGEX REPLACE "--prefix=([^\t\n ]+)" "\\1"
-      GFORTRAN_PREFIX_DIR "${GFORTRAN_PREFIX_STR}")
-  endif()
-  message(STATUS "Detected gfortran prefix: ${GFORTRAN_PREFIX_DIR}")
-  unset(GFORTRAN_PREFIX_STR)
-
-  # Find install exec-prefix, if it exists; if not, use default
-  string(REGEX MATCH "--exec-prefix=[^\t\n ]+[\t\n ]+" "\\1"
-    GFORTRAN_EXEC_PREFIX_STR "${GFORTRAN_VERBOSE_STR}")
-  if(NOT GFORTRAN_EXEC_PREFIX_STR)
-    message(STATUS "Detected default gfortran exec-prefix")
-    set(GFORTRAN_EXEC_PREFIX_DIR "${GFORTRAN_PREFIX_DIR}")
-  else()
-    string(REGEX REPLACE "--exec-prefix=([^\t\n ]+)" "\\1"
-      GFORTRAN_EXEC_PREFIX_DIR "${GFORTRAN_EXEC_PREFIX_STR}")
-  endif()
-  message(STATUS "Detected gfortran exec-prefix: ${GFORTRAN_EXEC_PREFIX_DIR}")
-  UNSET(GFORTRAN_EXEC_PREFIX_STR)
-
-  # Find library directory and include directory, if library directory specified
-  string(REGEX MATCH "--libdir=[^\t\n ]+"
-    GFORTRAN_LIB_DIR_STR "${GFORTRAN_VERBOSE_STR}")
-  if(NOT GFORTRAN_LIB_DIR_STR)
-    message(STATUS "Found --libdir flag -- not found")
-    message(STATUS "Using default gfortran library & include directory paths")
-    set(GFORTRAN_LIBRARIES_DIR
-      "${GFORTRAN_EXEC_PREFIX_DIR}/lib/gcc/${GFORTRAN_ARCH}/${GFORTRAN_VERSION_STRING}")
-    string(CONCAT GFORTRAN_INCLUDE_DIR "${GFORTRAN_LIBRARIES_DIR}" "/include")
-  else()
-    message(STATUS "Found --libdir flag -- yes")
-    string(REGEX REPLACE "--libdir=([^\t\n ]+)" "\\1"
-      GFORTRAN_LIBRARIES_DIR "${GFORTRAN_LIB_DIR_STR}")
-    string(CONCAT GFORTRAN_INCLUDE_DIR "${GFORTRAN_LIBRARIES_DIR}" "/gcc/" "${GFORTRAN_ARCH}" "/" "${GFORTRAN_VERSION_STRING}" "/include")
-  endif()
-  message(STATUS "gfortran libraries path: ${GFORTRAN_LIBRARIES_DIR}")
-  message(STATUS "gfortran include path dir: ${GFORTRAN_INCLUDE_DIR}")
-  unset(GFORTRAN_LIB_DIR_STR)
-
-  # There are lots of other build options for gcc & gfortran. For now, the
-  # options implemented above should cover a lot of common use cases.
-
-  # Clean up be deleting the output string from "gfortran -v"
-  unset(GFORTRAN_VERBOSE_STR)
-
-  # Find paths for libgfortran, libquadmath, libgomp
-  # libgomp needed for OpenMP support without Clang
-  find_library(LIBGFORTRAN_LIBRARIES NAMES gfortran libgfortran
-    HINTS ${GFORTRAN_LIBRARIES_DIR})
-  find_library(LIBQUADMATH_LIBRARIES NAMES quadmath libquadmath
-    HINTS ${GFORTRAN_LIBRARIES_DIR})
-  find_library(LIBGOMP_LIBRARIES NAMES gomp libgomp
-    HINTS ${GFORTRAN_LIBRARIES_DIR})
-
-  # Find OpenMP headers
-  find_path(LIBGOMP_INCLUDE_DIR NAMES omp.h HINTS ${GFORTRAN_INCLUDE_DIR})
-
-else()
-  message(STATUS "CMAKE_Fortran_COMPILER_ID does not match 'GNU'!")
-endif()
-
-include(FindPackageHandleStandardArgs)
-
-# Required: libgfortran, libquadmath, path for gfortran libraries
-# Optional: libgomp, path for OpenMP headers, path for gcc/gfortran headers
-find_package_handle_standard_args(GFortranLibs
-  REQUIRED_VARS LIBGFORTRAN_LIBRARIES LIBQUADMATH_LIBRARIES GFORTRAN_LIBRARIES_DIR
-  VERSION_VAR GFORTRAN_VERSION_STRING)
-
-if(GFORTRANLIBS_FOUND)
-  message(STATUS "Looking for gfortran libraries -- found")
-  message(STATUS "gfortran version: ${GFORTRAN_VERSION_STRING}")
-else()
-  message(STATUS "Looking for gfortran libraries -- not found")
-endif()
-
-mark_as_advanced(LIBGFORTRAN_LIBRARIES LIBQUADMATH_LIBRARIES
-  LIBGOMP_LIBRARIES LIBGOMP_INCLUDE_DIR
-  GFORTRAN_LIBRARIES_DIR GFORTRAN_INCLUDE_DIR)
-# FindGFortranLIBS.cmake ends here
+# if (NOT "$ENV{FC}" STREQUAL "")
+# set(FC "$ENV{FC}" CACHE INTERNAL "Got from environment variable")
+# set(CMAKE_Fortran_COMPILER "$ENV{FC}" CACHE INTERNAL "Got from environment variable")
+# set(GFORTRANLIBS_FOUND "1" CACHE INTERNAL "found environment variable")
+# endif()
diff --git a/CMakeModules/FindWcsLib.cmake b/CMakeModules/FindWcsLib.cmake
index f024a1b..23d7e6f 100644
--- a/CMakeModules/FindWcsLib.cmake
+++ b/CMakeModules/FindWcsLib.cmake
@@ -32,20 +32,28 @@
 #  WCSLIB_LIBRARIES    - the WCSLIB libraries
 #                        (identical to WCSLIB_LIBRARY)
 
+# find paths
 if(NOT WCSLIB_FOUND)
 
-  find_path(WCSLIB_INCLUDE_DIR wcslib/wcs.h
-    HINTS ${WCSLIB_ROOT_DIR} PATH_SUFFIXES include)
-  find_library(WCSLIB_LIBRARY wcs
-    HINTS ${WCSLIB_ROOT_DIR} PATH_SUFFIXES lib)
+  if (NOT "$ENV{WCSLIB_ROOT_DIR}" STREQUAL "")
+  set(WCSLIB_ROOT "$ENV{WCSLIB_ROOT_DIR}" CACHE INTERNAL "Got from environment variable")
+  endif()
+
+  find_path(WCSLIB_INCLUDE wcslib/wcs.h
+    HINTS ${WCSLIB_ROOT} PATH_SUFFIXES include)
+  find_library(WCSLIB_LIB wcs
+    HINTS ${WCSLIB_ROOT} PATH_SUFFIXES lib)
   find_library(M_LIBRARY m)
-  mark_as_advanced(WCSLIB_INCLUDE_DIR WCSLIB_LIBRARY M_LIBRARY)
+  mark_as_advanced(WCSLIB_INCLUDE WCSLIB_LIB M_LIBRARY)
 
   include(FindPackageHandleStandardArgs)
   find_package_handle_standard_args(WCSLIB DEFAULT_MSG
-    WCSLIB_LIBRARY M_LIBRARY WCSLIB_INCLUDE_DIR)
+    WCSLIB_LIB M_LIBRARY WCSLIB_INCLUDE)
 
-  set(WCSLIB_INCLUDE_DIRS ${WCSLIB_INCLUDE_DIR})
-  set(WCSLIB_LIBRARIES ${WCSLIB_LIBRARY} ${M_LIBRARY})
+  set(WCSLIB_INCLUDE_DIRS ${WCSLIB_INCLUDE})
+  set(WCSLIB_LIBRARIES ${WCSLIB_LIB} ${M_LIBRARY})
+
+  set(WCSLIB_INCLUDE ${WCSLIB_INCLUDE})
+  set(WCSLIB_LIB ${WCSLIB_LIB} ${M_LIBRARY})
 
 endif(NOT WCSLIB_FOUND)
diff --git a/Docker/build_images.sh b/Docker/build_images.sh
new file mode 100755
index 0000000..513ddfa
--- /dev/null
+++ b/Docker/build_images.sh
@@ -0,0 +1,36 @@
+#!/usr/bin/env bash
+
+version=$(date +"%d%m%Y")
+
+for folder in \
+	sl7-test \
+	ubuntu1604-test
+
+do
+    name=${folder%%/}
+    echo
+    echo "==========================================="
+    echo "Building: "$name
+    echo "==========================================="
+
+    cmd="docker build -t fdiblen/$name:$version -t fdiblen/$name:latest ./$name"
+    echo
+    echo "Running:"
+    echo $cmd
+    echo
+    echo
+    eval $cmd
+
+    if [ $? -eq 0 ]
+    then
+      echo
+      echo "Successfully built the $name image!"
+    else
+      echo
+      echo "Could not build the $name image" >&2
+      exit $?
+    fi
+
+    #docker push fdiblen/$name:$version
+    #docker push fdiblen/$name:latest
+done
diff --git a/Docker/sl7-test/Dockerfile b/Docker/sl7-test/Dockerfile
new file mode 100644
index 0000000..11e4cfc
--- /dev/null
+++ b/Docker/sl7-test/Dockerfile
@@ -0,0 +1,42 @@
+FROM scientificlinux/sl:7
+MAINTAINER f.diblen@esciencecenter.nl
+
+# add EPEL repository for openblas
+RUN yum -y \
+    install https://dl.fedoraproject.org/pub/epel/epel-release-latest-7.noarch.rpm
+
+# install dependencies
+RUN yum -y install \
+    wget git pkgconfig make cmake3 cmake3-gui gcc-gfortran gcc-c++ flex bison \
+    openblas openblas-devel glib2-devel lapack lapack-devel cfitsio cfitsio-devel \
+    wcslib wcslib-devel ncurses ncurses-devel readline readline-devel \
+    python-devel boost boost-devel fftw fftw-devel hdf5 hdf5-devel \
+    numpy boost-python mpich mpich-devel fftw fftw-libs fftw-devel
+
+RUN mkdir /build && cd /build
+
+# compile casacore
+RUN git clone --progress --verbose https://github.com/casacore/casacore.git casacore_src && \
+    cd casacore_src && \
+    mkdir build && cd build && \
+    cmake3 .. -DUSE_FFTW3=ON \
+    -DCMAKE_INSTALL_PREFIX=/opt/casacore \
+    -DDATA_DIR=/opt/casacore/data -DUSE_OPENMP=ON \
+    -DUSE_HDF5=ON \
+    -DBUILD_PYTHON=ON \
+    -DUSE_THREADS=ON && \
+    make -j4 && \
+    make install
+
+
+## compile sagecal
+#RUN cd /build && \
+#    mkdir build-sl && cd build-sl && \
+#    cmake3 .. -DCMAKE_INSTALL_PREFIX=/opt/sagecal \
+#    -DCASACORE_ROOT_DIR=/opt/casacore \
+#    -DCASACORE_INCLUDE=/opt/casacore/include/casacore
+#    make -j4 && \
+#    make install && \
+#
+#RUN ls -alsrt /opt/sagecal && \
+#    /opt/sagecal/bin/sagecal
diff --git a/Docker/ubuntu1604-test/Dockerfile b/Docker/ubuntu1604-test/Dockerfile
new file mode 100644
index 0000000..1760ab1
--- /dev/null
+++ b/Docker/ubuntu1604-test/Dockerfile
@@ -0,0 +1,28 @@
+FROM ubuntu:xenial
+MAINTAINER f.diblen@esciencecenter.nl
+
+RUN apt-get update -y && \
+    apt-get install software-properties-common -y && \
+    add-apt-repository -s ppa:kernsuite/kern-3 -y && \
+    apt-add-repository multiverse && \
+    apt-get update -y
+
+RUN apt-get install -y \
+    git cmake g++ pkg-config \
+    libcfitsio-bin libcfitsio-dev \
+    libopenblas-base libopenblas-dev \
+    wcslib-dev wcslib-tools \
+    libglib2.0-dev \
+    libcasa-casa2 casacore-dev casacore-data casacore-tools \
+    mpich \
+    fftw-dev libfftw3-mpi3 libfftw3-bin
+
+
+## compile sagecal
+#RUN mkdir /build && cd /build \
+#    mkdir build-ubuntu && cd build-ubuntu && \
+#    cmake .. -DCMAKE_INSTALL_PREFIX=/opt/sagecal && \
+#    make -j4 && \
+#    make install && \
+#    ls -alsrt /opt/sagecal && \
+#    /opt/sagecal/bin/sagecal
diff --git a/build-tests/arch/compile-arch.sh b/build-tests/arch/compile-arch.sh
deleted file mode 100755
index 3d1260c..0000000
--- a/build-tests/arch/compile-arch.sh
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/bin/bash
-
-echo "Building SageCal" && \
-echo "Branch --> $BRANCH" && \
-cd /travis/workdir && \
-mkdir build-arch && cd build-arch && \
-cmake .. -DENABLE_CUDA=OFF && \
-make -j4 && \
-ls -alsrt ./dist/bin && \
-./dist/bin/sagecal
-
diff --git a/build-tests/arch/prepare-arch.sh b/build-tests/arch/prepare-arch.sh
deleted file mode 100755
index cbd94ad..0000000
--- a/build-tests/arch/prepare-arch.sh
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/bin/bash
-
-pacman -Syy --needed --noconfirm \
-	base-devel cmake git glib2 cfitsio wcslib
diff --git a/build-tests/compile_sagecal.sh b/build-tests/compile_sagecal.sh
new file mode 100755
index 0000000..7b639a4
--- /dev/null
+++ b/build-tests/compile_sagecal.sh
@@ -0,0 +1,56 @@
+#!/usr/bin/env bash
+
+echo 'script: ' $0
+
+echo "Building SageCal" && \
+echo "Branch --> $BRANCH" && \
+echo "Image --> $IMAGE"
+
+BUILD_DIR=$IMAGE'-build'
+
+cd /travis/workdir && \
+    mkdir $BUILD_DIR && cd $BUILD_DIR
+
+CMAKE_EXE=''
+OPTS=''
+
+case $IMAGE in
+    ubuntu)
+        echo 'Building for Ubuntu'
+        CMAKE_EXE=$(which cmake)
+        OPTS=''
+        ;;
+    sl7)
+        echo 'Building for Scientific Linux'
+        CMAKE_EXE=$(which cmake3)
+        OPTS='-DCASACORE_ROOT_DIR=/opt/casacore'
+        ;;
+    arch)
+        OPTS=''
+        ;;
+    *)
+        echo 'Unknown image $IMAGE!'
+        exit 1
+        ;;
+esac
+
+echo 'CMAKE_EXE: ' $CMAKE_EXE
+echo 'CMake options: ' $OPTS
+echo 'pwd: ' $PWD
+
+echo 'ls -asl: '
+ls -asl
+
+echo 'ls -asl /travis/workdir: '
+ls -asl /travis/workdir
+
+echo 'ls -asl /travis/workdir/$BUILD_DIR: '
+ls -asl /travis/workdir/$BUILD_DIR
+
+
+$CMAKE_EXE /travis/workdir -DCMAKE_INSTALL_PREFIX=/opt/sagecal $OPTS
+
+make -j4 && \
+make install && \
+ls -alsrt /opt/sagecal && \
+/opt/sagecal/bin/sagecal
diff --git a/build-tests/sl/compile-sl.sh b/build-tests/sl/compile-sl.sh
deleted file mode 100755
index 448c0fe..0000000
--- a/build-tests/sl/compile-sl.sh
+++ /dev/null
@@ -1,32 +0,0 @@
-#!/bin/bash
-
-echo "Building SageCal for Scientific Linux" && \
-echo "Branch --> $BRANCH" && \
-cd /travis/workdir
-
-# compile casacore first
-
-# mkdir -p /opt/soft/casacore/data
-# cd /opt/soft/casacore/data
-# wget -c ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar
-# tar zxfv WSRT_Measures.ztar && rm -f WSRT_Measures.ztar
-
-cd /travis/workdir
-git clone --progress --verbose https://github.com/casacore/casacore.git casacore_src && cd casacore_src
-
-mkdir build && cd build
-cmake3 -DUSE_FFTW3=ON -DCMAKE_INSTALL_PREFIX=/opt/soft/casacore -DDATA_DIR=/opt/soft/casacore/data -DUSE_OPENMP=ON \
-    -DUSE_HDF5=ON -DBUILD_PYTHON=ON -DUSE_THREADS=ON ..
-make -j4
-make install
-
-# compile sagecal
-cd /travis/workdir && \
-mkdir build-sl && cd build-sl
-
-cmake3 .. -DENABLE_CUDA=OFF && \
-    -DCASACORE_ROOT_DIR=/opt/soft/casacore -DCASACORE_INCLUDE=/opt/soft/casacore/include/casacore
-make -j4 && \
-ls -alsrt ./dist/bin && \
-./dist/bin/sagecal
-
diff --git a/build-tests/sl/prepare-sl.sh b/build-tests/sl/prepare-sl.sh
deleted file mode 100755
index 7b26817..0000000
--- a/build-tests/sl/prepare-sl.sh
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/bin/bash
-
-# add EPEL repository for openblas
-yum -y install https://dl.fedoraproject.org/pub/epel/epel-release-latest-7.noarch.rpm
-
-# install dependencies
-yum -y install wget git pkgconfig make cmake3 cmake3-gui gcc-gfortran gcc-c++ flex bison \
-       openblas openblas-devel glib2-devel lapack lapack-devel cfitsio cfitsio-devel \
-       wcslib wcslib-devel ncurses ncurses-devel readline readline-devel\
-       python-devel boost boost-devel fftw fftw-devel hdf5 hdf5-devel\
-       numpy boost-python
\ No newline at end of file
diff --git a/build-tests/ubuntu/compile-ubuntu.sh b/build-tests/ubuntu/compile-ubuntu.sh
deleted file mode 100755
index 061a81a..0000000
--- a/build-tests/ubuntu/compile-ubuntu.sh
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/bin/bash
-
-echo "Building SageCal" && \
-echo "Branch --> $BRANCH" && \
-cd /travis/workdir && \
-mkdir build-ubuntu && cd build-ubuntu && \
-cmake .. -DENABLE_CUDA=OFF && \
-make -j4 && \
-ls -alsrt ./dist/bin && \
-./dist/bin/sagecal
-
diff --git a/build-tests/ubuntu/prepare-ubuntu.sh b/build-tests/ubuntu/prepare-ubuntu.sh
deleted file mode 100755
index 79cf219..0000000
--- a/build-tests/ubuntu/prepare-ubuntu.sh
+++ /dev/null
@@ -1,9 +0,0 @@
-#!/bin/bash
-
-apt-get update -y
-apt-get install software-properties-common -y
-add-apt-repository -s ppa:kernsuite/kern-3 -y
-apt-add-repository multiverse
-apt-get update -y
-apt-get install -y git cmake g++ pkg-config libcfitsio-bin libcfitsio-dev libopenblas-base libopenblas-dev wcslib-dev wcslib-tools libglib2.0-dev libcasa-casa2 casacore-dev casacore-data casacore-tools
-
diff --git a/src/MPI/CMakeLists.txt b/src/MPI/CMakeLists.txt
index 27b6fae..2b85a7b 100644
--- a/src/MPI/CMakeLists.txt
+++ b/src/MPI/CMakeLists.txt
@@ -18,7 +18,7 @@ add_executable(sagecal-mpi ${SRCFILES})
 
 target_link_libraries(sagecal-mpi
  ${CASACORE_LIBRARIES}
- ${CFITSIO_LIBRARIES}
+ ${CFITSIO_LIB}
  ${OpenBLAS_LIB}
  ${LAPACK_LIBRARIES}
  ${WCSLIB_LIBRARIES}
@@ -29,7 +29,6 @@ target_link_libraries(sagecal-mpi
  -lopenblas
  -lgfortran
  -lpthread
- -lcfitsio
  -lm
  -ldirac
  -lradio
diff --git a/src/MPI/data.cpp b/src/MPI/data.cpp
deleted file mode 120000
index 96bd930..0000000
--- a/src/MPI/data.cpp
+++ /dev/null
@@ -1 +0,0 @@
-../MS/data.cpp
\ No newline at end of file
diff --git a/src/MPI/data.cpp b/src/MPI/data.cpp
new file mode 100644
index 0000000..c606217
--- /dev/null
+++ b/src/MPI/data.cpp
@@ -0,0 +1,1166 @@
+/*
+ *
+ Copyright (C) 2006-2008 Sarod Yatawatta <sarod@users.sf.net>  
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ $Id$
+ */
+
+#include "data.h"
+#include "Radio.h"
+#include "Dirac.h"
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/UVWMachine.h>
+#include <casacore/casa/Quanta.h>
+#include <casacore/casa/Quanta/Quantum.h>
+
+/* speed of light */
+#ifndef CONST_C
+#define CONST_C 299792458.0
+#endif
+
+using namespace casacore;
+
+int Data::numChannels=1;
+unsigned long int Data::numRows;
+
+
+char *Data::TableName = NULL;
+float Data::min_uvcut=0.0f;
+float Data::max_uvcut=100e6;
+float Data::max_uvtaper=0.0f;
+String Data::DataField = "DATA";
+String Data::OutField = "CORRECTED_DATA";
+int Data::TileSize = 120;
+int Data::Nt= 6; 
+char *Data::SkyModel=NULL;
+char *Data::Clusters=NULL;
+int Data::format=0; /* old LSM */
+double Data::nulow=2.0;
+double Data::nuhigh=30.0;
+
+int Data::max_emiter=3;
+int Data::max_iter=2;
+int Data::max_lbfgs=10;
+int Data::lbfgs_m=7;
+int Data::gpu_threads=128;
+int Data::linsolv=1;
+int Data::randomize=1;
+int Data::whiten=0;
+int Data::DoSim=0;
+int Data::DoDiag=0;
+int Data::doChan=0; /* if 1, solve for each channel in multi channel data */
+int Data::doBeam=0; /* if >0, enable LOFAR beam model */
+int Data::phaseOnly=0; /* if >0, enable phase only correction */
+int Data::solver_mode=0;
+int Data::ccid=-99999;
+double Data::rho=1e-9;
+char *Data::solfile=NULL;
+char *Data::initsolfile=NULL;
+char *Data::ignorefile=NULL;
+char *Data::MSlist=NULL;
+char *Data::MSpattern=NULL;
+
+/* distributed sagecal parameters */
+int Data::Nadmm=1;
+int Data::Npoly=2;
+int Data::PolyType=2;
+double Data::admm_rho=5.0;
+char *Data::admm_rho_file=NULL;
+int Data::aadmm=0;
+
+/* no upper limit, solve for all timeslots */
+int Data::Nmaxtime=0;
+/* skip starting time slots if given */
+int Data::Nskip=0;
+int Data::verbose=0; /* no verbose output */
+int Data::mdl=0; /* no AIC/MDL calculation by default */
+int Data::GPUpredict=0; /* use CPU for model calculation, if GPU not specified */
+#ifdef HAVE_CUDA
+int Data::heapsize=GPU_HEAP_SIZE; /* heap size in GPU (MB) to be used in malloc() */
+#endif
+
+int Data::servermode=-1; /* by default, no client-server mode */
+char *Data::servername=NULL;
+char *Data::portnumber=NULL; 
+
+using namespace Data;
+
+void
+Data::readMSlist(char *fname, vector<string> *msnames) {
+  cout<<"Reading "<<Data::MSlist<<endl;
+     /* multiple MS */
+     ifstream infile(fname);
+    /* check if the file exists and readable */
+    if(!infile.good()) {
+     cout <<"File "<<Data::MSlist<<" does not exist."<<endl;
+     exit(1);
+    }
+     string buffer;
+     if (infile.is_open()) {
+      while(infile.good()) {
+       std::getline(infile,buffer);
+       if (buffer.length()>0) {
+       cout<<buffer<<endl;
+        msnames->push_back(buffer);
+       }
+      }
+     }
+}
+
+void 
+Data::readAuxData(const char *fname, Data::IOData *data) {
+
+    Table _t=Table(fname);
+    Table _ant = Table(_t.keywordSet().asTable("ANTENNA"));
+    ROScalarColumn<String> a1(_ant, "NAME");
+    data->N=a1.nrow();
+    data->Nbase=data->N*(data->N-1)/2;
+    cout <<"Stations: "<<data->N<<" Baselines: "<<data->Nbase<<endl;
+
+    ROScalarColumn<double> timeCol(_t, "INTERVAL"); 
+    data->deltat=timeCol.get(0);
+    data->totalt=(timeCol.nrow()+data->Nbase+data->N-1)/(data->Nbase+data->N);
+    cout<<"Integration Time: "<<data->deltat<<" s,"<<" Total timeslots: "<<data->totalt<<endl;
+
+    Table _field = Table(_t.keywordSet().asTable("FIELD"));
+    ROArrayColumn<double> ref_dir(_field, "REFERENCE_DIR");
+    Array<double> dir = ref_dir(0);
+    double *c = dir.data();
+    data->ra0=c[0];
+    data->dec0=c[1];
+    cout<<"Phase center ("<< c[0] << ", " << c[1] <<")"<<endl;
+
+    //obtain the chanel freq information
+    Table _freq = Table(_t.keywordSet().asTable("SPECTRAL_WINDOW"));
+    ROArrayColumn<double> chan_freq(_freq, "CHAN_FREQ"); 
+    data->Nchan=chan_freq.shape(0)[0];
+    data->Nms=1;
+   /* allocate memory */
+   try { 
+   data->u=new double[data->Nbase*data->tilesz];
+   data->v=new double[data->Nbase*data->tilesz];
+   data->w=new double[data->Nbase*data->tilesz];
+   data->x=new double[8*data->Nbase*data->tilesz];
+   data->xo=new double[8*data->Nbase*data->tilesz*data->Nchan];
+   data->freqs=new double[data->Nchan];
+   data->flag=new double[data->Nbase*data->tilesz];
+   data->NchanMS=new int[data->Nms];
+   } catch (const std::bad_alloc& e) {
+     cout<<"Allocating memory for data failed. Quitting."<< e.what() << endl;
+     exit(1);
+   }
+   data->NchanMS[0]=data->Nchan;
+
+   /* copy freq */
+   data->freq0=0.0;
+   for (int ci=0; ci<data->Nchan; ci++) {
+     data->freqs[ci]=chan_freq(0).data()[ci];
+     data->freq0+=data->freqs[ci];
+   }
+   data->freq0/=(double)data->Nchan;
+   /* need channel widths to calculate bandwidth */
+   ROArrayColumn<double> chan_width(_freq, "CHAN_WIDTH"); 
+   data->deltaf=(double)data->Nchan*(chan_width(0).data()[0]);
+}
+
+void 
+Data::readAuxData(const char *fname, Data::IOData *data, Data::LBeam *binfo) {
+
+    Table _t=Table(fname);
+    Table _ant = Table(_t.keywordSet().asTable("ANTENNA"));
+    ROScalarColumn<String> a1(_ant, "NAME");
+    data->N=a1.nrow();
+    data->Nbase=data->N*(data->N-1)/2;
+    cout <<"Stations: "<<data->N<<" Baselines: "<<data->Nbase<<endl;
+
+    ROScalarColumn<double> timeCol(_t, "INTERVAL"); 
+    data->deltat=timeCol.get(0);
+    data->totalt=(timeCol.nrow()+data->Nbase+data->N-1)/(data->Nbase+data->N);
+    cout<<"Integration Time: "<<data->deltat<<" s,"<<" Total timeslots: "<<data->totalt<<endl;
+
+    Table _field = Table(_t.keywordSet().asTable("FIELD"));
+    ROArrayColumn<double> ref_dir(_field, "PHASE_DIR"); /* old REFERENCE_DIR */
+    Array<double> dir = ref_dir(0);
+    double *c = dir.data();
+    data->ra0=c[0];
+    data->dec0=c[1];
+    cout<<"Phase center ("<< c[0] << ", " << c[1] <<")"<<endl;
+
+    //obtain the chanel freq information
+    Table _freq = Table(_t.keywordSet().asTable("SPECTRAL_WINDOW"));
+    ROArrayColumn<double> chan_freq(_freq, "CHAN_FREQ"); 
+    data->Nchan=chan_freq.shape(0)[0];
+    data->Nms=1;
+   /* allocate memory */
+   try {
+   data->u=new double[data->Nbase*data->tilesz];
+   data->v=new double[data->Nbase*data->tilesz];
+   data->w=new double[data->Nbase*data->tilesz];
+   data->x=new double[8*data->Nbase*data->tilesz];
+   data->xo=new double[8*data->Nbase*data->tilesz*data->Nchan];
+   data->freqs=new double[data->Nchan];
+   data->flag=new double[data->Nbase*data->tilesz];
+   data->NchanMS=new int[data->Nms];
+   } catch (const std::bad_alloc& e) {
+     cout<<"Allocating memory for data failed. Quitting."<< e.what() << endl;
+     exit(1);
+   }
+   data->NchanMS[0]=data->Nchan;
+
+   /* copy freq */
+   data->freq0=0.0;
+   for (int ci=0; ci<data->Nchan; ci++) {
+     data->freqs[ci]=chan_freq(0).data()[ci];
+     data->freq0+=data->freqs[ci];
+   }
+   data->freq0/=(double)data->Nchan;
+   /* need channel widths to calculate bandwidth */
+   ROArrayColumn<double> chan_width(_freq, "CHAN_WIDTH"); 
+   data->deltaf=(double)data->Nchan*(chan_width(0).data()[0]);
+
+   /* UTC time */
+   binfo->time_utc=new double[data->tilesz]; 
+   /* no of elements in each station */
+   binfo->Nelem=new int[data->N];
+   /* positions of stations */
+   binfo->sx=new double[data->N];
+   binfo->sy=new double[data->N];
+   binfo->sz=new double[data->N];
+   /* coordinates of elements */
+   binfo->xx=new double*[data->N];
+   binfo->yy=new double*[data->N];
+   binfo->zz=new double*[data->N];
+
+   Table antfield;
+   if(_t.keywordSet().fieldNumber("LOFAR_ANTENNA_FIELD") != -1) {
+    antfield = Table(_t.keywordSet().asTable("LOFAR_ANTENNA_FIELD"));
+   } else {
+    char buff[2048]={0};
+    sprintf(buff, "%s/LOFAR_ANTENNA_FIELD", fname);
+    antfield=Table(buff);
+   }
+   ROArrayColumn<double> position(antfield, "POSITION"); 
+   ROArrayColumn<double> offset(antfield, "ELEMENT_OFFSET"); 
+   ROArrayColumn<double> coord(antfield, "COORDINATE_AXES"); 
+   ROArrayColumn<bool> eflag(antfield, "ELEMENT_FLAG");
+   ROArrayColumn<double> tileoffset(antfield, "TILE_ELEMENT_OFFSET");
+   /* check if TILE_ELEMENT_OFFSET has any rows, of no rows present, 
+      we know this is LBA */
+   bool isHBA=tileoffset.hasContent(0);
+
+   /* read positions, also setup memory for element coords */
+   for (int ci=0; ci<data->N; ci++) {
+     Array<double> _pos=position(ci);
+     double *tx=_pos.data();
+     binfo->sz[ci]=tx[2];
+
+     MPosition stnpos(MVPosition(tx[0],tx[1],tx[2]),MPosition::ITRF);
+     Array<double> _radpos=stnpos.getAngle("rad").getValue();
+     tx=_radpos.data();
+
+     binfo->sx[ci]=tx[0];
+     binfo->sy[ci]=tx[1];
+//cout<<"position "<<binfo->sx[ci]<<" "<<binfo->sy[ci]<<" "<<binfo->sz[ci]<<endl;
+     binfo->Nelem[ci]=offset.shape(ci)[1];
+   }
+
+
+   if (isHBA) {
+    double cones[16];
+    for (int ci=0; ci<16; ci++) {
+      cones[ci]=1.0;
+    }
+    double tempT[3*16];
+    /* now read in element offsets, also transform them to local coordinates */
+    for (int ci=0; ci<data->N; ci++) {
+      Array<double> _off=offset(ci);
+      double *off=_off.data();
+      Array<double> _coord=coord(ci);
+      double *coordmat=_coord.data();
+      Array<bool> _eflag=eflag(ci);
+      bool *ef=_eflag.data();
+      Array<double> _toff=tileoffset(ci);
+      double *toff=_toff.data();
+
+/*
+cout<<"A=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<off[3*cj]<<","<<off[3*cj+1]<<","<<off[3*cj+2]<<endl;
+     }
+cout<<"];"<<endl;
+cout<<"T0=["<<endl;
+     for (int cj=0; cj<16; cj++) {
+cout<<toff[3*cj]<<","<<toff[3*cj+1]<<","<<toff[3*cj+2]<<endl;
+     }
+cout<<"];"<<endl;
+cout<<"BT=["<<endl;
+     for (int cj=0; cj<3; cj++) {
+cout<<coordmat[3*cj]<<","<<coordmat[3*cj+1]<<","<<coordmat[3*cj+2]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+
+      double *tempC=new double[3*binfo->Nelem[ci]];
+      my_dgemm('T', 'N', binfo->Nelem[ci], 3, 3, 1.0, off, 3, coordmat, 3, 0.0, tempC, binfo->Nelem[ci]); 
+      my_dgemm('T', 'N', 16, 3, 3, 1.0, toff, 3, coordmat, 3, 0.0, tempT, 16); 
+  
+/*
+cout<<"C=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<tempC[cj]<<","<<tempC[cj+binfo->Nelem[ci]]<<","<<tempC[cj+2*binfo->Nelem[ci]]<<endl;
+     }
+cout<<"];"<<endl;
+cout<<"T1=["<<endl;
+     for (int cj=0; cj<16; cj++) {
+cout<<tempT[cj]<<","<<tempT[cj+16]<<","<<tempT[cj+2*16]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+
+      /* now inspect the element flag table to see if any of the dipoles are flagged */
+      int fcount=0;
+      for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+       if (ef[2*cj]==1 || ef[2*cj+1]==1) {
+        fcount++;
+       } 
+      }
+
+//cout<<"%%Flagged "<<fcount<<endl;
+
+      binfo->xx[ci]=new double[16*(binfo->Nelem[ci]-fcount)];
+      binfo->yy[ci]=new double[16*(binfo->Nelem[ci]-fcount)];
+      binfo->zz[ci]=new double[16*(binfo->Nelem[ci]-fcount)];
+      /* copy unflagged coords, 16 times for each dipole */
+      fcount=0;
+      for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+       if (!(ef[2*cj]==1 || ef[2*cj+1]==1)) {
+        //binfo->xx[ci][fcount]=tempC[cj];
+        //binfo->yy[ci][fcount]=tempC[cj+binfo->Nelem[ci]];
+        //binfo->zz[ci][fcount]=tempC[cj+2*binfo->Nelem[ci]];
+        my_dcopy(16,&tempT[0],1,&(binfo->xx[ci][fcount]),1);
+        my_daxpy(16,cones,tempC[cj],&(binfo->xx[ci][fcount]));
+        my_dcopy(16,&tempT[16],1,&(binfo->yy[ci][fcount]),1);
+        my_daxpy(16,cones,tempC[cj+binfo->Nelem[ci]],&(binfo->yy[ci][fcount]));
+        my_dcopy(16,&tempT[24],1,&(binfo->zz[ci][fcount]),1);
+        my_daxpy(16,cones,tempC[cj+2*binfo->Nelem[ci]],&(binfo->zz[ci][fcount]));
+        fcount+=16;
+       }
+      }
+      binfo->Nelem[ci]=fcount;
+/*
+cout<<"%%Copied "<<fcount<<endl;
+cout<<"D=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<binfo->xx[ci][cj]<<","<<binfo->yy[ci][cj]<<","<<binfo->zz[ci][cj]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+      delete [] tempC;
+    }
+   } else { /* LBA */
+    /* now read in element offsets, also transform them to local coordinates */
+    for (int ci=0; ci<data->N; ci++) {
+      Array<double> _off=offset(ci);
+      double *off=_off.data();
+      Array<double> _coord=coord(ci);
+      double *coordmat=_coord.data();
+      Array<bool> _eflag=eflag(ci);
+      bool *ef=_eflag.data();
+
+/*
+cout<<"A=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<off[3*cj]<<","<<off[3*cj+1]<<","<<off[3*cj+2]<<endl;
+     }
+cout<<"];"<<endl;
+cout<<"BT=["<<endl;
+     for (int cj=0; cj<3; cj++) {
+cout<<coordmat[3*cj]<<","<<coordmat[3*cj+1]<<","<<coordmat[3*cj+2]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+
+      double *tempC=new double[3*binfo->Nelem[ci]];
+      my_dgemm('T', 'N', binfo->Nelem[ci], 3, 3, 1.0, off, 3, coordmat, 3, 0.0, tempC, binfo->Nelem[ci]); 
+  
+/*
+cout<<"C=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<tempC[cj]<<","<<tempC[cj+binfo->Nelem[ci]]<<","<<tempC[cj+2*binfo->Nelem[ci]]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+
+      /* now inspect the element flag table to see if any of the dipoles are flagged */
+      int fcount=0;
+      for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+       if (ef[2*cj]==1 || ef[2*cj+1]==1) {
+        fcount++;
+       } 
+      }
+
+//cout<<"%%Flagged "<<fcount<<endl;
+
+      binfo->xx[ci]=new double[(binfo->Nelem[ci]-fcount)];
+      binfo->yy[ci]=new double[(binfo->Nelem[ci]-fcount)];
+      binfo->zz[ci]=new double[(binfo->Nelem[ci]-fcount)];
+      /* copy unflagged coords for each dipole */
+      fcount=0;
+      for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+       if (!(ef[2*cj]==1 || ef[2*cj+1]==1)) {
+        binfo->xx[ci][fcount]=tempC[cj];
+        binfo->yy[ci][fcount]=tempC[cj+binfo->Nelem[ci]];
+        binfo->zz[ci][fcount]=tempC[cj+2*binfo->Nelem[ci]];
+        fcount++;
+       }
+      }
+      binfo->Nelem[ci]=fcount;
+/*
+cout<<"%%Copied "<<fcount<<endl;
+cout<<"D=["<<endl;
+     for (int cj=0; cj<binfo->Nelem[ci]; cj++) {
+cout<<binfo->xx[ci][cj]<<","<<binfo->yy[ci][cj]<<","<<binfo->zz[ci][cj]<<endl;
+     }
+cout<<"];"<<endl;
+*/
+      delete [] tempC;
+    }
+   }
+
+   /* read beam pointing direction (use Tile beam info: LBA also has it) */
+   ROArrayColumn<double> point_dir(_field, "LOFAR_TILE_BEAM_DIR");
+   Array<double> pdir = point_dir(0);
+   double *pc = pdir.data();
+   binfo->p_ra0=pc[0];
+   binfo->p_dec0=pc[1];
+
+   /* convert positions from xyz to longitude,latitude,height */
+/*   double *longitude=new double[data->N];
+   double *latitude=new double[data->N];
+   double *height=new double[data->N];
+   xyz2llh(binfo->sx,binfo->sy,binfo->sz,longitude,latitude,height,data->N);
+   delete [] binfo->sx;
+   delete [] binfo->sy;
+   delete [] binfo->sz;
+   binfo->sx=longitude;
+   binfo->sy=latitude;
+   binfo->sz=height;
+*/
+}
+
+
+void 
+Data::readAuxDataList(vector<string> msnames, Data::IOData *data) {
+    /* read first filename */
+    const char *fname=msnames[0].c_str();
+    Table _t=Table(fname);
+    //char buff[2048] = {0};
+    //sprintf(buff, "%s/ANTENNA", fname);
+    //Table _ant=Table(buff);
+    Table _ant = Table(_t.keywordSet().asTable("ANTENNA"));
+    ROScalarColumn<String> a1(_ant, "NAME");
+    data->N=a1.nrow();
+    data->Nbase=data->N*(data->N-1)/2;
+    cout <<"Stations: "<<data->N<<" Baselines: "<<data->Nbase<<endl;
+
+    ROScalarColumn<double> timeCol(_t, "INTERVAL"); 
+    data->deltat=timeCol.get(0);
+    data->totalt=(timeCol.nrow()+data->Nbase+data->N-1)/(data->Nbase+data->N);
+    cout<<"Integration Time: "<<data->deltat<<" s,"<<" Total timeslots: "<<data->totalt<<endl;
+
+    //sprintf(buff, "%s/FIELD", fname);
+    //Table _field = Table(buff);
+    Table _field = Table(_t.keywordSet().asTable("FIELD"));
+    ROArrayColumn<double> ref_dir(_field, "REFERENCE_DIR");
+    Array<double> dir = ref_dir(0);
+    double *c = dir.data();
+    data->ra0=c[0];
+    data->dec0=c[1];
+    cout<<"Phase center ("<< c[0] << ", " << c[1] <<")"<<endl;
+
+    data->Nchan=0;
+    data->Nms=msnames.size();
+    data->NchanMS=new int[data->Nms];
+    for (int cm=0; cm<data->Nms; cm++) {
+     //obtain the chanel freq information
+     fname=msnames[cm].c_str();
+     Table _t1=Table(fname);
+     //sprintf(buff, "%s/SPECTRAL_WINDOW", fname);
+     //Table _freq = Table(buff);
+     Table _freq = Table(_t1.keywordSet().asTable("SPECTRAL_WINDOW"));
+     ROArrayColumn<double> chan_freq(_freq, "CHAN_FREQ"); 
+     data->Nchan+=chan_freq.shape(0)[0];
+     data->NchanMS[cm]=chan_freq.shape(0)[0];
+    }
+   cout<<"Total channels="<<data->Nchan<<endl;
+ 
+   /* allocate memory */
+   data->u=new double[data->Nbase*data->tilesz];
+   data->v=new double[data->Nbase*data->tilesz];
+   data->w=new double[data->Nbase*data->tilesz];
+   data->x=new double[8*data->Nbase*data->tilesz];
+   data->xo=new double[8*data->Nbase*data->tilesz*data->Nchan];
+   data->freqs=new double[data->Nchan];
+   data->flag=new double[data->Nbase*data->tilesz];
+
+   /* copy freq */
+   data->freq0=0.0;
+   data->deltaf=0.0;
+   int ck=0;
+   for (int cm=0; cm<data->Nms; cm++) {
+     fname=msnames[cm].c_str();
+     Table _t1=Table(fname);
+     //sprintf(buff, "%s/SPECTRAL_WINDOW", fname);
+     //Table _freq = Table(buff);
+     Table _freq = Table(_t1.keywordSet().asTable("SPECTRAL_WINDOW"));
+     ROArrayColumn<double> chan_freq(_freq, "CHAN_FREQ"); 
+     /* need channel widths to calculate bandwidth */
+     ROArrayColumn<double> chan_width(_freq, "CHAN_WIDTH"); 
+     for (int ci=0; ci<chan_freq.shape(0)[0]; ci++) {
+     data->freqs[ck]=chan_freq(0).data()[ci];
+     data->freq0+=data->freqs[ck++];
+     data->deltaf+=(chan_width(0).data()[ci]);
+    }
+   }
+   data->freq0/=(double)data->Nchan;
+   cout<<"freq0="<<data->freq0/1e6<<endl;
+   cout<<"deltaf="<<data->deltaf/1e6<<endl;
+   for (ck=0; ck<data->Nchan;ck++){
+    cout<<ck<<" "<<data->freqs[ck]<<endl;
+   }
+}
+
+
+/* each time this is called read in data from MS, and format them as
+  u,v,w: u,v,w coordinates (wavelengths) size Nbase*tilesz x 1 
+  u,v,w are ordered with baselines, timeslots
+  x: data to write size Nbase*8*tilesz x 1
+  ordered by XX(re,im),XY(re,im),YX(re,im), YY(re,im), baseline, timeslots
+  fratio: flagged data as a ratio to all available data
+*/
+void 
+Data::loadData(Table ti, Data::IOData iodata, double *fratio) {
+
+    /* sort input table by ant1 and ant2 */
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=ti.sort(iv1,Sort::Ascending);
+
+    ROScalarColumn<int> a1(t, "ANTENNA1"), a2(t, "ANTENNA2");
+    /* only read only access for input */
+    ROArrayColumn<Complex> dataCol(t, Data::DataField);
+    ROArrayColumn<double> uvwCol(t, "UVW"); 
+    ROArrayColumn<bool> flagCol(t, "FLAG");
+
+    /* check we get correct rows */
+    int nrow=t.nrow();
+    if(nrow<iodata.tilesz*iodata.Nbase-iodata.tilesz*iodata.N) {
+      cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
+    }
+    int row0=0;
+    /* tapering */
+    bool dotaper=false;
+    double invtaper=1.0;
+    if (max_uvtaper>0.0f) {
+      dotaper=true;
+      /* taper in m */
+      invtaper=iodata.freq0/((double)max_uvtaper*CONST_C);
+    }
+    /* counters for finding flagged data ratio */
+    int countgood=0; int countbad=0;
+    for(int row = 0; row < nrow && row0<iodata.tilesz*iodata.Nbase; row++) {
+        uInt i = a1(row); //antenna1 
+        uInt j = a2(row); //antenna2
+        /* only work with cross correlations */
+        if (i!=j) {
+        Array<Complex> data = dataCol(row);
+        Matrix<double> uvw = uvwCol(row);
+        Array<bool> flag = flagCol(row);
+
+        Complex cxx(0, 0);
+        Complex cxy(0, 0);
+        Complex cyx(0, 0);
+        Complex cyy(0, 0);
+        /* calculate sqrt(u^2+v^2) to select uv cuts */
+        double *c = uvw.data();
+        double uvd=sqrt(c[0]*c[0]+c[1]*c[1]);
+        bool flag_uvcut=0;
+        if (uvd<min_uvcut || uvd>max_uvcut) {
+          flag_uvcut=true;
+        } 
+        double uvtaper=1.0;
+        if (dotaper) {
+         uvtaper=uvd*invtaper;
+         if (uvtaper>1.0) {
+          uvtaper=1.0;
+         }
+        }
+        int nflag=0;
+        for(int k = 0; k < iodata.Nchan; k++) {
+           Complex *ptr = data[k].data();
+           bool *flgptr=flag[k].data();
+           if (!flgptr[0] && !flgptr[1] && !flgptr[2] && !flgptr[3]){
+             cxx+=ptr[0];
+             cxy+=ptr[1];
+             cyx+=ptr[2];
+             cyy+=ptr[3];
+             nflag++; /* remeber unflagged datapoints */ 
+           } 
+        
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8]=ptr[0].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+1]=ptr[0].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+2]=ptr[1].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+3]=ptr[1].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+4]=ptr[2].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+5]=ptr[2].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+6]=ptr[3].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+7]=ptr[3].imag();
+        }
+        if (nflag>iodata.Nchan/2) { /* at least half channels should have good data */
+         double invnflag=1.0/(double)nflag;
+         cxx*=invnflag;
+         cxy*=invnflag;
+         cyx*=invnflag;
+         cyy*=invnflag;
+         if (dotaper) {
+          cxx*=uvtaper;
+          cxy*=uvtaper;
+          cyx*=uvtaper;
+          cyy*=uvtaper;
+         }
+         iodata.flag[row0]=0;
+         countgood++;
+        } else {
+         if (!nflag) {
+         /* all channels flagged, flag this row */
+          iodata.flag[row0]=1;
+          countbad++;
+         } else {
+          iodata.flag[row0]=2;
+         }
+        }
+        iodata.u[row0]=c[0];
+        iodata.v[row0]=c[1];
+        iodata.w[row0]=c[2];
+        if (flag_uvcut) {
+            iodata.flag[row0]=2;
+        }
+        iodata.x[row0*8]=cxx.real();
+        iodata.x[row0*8+1]=cxx.imag();
+        iodata.x[row0*8+2]=cxy.real();
+        iodata.x[row0*8+3]=cxy.imag();
+        iodata.x[row0*8+4]=cyx.real();
+        iodata.x[row0*8+5]=cyx.imag();
+        iodata.x[row0*8+6]=cyy.real();
+        iodata.x[row0*8+7]=cyy.imag();
+
+       row0++;
+      }
+    }
+    /* now if there is a tail of empty data remaining, flag them */
+    if (row0<iodata.tilesz*iodata.Nbase) {
+      for(int row = row0; row<iodata.tilesz*iodata.Nbase; row++) {
+        iodata.flag[row]=1;
+
+      }
+      /* set uvw and data to 0 to eliminate any funny business */
+      memset(&iodata.u[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.v[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.w[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.x[8*row0],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
+
+      for(int k = 0; k < iodata.Nchan; k++) {
+       memset(&iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
+      }
+    }
+    /* flagged data / total usable data, not counting excluded baselines */
+    if (countgood+countbad>0) {
+     *fratio=(double)countbad/(double)(countgood+countbad);
+    } else {
+     *fratio=1.0;
+    }
+}
+
+/* each time this is called read in data from MS, and format them as
+  u,v,w: u,v,w coordinates (wavelengths) size Nbase*tilesz x 1 
+  u,v,w are ordered with baselines, timeslots
+  x: data to write size Nbase*8*tileze x 1
+  ordered by XX(re,im),XY(re,im),YX(re,im), YY(re,im), baseline, timeslots
+
+  time_utc: UTC time, 1 per each Nbase
+  fratio: flagged data as a ratio to all available data
+*/
+void 
+Data::loadData(Table ti, Data::IOData iodata, LBeam binfo, double *fratio) {
+
+    /* sort input table by ant1 and ant2 */
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=ti.sort(iv1,Sort::Ascending);
+
+    ROScalarColumn<int> a1(t, "ANTENNA1"), a2(t, "ANTENNA2");
+    /* only read only access for input */
+    ROArrayColumn<Complex> dataCol(t, Data::DataField);
+    ROArrayColumn<double> uvwCol(t, "UVW"); 
+    ROArrayColumn<bool> flagCol(t, "FLAG");
+    ROScalarColumn<double> tut(t,"TIME");
+
+    /* check we get correct rows */
+    int nrow=t.nrow();
+    if(nrow<iodata.tilesz*iodata.Nbase-iodata.tilesz*iodata.N) {
+      cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
+    }
+    int row0=0;
+    int rowt=0;
+    /* tapering */
+    bool dotaper=false;
+    double invtaper=1.0;
+    if (max_uvtaper>0.0f) {
+      dotaper=true;
+      /* taper in m */
+      invtaper=iodata.freq0/((double)max_uvtaper*CONST_C);
+    }
+    /* counters for finding flagged data ratio */
+    int countgood=0; int countbad=0;
+    for(int row = 0; row < nrow && row0<iodata.tilesz*iodata.Nbase; row++) {
+        uInt i = a1(row); //antenna1 
+        uInt j = a2(row); //antenna2
+        if (!i && !j) {/* use baseline 0-0 to extract time */
+         double tt=tut(row);
+         /* convert MJD (s) to JD (days) */
+         binfo.time_utc[rowt++]=(tt/86400.0+2400000.5); /* no +0.5 added */
+        }
+        /* only work with cross correlations */
+        if (i!=j) {
+        Array<Complex> data = dataCol(row);
+        Matrix<double> uvw = uvwCol(row);
+        Array<bool> flag = flagCol(row);
+
+        Complex cxx(0, 0);
+        Complex cxy(0, 0);
+        Complex cyx(0, 0);
+        Complex cyy(0, 0);
+        /* calculate sqrt(u^2+v^2) to select uv cuts */
+        double *c = uvw.data();
+        double uvd=sqrt(c[0]*c[0]+c[1]*c[1]);
+        bool flag_uvcut=0;
+        if (uvd<min_uvcut || uvd>max_uvcut) {
+          flag_uvcut=true;
+        } 
+        double uvtaper=1.0;
+        if (dotaper) {
+         uvtaper=uvd*invtaper;
+         if (uvtaper>1.0) {
+          uvtaper=1.0;
+         }
+        }
+        int nflag=0;
+        for(int k = 0; k < iodata.Nchan; k++) {
+           Complex *ptr = data[k].data();
+           bool *flgptr=flag[k].data();
+           if (!flgptr[0] && !flgptr[1] && !flgptr[2] && !flgptr[3]){
+             cxx+=ptr[0];
+             cxy+=ptr[1];
+             cyx+=ptr[2];
+             cyy+=ptr[3];
+             nflag++; /* remeber unflagged datapoints */ 
+           } 
+        
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8]=ptr[0].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+1]=ptr[0].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+2]=ptr[1].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+3]=ptr[1].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+4]=ptr[2].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+5]=ptr[2].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+6]=ptr[3].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+7]=ptr[3].imag();
+        }
+        if (nflag>iodata.Nchan/2) { /* at least half channels should have good data */
+         double invnflag=1.0/(double)nflag;
+         cxx*=invnflag;
+         cxy*=invnflag;
+         cyx*=invnflag;
+         cyy*=invnflag;
+         if (dotaper) {
+          cxx*=uvtaper;
+          cxy*=uvtaper;
+          cyx*=uvtaper;
+          cyy*=uvtaper;
+         }
+         iodata.flag[row0]=0;
+         countgood++;
+        } else {
+         if (!nflag) {
+         /* all channels flagged, flag this row */
+          iodata.flag[row0]=1;
+          countbad++;
+         } else {
+          iodata.flag[row0]=2;
+         }
+        }
+        iodata.u[row0]=c[0];
+        iodata.v[row0]=c[1];
+        iodata.w[row0]=c[2];
+        if (flag_uvcut) {
+            iodata.flag[row0]=2;
+        }
+        iodata.x[row0*8]=cxx.real();
+        iodata.x[row0*8+1]=cxx.imag();
+        iodata.x[row0*8+2]=cxy.real();
+        iodata.x[row0*8+3]=cxy.imag();
+        iodata.x[row0*8+4]=cyx.real();
+        iodata.x[row0*8+5]=cyx.imag();
+        iodata.x[row0*8+6]=cyy.real();
+        iodata.x[row0*8+7]=cyy.imag();
+
+       row0++;
+      }
+    }
+    /* now if there is a tail of empty data remaining, flag them */
+    if (row0<iodata.tilesz*iodata.Nbase) {
+      for(int row = row0; row<iodata.tilesz*iodata.Nbase; row++) {
+        iodata.flag[row]=1;
+
+      }
+      /* set uvw and data to 0 to eliminate any funny business */
+      memset(&iodata.u[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.v[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.w[row0],0,sizeof(double)*(size_t)(iodata.tilesz*iodata.Nbase-row0));
+      memset(&iodata.x[8*row0],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
+
+      for(int k = 0; k < iodata.Nchan; k++) {
+       memset(&iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8],0,sizeof(double)*(size_t)8*(iodata.tilesz*iodata.Nbase-row0));
+      }
+    }
+    /* flagged data / total usable data, not counting excluded baselines */
+    if (countgood+countbad>0) {
+     *fratio=(double)countbad/(double)(countgood+countbad);
+    } else {
+     *fratio=1.0;
+    }
+}
+
+
+/* each time this is called read in data from MS, and format them as
+  u,v,w: u,v,w coordinates (wavelengths) size Nbase*tilesz x 1 
+  u,v,w are ordered with baselines, timeslots
+  x: data to write size Nbase*8*tileze x 1
+  ordered by XX(re,im),XY(re,im),YX(re,im), YY(re,im), baseline, timeslots
+  fratio: flagged data as a ratio to all available data
+*/
+void 
+Data::loadDataList(vector<MSIter*> msitr, Data::IOData iodata, double *fratio) {
+    Table ti=msitr[0]->table();
+    /* sort input table by ant1 and ant2 */
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=ti.sort(iv1,Sort::Ascending);
+
+    ROScalarColumn<int> a1(t, "ANTENNA1"), a2(t, "ANTENNA2");
+    /* only read only access for input */
+    ROArrayColumn<double> uvwCol(t, "UVW"); 
+
+    /* check we get correct rows */
+    int nrow=t.nrow();
+    if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
+      cout<<"Error in rows"<<endl;
+    }
+  vector<ROArrayColumn<Complex>* > dataCols(iodata.Nms);
+  vector<ROArrayColumn<bool>* > flagCols(iodata.Nms);
+  for (int cm=0; cm<iodata.Nms;cm++) { 
+    Table tti=(msitr[cm]->table());
+    Table *tt=new Table(tti.sort(iv1,Sort::Ascending));
+    dataCols[cm] = new  ROArrayColumn<Complex>(*tt,Data::DataField);
+    flagCols[cm] = new  ROArrayColumn<bool>(*tt,"FLAG");
+  }
+    /* tapering */
+    bool dotaper=false;
+    double invtaper=1.0;
+    if (max_uvtaper>0.0f) {
+      dotaper=true;
+      /* taper in m */
+      invtaper=iodata.freq0/((double)max_uvtaper*CONST_C);
+    }
+
+    /* counters for finding flagged data ratio */
+    int countgood=0; int countbad=0;
+
+    int row0=0;
+    for(int row = 0; row < nrow; row++) {
+        uInt i = a1(row); //antenna1 
+        uInt j = a2(row); //antenna2
+        /* only work with cross correlations */
+        if (i!=j) {
+        Matrix<double> uvw = uvwCol(row);
+
+        Complex cxx(0, 0);
+        Complex cxy(0, 0);
+        Complex cyx(0, 0);
+        Complex cyy(0, 0);
+        /* calculate sqrt(u^2+v^2) to select uv cuts */
+        double *c = uvw.data();
+        double uvd=sqrt(c[0]*c[0]+c[1]*c[1]);
+        bool flag_uvcut=0;
+        if (uvd<min_uvcut || uvd>max_uvcut) {
+          flag_uvcut=true;
+        } 
+        int nflag=0;
+        double uvtaper=1.0;
+        if (dotaper) {
+         uvtaper=uvd*invtaper;
+         if (uvtaper>1.0) {
+          uvtaper=1.0;
+         }
+        }
+
+  int chanoff=0;
+  for (int cm=0; cm<iodata.Nms;cm++) { 
+        Array<Complex> data = (*(dataCols[cm]))(row);
+        Array<bool> flag = (*(flagCols[cm]))(row);
+        for(int k = 0; k < iodata.NchanMS[cm]; k++) {
+           Complex *ptr = data[k].data();
+           bool *flgptr=flag[k].data();
+           if (!flgptr[0] && !flgptr[1] && !flgptr[2] && !flgptr[3]){
+             cxx+=ptr[0];
+             cxy+=ptr[1];
+             cyx+=ptr[2];
+             cyy+=ptr[3];
+             nflag++; /* remeber unflagged datapoints */ 
+           } 
+        
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8]=ptr[0].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+1]=ptr[0].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+2]=ptr[1].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+3]=ptr[1].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+4]=ptr[2].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+5]=ptr[2].imag();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+6]=ptr[3].real();
+           iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+7]=ptr[3].imag();
+           chanoff++;
+        }
+   }
+
+        if (nflag>iodata.Nchan/2) {/* at least half channels should have good data */
+         double invnflag=1.0/(double)nflag;
+         cxx*=invnflag;
+         cxy*=invnflag;
+         cyx*=invnflag;
+         cyy*=invnflag;
+         if (dotaper) {
+          cxx*=uvtaper;
+          cxy*=uvtaper;
+          cyx*=uvtaper;
+          cyy*=uvtaper;
+         }
+         iodata.flag[row0]=0;
+         countgood++;
+        } else {
+         if (!nflag) {
+         /* all channels flagged, flag this row */
+          iodata.flag[row0]=1;
+          countbad++;
+         } else {
+          iodata.flag[row0]=2;
+         }
+        }
+        iodata.u[row0]=c[0];
+        iodata.v[row0]=c[1];
+        iodata.w[row0]=c[2];
+        if (flag_uvcut) {
+            iodata.flag[row0]=2;
+        }
+        iodata.x[row0*8]=cxx.real();
+        iodata.x[row0*8+1]=cxx.imag();
+        iodata.x[row0*8+2]=cxy.real();
+        iodata.x[row0*8+3]=cxy.imag();
+        iodata.x[row0*8+4]=cyx.real();
+        iodata.x[row0*8+5]=cyx.imag();
+        iodata.x[row0*8+6]=cyy.real();
+        iodata.x[row0*8+7]=cyy.imag();
+
+       row0++;
+      }
+    }
+    /* now if there is a tail of empty data remaining, flag them */
+    if (row0<iodata.tilesz*iodata.Nbase) {
+      for(int row = row0; row<iodata.tilesz*iodata.Nbase; row++) {
+        iodata.flag[row]=1;
+      }
+    }
+    for (int cm=0; cm<iodata.Nms;cm++) {
+     delete dataCols[cm];
+     delete flagCols[cm];
+    }
+    /* flagged data / total usable data, not counting excluded baselines */
+    if (countgood+countbad>0) {
+     *fratio=(double)countbad/(double)(countgood+countbad);
+    } else {
+     *fratio=1.0;
+    }
+}
+
+
+void 
+Data::writeData(Table ti, Data::IOData iodata) {
+
+    /* sort input table by ant1 and ant2 */
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=ti.sort(iv1,Sort::Ascending);
+
+    ROScalarColumn<int> a1(t, "ANTENNA1"), a2(t, "ANTENNA2");
+    /* writable access for output */
+    ArrayColumn<Complex> dataCol(t, Data::OutField);
+
+    /* check we get correct rows */
+    int nrow=t.nrow();
+    if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
+      cout<<"Warning: Missing rows, got "<<nrow<<" expect "<<iodata.tilesz*iodata.Nbase<<" +- "<<iodata.tilesz*iodata.N<<". (probably the last time interval, so not a big issue)."<<endl;
+    }
+    //cout<<"Table rows "<<nrow<<" Data rows "<<iodata.tilesz*iodata.Nbase+iodata.tilesz*iodata.N<<endl;
+    int row0=0;
+    IPosition pos(2,4,iodata.Nchan);
+    for(int row = 0; row < nrow; row++) {
+        uInt i = a1(row); //antenna1 
+        uInt j = a2(row); //antenna2
+        /* only work with cross correlations */
+        if (i!=j) {
+        Array<Complex> data = dataCol(row);
+        for(int k = 0; k < iodata.Nchan; k++) {
+           pos(0)=0;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8],iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+1]);
+           pos(0)=1;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+2],iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+3]);
+           pos(0)=2;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+4],iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+5]);
+           pos(0)=3;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+6],iodata.xo[iodata.Nbase*iodata.tilesz*8*k+row0*8+7]);
+       }
+
+       row0++;
+       dataCol.put(row,data); // copy to output
+      }
+    }
+}
+
+void 
+Data::writeDataList(vector<MSIter*> msitr, Data::IOData iodata) {
+    Table ti=msitr[0]->table();
+    /* sort input table by ant1 and ant2 */
+    Block<String> iv1(3);
+    iv1[0] = "TIME";
+    iv1[1] = "ANTENNA1";
+    iv1[2] = "ANTENNA2";
+    Table t=ti.sort(iv1,Sort::Ascending);
+
+    ROScalarColumn<int> a1(t, "ANTENNA1"), a2(t, "ANTENNA2");
+
+    /* check we get correct rows */
+    int nrow=t.nrow();
+    if(nrow-iodata.N*iodata.tilesz>iodata.tilesz*iodata.Nbase) {
+      cout<<"Error in rows"<<endl;
+    }
+    vector<ArrayColumn<Complex>* > dataCols(iodata.Nms);
+    for (int cm=0; cm<iodata.Nms;cm++) {
+     Table tti=(msitr[cm]->table());
+     Table *tt=new Table(tti.sort(iv1,Sort::Ascending));
+     dataCols[cm] = new  ArrayColumn<Complex>(*tt,Data::OutField);
+   }
+
+    int row0=0;
+    for(int row = 0; row < nrow; row++) {
+        uInt i = a1(row); //antenna1 
+        uInt j = a2(row); //antenna2
+        /* only work with cross correlations */
+        if (i!=j) {
+  int chanoff=0;
+  for (int cm=0; cm<iodata.Nms;cm++) {
+        Array<Complex> data = (*(dataCols[cm]))(row);
+        IPosition pos(2,4,iodata.NchanMS[cm]);
+        //Array<Complex> data = dataCol(row);
+        for(int k = 0; k < iodata.NchanMS[cm]; k++) {
+           pos(0)=0;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8],iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+1]);
+           pos(0)=1;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+2],iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+3]);
+           pos(0)=2;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+4],iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+5]);
+           pos(0)=3;pos(1)=k;
+           data(pos)=Complex(iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+6],iodata.xo[iodata.Nbase*iodata.tilesz*8*chanoff+row0*8+7]);
+           chanoff++;
+       }
+
+       (*(dataCols[cm])).put(row,data);
+       //dataCol.put(row,data); // copy to output
+
+   }
+       row0++;
+
+      }
+
+    }
+    for (int cm=0; cm<iodata.Nms;cm++) {
+     delete dataCols[cm];
+    }
+
+}
+
+
+void Data::freeData(Data::IOData data)
+{
+   delete [] data.u;
+   delete [] data.v;
+   delete [] data.w;
+   delete [] data.x;
+   delete [] data.xo;
+   delete [] data.freqs;
+   delete [] data.flag;
+   delete [] data.NchanMS;
+}
+
+
+void Data::freeData(Data::IOData data, Data::LBeam binfo)
+{
+   delete [] data.u;
+   delete [] data.v;
+   delete [] data.w;
+   delete [] data.x;
+   delete [] data.xo;
+   delete [] data.freqs;
+   delete [] data.flag;
+   delete [] data.NchanMS;
+
+   delete [] binfo.time_utc;
+   delete [] binfo.Nelem;
+   delete [] binfo.sx;
+   delete [] binfo.sy;
+   delete [] binfo.sz;
+   for (int ci=0; ci<data.N; ci++) {
+     delete [] binfo.xx[ci];
+     delete [] binfo.yy[ci];
+     delete [] binfo.zz[ci];
+   }
+   delete [] binfo.xx;
+   delete [] binfo.yy;
+   delete [] binfo.zz;
+}
diff --git a/src/MPI/data.h b/src/MPI/data.h
deleted file mode 120000
index 87a149c..0000000
--- a/src/MPI/data.h
+++ /dev/null
@@ -1 +0,0 @@
-../MS/data.h
\ No newline at end of file
diff --git a/src/MPI/data.h b/src/MPI/data.h
new file mode 100644
index 0000000..925aaad
--- /dev/null
+++ b/src/MPI/data.h
@@ -0,0 +1,174 @@
+/*
+ *
+ Copyright (C) 2006-2008 Sarod Yatawatta <sarod@users.sf.net>  
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ $Id$
+ */
+
+#ifndef __DATA_H__
+#define __DATA_H__
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <casacore/ms/MeasurementSets/MSIter.h>
+#include <casacore/tables/Tables/Table.h>
+#include <casacore/tables/Tables/TableVector.h>
+#include <casacore/tables/Tables/TableRecord.h>
+#include <casacore/tables/Tables/TableColumn.h>
+#include <casacore/tables/Tables/TableIter.h>
+#include <casacore/tables/Tables/ScalarColumn.h>
+#include <casacore/tables/Tables/ArrayColumn.h>
+#include <casacore/casa/Arrays/Array.h>
+#include <casacore/casa/Arrays/Cube.h>
+#include <fstream>
+#include <math.h>
+#include <complex>
+
+
+using namespace casacore;
+
+namespace Data
+{
+
+    struct IOData {
+       int N; /* no of stations */
+       int Nbase; /* baselines, exclude autocorrelations */
+       int tilesz;
+       int Nchan; /* total no of channels */
+       int Nms; /* no. of MS */
+       int *NchanMS; /* total channels per MS : Nms x 1 vector */
+       double deltat; /* integration time (s) */
+       int totalt; /* total no of time slots */
+       double ra0; /* phase center */
+       double dec0; /* phase center */
+       double *u; /* uvw coords, size Nbase*tilesz x 1 */
+       double *v;
+       double *w;
+       double *x; /* averaged data, size Nbase*8*tilez x 1
+       ordered by XX(re,im),XY(re,im),YX(re,im),YY(re,im), baseline, timeslots */
+       double *xo; /* unaveraged data, size Nbase*8*tilesz*Nchan x 1
+       ordered by  XX(re,im),XY(re,im),YX(re,im),YY(re,im), baseline, timeslots, channel */
+       double *freqs; /* channel freqs, size Nchan x 1 */
+       double freq0; /* averaged freq */
+       double *flag; /* double for conforming with old routines size Nbase*tilesz x 1 */
+       double deltaf; /* total bandwidth for freq. smearing */
+
+       double fratio; /* flagged data ratio = flagged/total, not counting data excluded from uv cut */
+                      /* if 1, all usable data are flagged */
+    };
+
+    /* Station beam info */
+    struct LBeam {
+      double *time_utc; /* time coord UTC (s), size tileszx1, 
+                         convert from MJD (s) to JD (days) */
+      int *Nelem; /* no of elements in each station, size Nx1 */
+      /* position (ITRF) of stations (m)
+       later changed to logitude,latitude,height (rad,rad,m) */
+      double *sx; /* x: size Nx1 */
+      double *sy; /* y: ... */
+      double *sz; /* z: ... */
+      /* x,y,z coords of elements, projected, converted to ITRF (m) */
+      double **xx; /* x coord pointer, size Nx1, each *x: x coord of station, size Nelem[]x1 */
+      double **yy; /* y ... */
+      double **zz; /* z ... */
+      /* pointing center of beams (only one) (could be different from phase center) */
+      double p_ra0;
+      double p_dec0;
+    };
+
+
+    /* read Auxilliary info and setup memory */
+    void readAuxData(const char *fname, IOData *data);
+    void readAuxData(const char *fname, IOData *data, LBeam *binfo);
+
+    void readAuxDataList(vector<string> msnames, IOData *data);
+
+    void readMSlist(char *fname, vector<string> *msnames);
+    /* load data using MS Iterator */
+    void loadData(Table t, IOData iodata, double *fratio);
+    void loadData(Table t, IOData iodata, LBeam binfo, double *fratio);
+
+    void loadDataList(vector<MSIter*> msitr, Data::IOData iodata, double *fratio);
+    /* write back data using MS Iterator */
+    void writeData(Table t, IOData iodata);
+    void writeDataList(vector<MSIter*> msitr, IOData iodata);
+    void freeData(IOData data);
+    void freeData(IOData data, LBeam binfo);
+
+    extern int numChannels; 
+    extern unsigned long int numRows;
+
+    struct float2 {
+      float x,y;
+    };
+
+
+    extern char *TableName; /* MS name */
+    extern char *MSlist; /* text file with MS names */
+    extern char *MSpattern; /* pattern to match all MS names used in calibration */
+    extern float min_uvcut;
+    extern float max_uvcut;
+    extern float max_uvtaper;
+    extern casacore::String DataField; /* input column DATA/CORRECTED_DATA */
+    extern casacore::String OutField; /* output column DATA/CORRECTED_DATA */
+    extern int TileSize; //Tile size
+    extern int Nt; /* no of worker threads */
+    extern char *SkyModel; /* sky model file */
+    extern char *Clusters; /* cluster file */
+    extern int format; /* sky model format 0: LSM, 1: LSM with 3 order spec idx*/
+ 
+    /* sagecal paramters */
+    extern int max_emiter;
+    extern int max_iter;
+    extern int max_lbfgs;
+    extern int lbfgs_m;
+    extern int gpu_threads;
+    extern int linsolv;
+    extern int solver_mode;
+    extern int ccid;
+    extern double rho;
+    extern char *solfile;
+    extern char *initsolfile;
+    extern char *ignorefile;
+    extern double nulow,nuhigh;
+    extern int randomize;
+    extern int whiten;
+    extern int DoSim; /* if 1, simulation mode */
+    extern int doChan; /* if 1, solve for each channel in multi channel data */
+    extern int doBeam; /* if 1, predict (LOFAR) beam array factor */
+    extern int DoDiag; /* if >0, enables diagnostics (Leverage) 1: write leverage as output (no residual), 2: only calculate fractions of leverage/noise */
+    extern int phaseOnly; /* if >0, and if any correction is done, extract phase and do phase only correction */
+
+    /* distributed sagecal parameters */
+    extern int Nadmm; /* ADMM iterations >=1 */
+    extern int Npoly; /* polynomial order >=1 */
+    extern int PolyType; /* what kind on polynomials to use 0,1,2,3 */
+    extern double admm_rho; /* regularization */
+    extern char *admm_rho_file; /* text file for regularization of each cluster */
+    extern int aadmm; /* if >0, enable adaptive update of rho */
+    /* for debugging, upper limit on time slots */
+    extern int Nmaxtime;
+    /* skipping initial timeslots */
+    extern int Nskip;
+    extern int verbose; /* if >0, enable verbose output */
+    extern int mdl; /* if given, calculate AIC/MDL for different poly configs and find minimum */
+    extern int GPUpredict; /* if given, use GPU for model calculation */
+    extern int heapsize; /* heap size in GPU (MB), for using malloc() */
+    /* for client server mode */
+    extern int servermode; /* 0: client, 1: server, else default operation */
+    extern char *servername; /* server host name or ip address */
+    extern char *portnumber; /* which port number to use for communication */
+}
+#endif //__DATA_H__
diff --git a/src/MS/CMakeLists.txt b/src/MS/CMakeLists.txt
index e7e79d6..8b32689 100644
--- a/src/MS/CMakeLists.txt
+++ b/src/MS/CMakeLists.txt
@@ -23,7 +23,7 @@ add_executable(sagecal ${SRCFILES})
 
 target_link_libraries(sagecal
  ${CASACORE_LIBRARIES}
- ${CFITSIO_LIBRARIES}
+ ${CFITSIO_LIB}
  ${OpenBLAS_LIB}
  ${LAPACK_LIBRARIES}
  ${WCSLIB_LIBRARIES}
@@ -32,7 +32,6 @@ target_link_libraries(sagecal
  -lopenblas
  -lgfortran
  -lpthread
- -lcfitsio
  -lm
  -ldirac
  -lradio
diff --git a/src/buildsky/CMakeLists.txt b/src/buildsky/CMakeLists.txt
index 47ceb0f..130ff45 100644
--- a/src/buildsky/CMakeLists.txt
+++ b/src/buildsky/CMakeLists.txt
@@ -1,5 +1,6 @@
 find_package(WcsLib REQUIRED)
-include_directories(${WCSLIB_INCLUDE_DIR}/wcslib)
+include_directories(${WCSLIB_INCLUDE_DIRS}/wcslib)
+include_directories(${WCSLIB_INCLUDE_DIRS})
 include_directories(./)
 
 link_directories(${LIBRARY_OUTPUT_PATH})
@@ -8,7 +9,7 @@ FILE(GLOB SRCFILES *.c)
 add_executable(buildsky ${SRCFILES})
 
 target_link_libraries(buildsky
- ${CFITSIO_LIBRARIES}
+ ${CFITSIO_LIB}
  ${OpenBLAS_LIB}
  ${LAPACK_LIBRARIES}
  ${WCSLIB_LIBRARIES}
@@ -17,7 +18,6 @@ target_link_libraries(buildsky
  -lopenblas
  -lgfortran
  -lpthread
- -lcfitsio
  -lm
 )
 
diff --git a/src/restore/CMakeLists.txt b/src/restore/CMakeLists.txt
index 85f5c95..b7f11ee 100644
--- a/src/restore/CMakeLists.txt
+++ b/src/restore/CMakeLists.txt
@@ -1,6 +1,7 @@
 find_package(WcsLib REQUIRED)
 find_package(FFTW REQUIRED)
-include_directories(${WCSLIB_INCLUDE_DIR}/wcslib)
+include_directories(${WCSLIB_INCLUDE_DIRS}/wcslib)
+include_directories(${WCSLIB_INCLUDE_DIRS})
 include_directories(./)
 
 link_directories(${LIBRARY_OUTPUT_PATH})
@@ -9,7 +10,7 @@ FILE(GLOB SRCFILES *.c)
 add_executable(restore ${SRCFILES})
 
 target_link_libraries(restore
- ${CFITSIO_LIBRARIES}
+ ${CFITSIO_LIB}
  ${OpenBLAS_LIB}
  ${LAPACK_LIBRARIES}
  ${WCSLIB_LIBRARIES}
@@ -20,7 +21,6 @@ target_link_libraries(restore
  -lopenblas
  -lgfortran
  -lpthread
- -lcfitsio
  -lm
 )