diff --git a/mssrc.py b/mssrc.py
index 8f8f44b93e87115399ed3b9306b317e1709a29bb..b02565f8c7248c7ce0e6d4d63d03bf633dd00dca 100644
--- a/mssrc.py
+++ b/mssrc.py
@@ -9,8 +9,11 @@ import numpy as np
 from numpy import linalg
 import pandas as pd
 import time
+import math
 import matplotlib as plt
 import sys
+import decimal
+from decimal import *
 try:
     from dataio.H5Writer import H5Writer
     from dataio.H5Reader import H5Reader
@@ -33,7 +36,34 @@ except ImportError:
 """
 
 # Define global variables
-COEFF_0 = 0.3   # Eq. (43) in Shen et al. (2005): low limit constant
+COEFF_0 = Decimal(0.3)   # Eq. (43) in Shen et al. (2005): low limit constant
+
+def inverse_den(x):
+    """
+    Calculate inverse of a number.
+    """
+    if abs(x) < Decimal(1.e-24):
+        return Decimal(0.0)
+    else:
+        return Decimal(1.0) / x
+
+def calc_det(mat):
+    """
+        Takes a matrix as and argument and returns the determinant
+
+        In this function, the configuration JSON-File and the the signal data are parsed. The from the signal time series, the local velocities, the void fraction, the bubble diameters and the interfacial area density are recovered by a reconstruction algorithm.
+    """
+    a = mat[0][0]
+    b = mat[0][1]
+    c = mat[0][2]
+    d = mat[1][0]
+    e = mat[1][1]
+    f = mat[1][2]
+    g = mat[2][0]
+    h = mat[2][1]
+    i = mat[2][2]
+    det = a*e*i + b*f*g + c*d*h - c*e*g - b*d*i - a*f*h
+    return det
 
 def main():
     """
@@ -56,7 +86,7 @@ def main():
     args = parser.parse_args()
 
     # Eq. (43) in Shen et al. (2005): estimated gas velocity
-    V_GAS = float(args.velocity)
+    V_GAS = Decimal(args.velocity)
 
     # Create Posix path for OS indepency
     path = pathlib.Path(args.path)
@@ -86,13 +116,16 @@ def main():
         if s_id != 0:
             aux_sensor_ids.append(s_id)
             # Set the distance vectors
-            S_0k[s_id] = np.asarray(sensor['relative_location'])
+            S_0k[s_id] = np.asarray(sensor['relative_location'],dtype='str')
             # Calculate magnitudes
-            S_0k_mag[s_id] = np.sqrt(S_0k[s_id][0]**2
-                                  + S_0k[s_id][1]**2
-                                  + S_0k[s_id][2]**2)
+            S_0k_mag[s_id] = (Decimal(S_0k[s_id][0]) * Decimal(S_0k[s_id][0])
+                            + Decimal(S_0k[s_id][1]) * Decimal(S_0k[s_id][1])
+                            + Decimal(S_0k[s_id][2]) * Decimal(S_0k[s_id][2])).sqrt()
             # Calculate cos of angles between vector S0k and (x,y,z)-axis
-            cos_eta_0k[s_id] = S_0k[s_id] / S_0k_mag[s_id]
+            tmp = []
+            for ii in range(0,3):
+                tmp.append(Decimal(S_0k[s_id][ii]) / S_0k_mag[s_id])
+            cos_eta_0k[s_id] = tmp
             # Determine valid range for lag between signal of sensors 0
             # and signal of sensor k from Eq. (43) in Shen et al. (2005)
             max_t_k[s_id] = S_0k_mag[s_id] / (COEFF_0 * V_GAS)
@@ -100,6 +133,7 @@ def main():
     # Create a H5-file reader
     reader = H5Reader(path / 'binary_signal.h5')
     # Read the time vector
+    t_signal_str = (reader.getDataSet('time')[:]).astype(str)
     t_signal = reader.getDataSet('time')[:]
     # Read the signal time series
     ds_signal = reader.getDataSet('signal')
@@ -115,6 +149,10 @@ def main():
     # of significance
     t_signal_ifd = t_signal[0:(len(t_signal)-1)] \
             + (t_signal[1:len(t_signal)]-t_signal[0:(len(t_signal)-1)])/2.0
+    t_signal_ifd_list = []
+    for ii in range(0,len(t_signal_ifd)):
+        t_signal_ifd_list.append(str(Decimal(t_signal_str[ii]) \
+            + (Decimal(t_signal_str[ii+1])-Decimal(t_signal_str[ii]))/Decimal(2.0)))
 
     # Interface-pairing signal-processing scheme
     # Shen et al. (2005)
@@ -131,10 +169,10 @@ def main():
     for ii,idx_rise_0 in enumerate(signal_ifd_rise_0):
         # DataFrame with time of IFD of front (2h) and rear (2h+1) bubble
         # interface for each sensor
-        ifd_times = pd.DataFrame(index=sensor_ids, columns=['t_2h','t_2h+1'],dtype='float64')
+        ifd_times = pd.DataFrame(index=sensor_ids, columns=['t_2h','t_2h+1'])
         # Main Sensor 0
         # Write t_2h of sensor 0
-        ifd_times['t_2h'][0] = t_signal_ifd[idx_rise_0]
+        ifd_times['t_2h'][0] = t_signal_ifd_list[idx_rise_0]
 
         # Search for t_2h+1
         # Indices of the falling IFD signals for main sensor 0
@@ -142,8 +180,7 @@ def main():
         # Index of the falling IFD signal for main sensor 0
         idx_fall_0 = min(signal_ifd_fall_0[signal_ifd_fall_0 > idx_rise_0])
         # Write t_b+1 of sensor 0
-        ifd_times['t_2h+1'][0] = t_signal_ifd[idx_fall_0]
-
+        ifd_times['t_2h+1'][0] = t_signal_ifd_list[idx_fall_0]
         # Auxillary sensors k
         for k in aux_sensor_ids:
             idk = np.where(sensor_ids == k)[0][0]
@@ -194,12 +231,12 @@ def main():
             # interface for each sensor
             # Write t_2h of sensor k
             if not np.isnan(idx_rise_k):
-                ifd_times['t_2h'][k] = t_signal_ifd[idx_rise_k]
+                ifd_times['t_2h'][k] = t_signal_ifd_list[idx_rise_k]
             else:
                 ifd_times['t_2h'][k] = np.nan
             # Write t_2h+1 of sensor k
             if not np.isnan(idx_fall_k):
-                ifd_times['t_2h+1'][k] = t_signal_ifd[idx_fall_k]
+                ifd_times['t_2h+1'][k] = t_signal_ifd_list[idx_fall_k]
             else:
                 ifd_times['t_2h+1'][k] = np.nan
 
@@ -208,12 +245,12 @@ def main():
         # and signal of sensor k from Eq. (42) in Shen et al. (2005)
         for k in aux_sensor_ids:
             # first interface (2h)
-            lag_2h = abs(ifd_times['t_2h'][0] - ifd_times['t_2h'][k])
+            lag_2h = abs(Decimal(ifd_times['t_2h'][0]) - Decimal(ifd_times['t_2h'][k]))
             if lag_2h > max_t_k[k]:
                 # lag too large
                 ifd_times['t_2h'][k] = np.nan
             # second interface (2h+1)
-            lag_2hp1 = abs(ifd_times['t_2h+1'][0] - ifd_times['t_2h+1'][k])
+            lag_2hp1 = abs(Decimal(ifd_times['t_2h+1'][0]) - Decimal(ifd_times['t_2h+1'][k]))
             if lag_2hp1 > max_t_k[k]:
                 ifd_times['t_2h+1'][k] = np.nan
 
@@ -233,7 +270,7 @@ def main():
     print(f"Detected {len(bubbles_incomplete)} incomplete bubble signals.\n")
     print('\nRunning reconstruction algorithm....\n')
     # Initialize the Interfacial Area Concentration (IAC)
-    iac = 0.0
+    iac = Decimal(0.0)
     # Loop over complete bubbles
     for ii,bubble_props in enumerate(bubbles_complete):
         ifd_times = bubble_props['ifd_times']
@@ -242,8 +279,7 @@ def main():
         delta_t_0k_2hp1 = {}
         # Time differences k=1,2,3,4 for interfaces 2h and 2h+1
         delta_t_k_h = {}
-        delta_t_k_h[0] = ifd_times['t_2h+1'][0] - ifd_times['t_2h'][0]
-
+        delta_t_k_h[0] = Decimal(ifd_times['t_2h+1'][0]) - Decimal(ifd_times['t_2h'][0])
         # 1. Reconstruct instantaneous local velocity
         # Inverse of measurable displacement velocity of the h-th bubble
         # for sensor 0 and k
@@ -264,17 +300,17 @@ def main():
         for k in aux_sensor_ids:
             # Time differences
             # Eq. (2) from Shen and Nakamura (2014)
-            delta_t_0k_2h[k] = ifd_times['t_2h'][k] \
-                             - ifd_times['t_2h'][0]
-            delta_t_0k_2hp1[k] = ifd_times['t_2h+1'][k] \
-                               - ifd_times['t_2h+1'][0]
+            delta_t_0k_2h[k] = Decimal(ifd_times['t_2h'][k]) \
+                             - Decimal(ifd_times['t_2h'][0])
+            delta_t_0k_2hp1[k] = Decimal(ifd_times['t_2h+1'][k]) \
+                               - Decimal(ifd_times['t_2h+1'][0])
             # Eq. (3) from Shen and Nakamura (2014)
-            delta_t_k_h[k] = ifd_times['t_2h+1'][k] \
-                           - ifd_times['t_2h'][k]
+            delta_t_k_h[k] = Decimal(ifd_times['t_2h+1'][k]) \
+                           - Decimal(ifd_times['t_2h'][k])
             # Measurable displacement velocity
             # Eq. (30) from Shen and Nakamura (2014)
             V_m0k_h_inv[k] = (delta_t_0k_2h[k] + delta_t_0k_2hp1[k]) \
-                        / (2.0 * S_0k_mag[k])
+                        / (Decimal(2.0) * Decimal(S_0k_mag[k]))
             A_0.append(cos_eta_0k[k])
             A_01_h.append(np.array([V_m0k_h_inv[k], cos_eta_0k[k][1], cos_eta_0k[k][2]]))
             A_02_h.append(np.array([cos_eta_0k[k][0], V_m0k_h_inv[k], cos_eta_0k[k][2]]))
@@ -282,13 +318,12 @@ def main():
 
         # Calculate matrix determinants
         # Eq. (32) to (35) from Shen and Nakamura (2014)
-        A_0_det = linalg.det(A_0)
-        A_01_h_det = linalg.det(A_01_h)
-        A_02_h_det = linalg.det(A_02_h)
-        A_03_h_det = linalg.det(A_03_h)
-
+        A_0_det = calc_det(A_0)
+        A_01_h_det = calc_det(A_01_h)
+        A_02_h_det = calc_det(A_02_h)
+        A_03_h_det = calc_det(A_03_h)
         # Calculate instantaneous velocity magnitude and components
-        dummy = np.sqrt(A_01_h_det**2 + A_02_h_det**2 + A_03_h_det**2)
+        dummy = (A_01_h_det* A_01_h_det + A_02_h_det*A_02_h_det + A_03_h_det*A_03_h_det).sqrt()
         # Eq. (37) from Shen and Nakamura (2014)
         V_bh_mag =  abs(A_0_det) * inverse_den(dummy)
         # Eq. (36) from Shen and Nakamura (2014)
@@ -310,10 +345,10 @@ def main():
         for k in aux_sensor_ids:
             # Chord lengths of interfaces 2h and 2h+1
             # Eq. (43) and (44) in Shen and Nakamura (2014)
-            E_0k_2h[k] = abs(V_bh_mag)**2 * delta_t_0k_2h[k] \
+            E_0k_2h[k] = V_bh_mag*V_bh_mag * delta_t_0k_2h[k] \
                         * (delta_t_0k_2hp1[k] + delta_t_k_h[0]) / S_0k_mag[k] \
                         - S_0k_mag[k]
-            E_0k_2hp1[k] = abs(V_bh_mag)**2 * delta_t_0k_2hp1[k] \
+            E_0k_2hp1[k] = V_bh_mag*V_bh_mag * delta_t_0k_2hp1[k] \
                         * (delta_t_0k_2h[k] - delta_t_k_h[0]) / S_0k_mag[k] \
                         - S_0k_mag[k]
             # Eq. (46) and (48) in Shen and Nakamura (2014)
@@ -325,17 +360,20 @@ def main():
             B_03_2hp1.append(np.array([cos_eta_0k[k][0], cos_eta_0k[k][1], E_0k_2hp1[k]]))
         # Calculate matrix determinants
         # Eq. (46) and (48) in Shen and Nakamura (2014)
-        B_01_2h_det = linalg.det(B_01_2h)
-        B_02_2h_det = linalg.det(B_02_2h)
-        B_03_2h_det = linalg.det(B_03_2h)
-        B_01_2hp1_det = linalg.det(B_01_2hp1)
-        B_02_2hp1_det = linalg.det(B_02_2hp1)
-        B_03_2hp1_det = linalg.det(B_03_2hp1)
+        B_01_2h_det = calc_det(B_01_2h)
+        B_02_2h_det = calc_det(B_02_2h)
+        B_03_2h_det = calc_det(B_03_2h)
+        B_01_2hp1_det = calc_det(B_01_2hp1)
+        B_02_2hp1_det = calc_det(B_02_2hp1)
+        B_03_2hp1_det = calc_det(B_03_2hp1)
 
         # Calculate the bubble diameter
-        dummy_2h = np.sqrt(B_01_2h_det**2 + B_02_2h_det**2 + B_03_2h_det**2)
-        dummy_2hp1 = np.sqrt(B_01_2hp1_det**2 + B_02_2hp1_det**2 \
-                           + B_03_2hp1_det**2)
+        dummy_2h = (B_01_2h_det*B_01_2h_det \
+                    + B_02_2h_det*B_02_2h_det \
+                    + B_03_2h_det*B_03_2h_det).sqrt()
+        dummy_2hp1 = (B_01_2hp1_det*B_01_2hp1_det \
+                    + B_02_2hp1_det*B_02_2hp1_det \
+                    + B_03_2hp1_det*B_03_2hp1_det).sqrt()
         # Eq. (50) from Shen and Nakamura (2014)
         D_h_2h = dummy_2h / A_0_det
         D_h_2hp1 = dummy_2hp1 / A_0_det
@@ -343,12 +381,12 @@ def main():
 
         # Calculate interfacial normal unit vectors
         # Eq. (51) in Shen and Nakamura (2014)
-        cos_eta_i_2h = np.ones(3) * A_0_det \
+        cos_eta_i_2h = np.array([Decimal(1.0), Decimal(1.0), Decimal(1.0)]) * A_0_det \
                          * inverse_den(abs(A_0_det)) \
                          * inverse_den(dummy_2h)
         cos_eta_i_2h = cos_eta_i_2h \
                     * np.array([B_01_2h_det, B_02_2h_det, B_03_2h_det])
-        cos_eta_i_2hp1 = np.ones(3) * A_0_det \
+        cos_eta_i_2hp1 = np.array([Decimal(1.0), Decimal(1.0), Decimal(1.0)]) * A_0_det \
                         * inverse_den(abs(A_0_det)) \
                         * inverse_den(dummy_2hp1)
         cos_eta_i_2hp1 = cos_eta_i_2hp1 \
@@ -359,14 +397,16 @@ def main():
         # Calculate the local Interfacial Area Concentration (IAC)
         # Eq. (54) in Shen and Nakamura (2014)
         # Sampling duration Omega
-        Omega = max(t_signal) - min(t_signal)
-        dummy_A = A_01_h_det**2 + A_02_h_det**2 + A_03_h_det**2
-        dummy_B = np.sqrt(B_01_2h_det**2 + B_02_2h_det**2 + B_03_2h_det**2)
+        Omega = Decimal(t_signal[-1]) - Decimal(t_signal[0])
+        dummy_A = A_01_h_det*A_01_h_det + A_02_h_det*A_02_h_det + A_03_h_det*A_03_h_det
+        dummy_B = (B_01_2h_det*B_01_2h_det \
+                    + B_02_2h_det*B_02_2h_det \
+                    + B_03_2h_det*B_03_2h_det).sqrt()
         dummy_AB = abs(A_0_det * (A_01_h_det*B_01_2h_det \
                                 + A_02_h_det*B_02_2h_det \
                                 + A_03_h_det*B_03_2h_det))
         iac_h = dummy_A * dummy_B * inverse_den(dummy_AB)
-        iac += 2.0 / Omega * iac_h
+        iac += Decimal(2.0) / Omega * iac_h
         # Display progress
         if args.progress:
             printProgressBar(ii, len(bubbles_complete), prefix = 'Progress:', suffix = 'Complete', length = 50)
@@ -436,7 +476,7 @@ def main():
     # Add the attributes
     ds_d.attrs['labels'] = ['D_h_2h', 'D_h_2hp1']
     # Create the IAC and void_fraction datasets
-    writer.writeDataSet('IAC', np.array([iac]), 'float64')
+    writer.writeDataSet('IAC', np.array([iac],dtype='float64'), 'float64')
     writer.writeDataSet('voidFraction', np.array([np.mean(signal)]), 'float64')
     writer.close()
     time2 = time.time()