diff --git a/Plot-Scripts/WoodBilayer_ExperimentComparison.py b/Plot-Scripts/WoodBilayer_ExperimentComparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..2a2d68d0700f6464ca2d763ce82df78b573f3c75
--- /dev/null
+++ b/Plot-Scripts/WoodBilayer_ExperimentComparison.py
@@ -0,0 +1,199 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import os
+import subprocess
+import fileinput
+import re
+import sys
+import matplotlib as mpl
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+import matplotlib.ticker as ticker
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import seaborn as sns
+import matplotlib.colors as mcolors
+# ----------------------------------------------------
+# --- Define Plot style:
+# plt.style.use("seaborn-white")
+# plt.style.use("seaborn-pastel")
+# plt.style.use("seaborn-colorblind")
+plt.style.use("seaborn")
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "14"
+mpl.rcParams['xtick.bottom'] = True
+mpl.rcParams['xtick.major.size'] = 3
+mpl.rcParams['xtick.minor.size'] = 1.5
+mpl.rcParams['xtick.major.width'] = 0.75
+mpl.rcParams['ytick.left'] = True
+mpl.rcParams['ytick.major.size'] = 3
+mpl.rcParams['ytick.minor.size'] = 1.5
+mpl.rcParams['ytick.major.width'] = 0.75
+
+#Adjust grid:
+mpl.rcParams.update({"axes.grid" : True}) # Add grid
+mpl.rcParams['axes.labelpad'] = 3
+mpl.rcParams['grid.linewidth'] = 0.25
+mpl.rcParams['grid.alpha'] = 0.9 # 0.75
+mpl.rcParams['grid.linestyle'] = '-'
+mpl.rcParams['grid.color']   = 'gray'#'black'
+mpl.rcParams['text.latex.preamble'] = r'\usepackage{amsfonts}' # Makes Use of \mathbb possible.
+# ----------------------------------------------------------------------------------------
+width = 5.79
+height = width / 1.618 # The golden ratio.
+
+
+
+
+#TODO read values from .txt output files....
+
+dataset_numbers = [0, 1, 2, 3, 4 ,5]
+
+for dataset_number in dataset_numbers:
+
+    # fig = plt.figure()      #main
+    fig, ax = plt.subplots(figsize=(width,height))
+
+    dataset = [
+
+        [   # Dataset Ratio r = 0.12
+        [14.70179844, 13.6246    , 12.42994508, 11.69773413,11.14159987, 9.500670278, 9.005046347] ,    # input (moisture-content in %)
+        [1.29323308,  1.83458647 , 2.40601504 , 2.76691729 , 3.03759398, 3.81954887 , 4.03007519],   # curvature kappa from Simulation
+        [1.140351217, 1.691038688, 2.243918105, 2.595732726, 2.945361006,4.001528043, 4.312080261] # curvature kappa from Experiment
+        ],
+
+        [   # Dataset Ratio r = 0.17
+        [14.75453569,13.71227639,12.54975012,11.83455959,11.29089521,9.620608917,9.101671742],   # input (moisture-content in %)
+        [1.26315789, 1.77443609, 2.34586466, 2.70676692, 2.97744361, 3.78947368,4.03007519],   # curvature kappa from Simulation
+        [1.02915975,1.573720805,2.407706364,2.790518802,3.173814476,4.187433094,4.511739072]    # curvature kappa from Experiment
+        ],
+
+        [   # Dataset Ratio r = 0.22
+        [14.72680026, 13.64338887, 12.41305478, 11.66482931, 11.09781471, 9.435795985, 8.959564147],   # input (moisture-content in %)
+        [1.20401338 , 1.72575251 , 2.28762542 , 2.64882943 , 2.92976589 , 3.73244147 , 3.93311037 ],   # curvature kappa from Simulation
+        [1.058078122, 1.544624544, 2.317033799, 2.686043143, 2.967694189, 3.913528418, 4.262750825]    # curvature kappa from Experiment
+        ],
+
+        [   # Dataset Ratio r = 0.34
+        [14.98380876,13.97154915,12.77309253,12.00959929,11.42001731,9.561447179,8.964704969],   # input (moisture-content in %)
+        [0.87218045, 1.26315789, 1.7443609,  2.07518797, 2.28571429, 3.03759398, 3.27819549],   # curvature kappa from Simulation
+        [0.789078472,1.1299263,1.738136936,2.159520896,2.370047499,3.088299431,3.18097558]    # curvature kappa from Experiment
+        ],
+
+        [   # Dataset Ratio r = 0.43
+        [15.11316339,14.17997082,13.05739844,12.32309209,11.74608518,9.812372466,9.10519385 ],   # input (moisture-content in %)
+        [0.63157895, 0.93233083, 1.29323308, 1.53383459, 1.71428571, 2.34586466,2.55639098],   # curvature kappa from Simulation
+        [0.577989364,0.829007544,1.094211707,1.325332511,1.400455154,1.832325697,2.047483977]  # curvature kappa from Experiment
+        ],
+
+        [   # Dataset Ratio r = 0.49
+        [15.30614414,14.49463867,13.46629742,12.78388234,12.23057715,10.21852839,9.341730605],   # input (moisture-content in %)
+        [0.60150376, 0.87218045, 1.23308271, 1.5037594 , 1.71428571, 2.46616541,2.79699248],   # curvature kappa from Simulation
+        [0.357615902,0.376287785,0.851008627,0.904475291,1.039744708,1.346405241,1.566568558]   # curvature kappa from Experiment
+        ]
+    ]
+
+
+
+    # input = [14.72680026, 13.64338887, 12.41305478, 11.66482931, 11.09781471, 9.435795985, 8.959564147]  #moisture-content (in %)
+    # kappa_sim = [1.20401338 , 1.72575251, 2.28762542 , 2.64882943 , 2.92976589 , 3.73244147 , 3.93311037 ]    #kappa from Simulation
+    # kappa_exp = [1.058078122,1.544624544, 2.317033799, 2.686043143, 2.967694189, 3.913528418, 4.262750825]    #kappa from Experiment
+
+    # compute difference:
+    # relative_error = (np.array(kappa_sim) - np.array(kappa_exp)) / np.array(kappa_exp)
+    # print('relative_error:', relative_error)
+
+    relative_error = (np.array(dataset[dataset_number][1]) - np.array(dataset[dataset_number][2])) / np.array(dataset[dataset_number][2])
+    print('relative_error:', relative_error)
+
+
+    # --------------- Plot Lines + Scatter -----------------------
+    line_1 = ax.plot(np.array(dataset[dataset_number][0]), np.array(dataset[dataset_number][1]),                    # data
+                #  color='forestgreen',              # linecolor
+                marker='D',                         # each marker will be rendered as a circle
+                markersize=5,                       # marker size
+                #   markerfacecolor='darkorange',      # marker facecolor
+                markeredgecolor='black',            # marker edgecolor
+                markeredgewidth=0.75,                  # marker edge width
+                # linestyle='dashdot',              # line style will be dash line
+                linewidth=1.5,                      # line width
+                zorder=3,
+                label = r"$\kappa_{sim}$")
+
+    line_2 = ax.plot(np.array(dataset[dataset_number][0]), np.array(dataset[dataset_number][2]),                     # data
+                # color='orangered',                # linecolor
+                marker='o',                         # each marker will be rendered as a circle
+                markersize=5,                       # marker size
+                #  markerfacecolor='cornflowerblue',   # marker facecolor
+                markeredgecolor='black',            # marker edgecolor
+                markeredgewidth=0.75,                  # marker edge width
+                # linestyle='--',                   # line style will be dash line
+                linewidth=1.5,                      # line width
+                zorder=3,
+                alpha=0.8,                           # Change opacity
+                label = r"$\kappa_{exp}$")
+
+    # --- Plot order line
+    # x = np.linspace(0.01,1/2,100)
+    # y = CC_L2[0]*x**2
+    # OrderLine = ax.plot(x,y,linestyle='--', label=r"$\mathcal{O}(h)$")
+
+
+
+    # Fix_value = 7.674124
+    # l3 = plt.axhline(y = Fix_value, color = 'black', linewidth=0.75, linestyle = 'dashed')
+    # --------------- Set Axes  -----------------------
+    # ax.set_title(r"ratio $r = 0.22$")   # Plot - Title
+
+    # Plot - Title
+    match dataset_number:
+        case 0:
+            ax.set_title(r"ratio $r = 0.12$") 
+        case 1: 
+            ax.set_title(r"ratio $r = 0.17$") 
+        case 2:
+            ax.set_title(r"ratio $r = 0.22$") 
+        case 3: 
+            ax.set_title(r"ratio $r = 0.34$") 
+        case 4:
+            ax.set_title(r"ratio $r = 0.43$") 
+        case 5: 
+            ax.set_title(r"ratio $r = 0.49$") 
+    
+    # plt.xscale('log') # Use Logarithmic-Scale
+    # plt.yscale('log')
+    ax.set_xlabel(r"Wood moisture content $\omega (\%)$", labelpad=4)
+    ax.set_ylabel(r"Curvature $\kappa$", labelpad=4)
+    plt.tight_layout()
+
+    # # --- Set Line labels
+    # line_labels = [r"$CC_{L_2}$",r"$CC_{H_1}$", r"$\mathcal{O}(h)$"]
+
+    # --- Set Legend
+    legend = ax.legend()
+    # legend = fig.legend([line_1 , line_2, OrderLine],
+    #                     labels = line_labels,
+    #                     bbox_to_anchor=[0.97, 0.50],
+    #                     # bbox_to_anchor=[0.97, 0.53],
+    #                     # loc='center',
+    #                     ncol=1,                  # Number of columns used for legend
+    #                     # borderaxespad=0.15,    # Small spacing around legend box
+    #                     frameon=True,
+    #                     prop={'size': 10})
+
+
+    frame = legend.get_frame()
+    frame.set_edgecolor('black')
+
+
+    # --- Adjust left/right spacing:
+    # plt.subplots_adjust(right=0.81)
+    # plt.subplots_adjust(left=0.11)
+
+    # ---------- Output Figure as pdf:
+    fig.set_size_inches(width, height)
+    fig.savefig('WoodBilayer_expComparison_' + str(dataset_number) + '.pdf')
+    # plt.show()
+
+
diff --git a/experiment/wood-bilayer/PolarPlotLocalEnergy.py b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
index 6e623b15cf0658cd74b5fcb3bbe78ed0100a8d2a..daa490a07d3db90a976ef6131c482192efffdd7f 100644
--- a/experiment/wood-bilayer/PolarPlotLocalEnergy.py
+++ b/experiment/wood-bilayer/PolarPlotLocalEnergy.py
@@ -48,48 +48,63 @@ def ReadEffectiveQuantities(QFilePath, BFilePath):
     B = np.array([X[0][2], X[1][2], X[2][2]])
     return Q, B
 
+# -------------------------------------------------------------------
+
+
 # Number of experiments / folders
 number=7
-kappa=np.zeros(number)
-for n in range(0,number):
-    #   Read from Date
-    print(str(n))
-    DataPath = './experiment/wood-bilayer/results/'+str(n)
-    DataPath = './experiment/wood-bilayer/results/'+str(n)
-    QFilePath = DataPath + '/QMatrix.txt'
-    BFilePath = DataPath + '/BMatrix.txt'
-    Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
-    # Q=0.5*(np.transpose(Q)+Q) # symmetrize
-    B=np.transpose([B])
-    # 
-    
-    N=300
-    length=12
-    r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
-    E=np.zeros(np.shape(r))
-    for i in range(0,N): 
-        for j in range(0,N):     
-            if theta[i,j]<np.pi:
-                E[i,j]=energy(r[i,j],theta[i,j],Q,B)
-            else:
-                E[i,j]=energy(-r[i,j],theta[i,j],Q,B)
-            
-    # Compute Minimizer
-    [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
-    kappamin=r[imin,jmin]
-    alphamin=theta[imin,jmin]
-    kappa[n]=kappamin
-    fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
-    levs=np.geomspace(E.min(),E.max(),400)
-    pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
-    ax.set_xticks([0,np.pi/2])
-    ax.set_yticks([kappamin])
-    colorbarticks=np.linspace(E.min(),E.max(),6)
-    plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
-    plt.show()
-
-f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
-f.write(str(kappa))         
-f.close()   
+show_plot = False
+
+dataset_numbers = [0, 1, 2, 3, 4, 5]
+
+for dataset_number in dataset_numbers:
+
+    kappa=np.zeros(number)
+    for n in range(0,number):
+        #   Read from Date
+        print(str(n))
+        DataPath = './experiment/wood-bilayer/results_'+ str(dataset_number) + '/' +str(n)
+        QFilePath = DataPath + '/QMatrix.txt'
+        BFilePath = DataPath + '/BMatrix.txt'
+        Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+        # Q=0.5*(np.transpose(Q)+Q) # symmetrize
+        B=np.transpose([B])
+        # 
+        
+        N=400
+        length=12
+        r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
+        E=np.zeros(np.shape(r))
+        for i in range(0,N): 
+            for j in range(0,N):     
+                if theta[i,j]<np.pi:
+                    E[i,j]=energy(r[i,j],theta[i,j],Q,B)
+                else:
+                    E[i,j]=energy(-r[i,j],theta[i,j],Q,B)
+                
+        # Compute Minimizer
+        [imin,jmin]=np.unravel_index(E.argmin(),(N,N))
+        kappamin=r[imin,jmin]
+        alphamin=theta[imin,jmin]
+        kappa[n]=kappamin
+        fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
+        levs=np.geomspace(E.min(),E.max(),400)
+        pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
+        ax.set_xticks([0,np.pi/2])
+        ax.set_yticks([kappamin])
+        colorbarticks=np.linspace(E.min(),E.max(),6)
+        plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
+        if (show_plot):
+            plt.show()
+        # Save Figure as .pdf
+        width = 5.79
+        height = width / 1.618 # The golden ratio.
+        fig.set_size_inches(width, height)
+        fig.savefig('Plot_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')
+
+    # f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
+    f = open("./experiment/wood-bilayer/results_" + str(dataset_number) +  "/kappa_simulation.txt", "w")
+    f.write(str(kappa))         
+    f.close()   
 
 
diff --git a/experiment/wood-bilayer/wood_european_beech.py b/experiment/wood-bilayer/wood_european_beech.py
index 135ca8b24ef61d2fa0255e4c567bf283fd22d875..89da18ae223f5c8e5a6d38843ca0237fd4a1c0e0 100644
--- a/experiment/wood-bilayer/wood_european_beech.py
+++ b/experiment/wood-bilayer/wood_european_beech.py
@@ -42,13 +42,13 @@ parameterSet.Phases=2
 
 # Parameters of the model
 # -- (thickness upper layer) / (thickness)
-param_r = 0.22
+param_r = 0.49
 # -- thickness [meter]
-param_h = 0.0053
+param_h = 0.008
 # -- moisture content in the flat state [%]
-param_omega_flat = 17.17547062
+param_omega_flat = 17.01520754
 # -- moisture content in the target state [%]
-param_omega_target = 8.959564147
+param_omega_target = 9.341730605
 # -- Drehwinkel
 param_theta = 0
 
diff --git a/experiment/wood-bilayer/wood_test.py b/experiment/wood-bilayer/wood_test.py
index ad7610cb6e6adc8dc51c9bfc6a8b5c216a01ce47..1e96127ac37496864160a9e4e1131ca8758c884c 100644
--- a/experiment/wood-bilayer/wood_test.py
+++ b/experiment/wood-bilayer/wood_test.py
@@ -101,124 +101,239 @@ def eval_energy(kappa,alpha,Q,B)  :
 ########################
 #### SET PARAMETERS ####
 ########################
-# ----- Setup Paths -----
-# write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
-# path='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer/results/'  
-# pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer'
-path = os.getcwd() + '/experiment/wood-bilayer/results/'
-pythonPath = os.getcwd() + '/experiment/wood-bilayer'
-pythonModule = "wood_european_beech"
-executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
-# ---------------------------------
-# Setup Experiment
-# ---------------------------------
-gamma = 1.0
-
-# ----- Define Parameters for Material Function  --------------------
-# [r, h, omega_flat, omega_target, theta, experimental_kappa]
-# r = (thickness upper layer)/(thickness)
-# h = thickness [meter]
-# omega_flat = moisture content in the flat state before drying [%]
-# omega_target = moisture content in the target state [%]
-# theta = rotation angle (not implemented and used)
-# experimental_kappa = curvature measure in experiment
-
-#First Experiment:
-# Ratio r = 0.22
-# materialFunctionParameter=[
-#    [0.22, 0.0053, 17.17547062, 14.72680026, 0, 1.058078122],
-#    [0.22, 0.0053, 17.17547062, 13.64338887, 0, 1.544624544],
-#    [0.22, 0.0053, 17.17547062, 12.41305478, 0, 2.317033799],
-#    [0.22, 0.0053, 17.17547062, 11.66482931, 0, 2.686043143],
-#    [0.22, 0.0053, 17.17547062, 11.09781471, 0, 2.967694189],
-#    [0.22, 0.0053, 17.17547062, 9.435795985, 0, 3.913528418],
-#    [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825]
-# ]
-
-# Ratio r = 0.49
-# materialFunctionParameter=[
-#    [0.49, 0.008,  17.01520754, 15.30614414, 0, 0.357615902],
-#    [0.49, 0.008,  17.01520754, 14.49463867, 0, 0.376287785],
-#    [0.49, 0.008,  17.01520754, 13.46629742, 0, 0.851008627],
-#    [0.49, 0.008,  17.01520754, 12.78388234, 0, 0.904475291],
-#    [0.49, 0.008,  17.01520754, 12.23057715, 0, 1.039744708],
-#    [0.49, 0.008,  17.01520754, 10.21852839, 0, 1.346405241],
-#    [0.49, 0.008,  17.01520754, 9.341730605, 0, 1.566568558]
-# ]
-
-#Second Experiment: Rotate "active" bilayer phase 
-# materialFunctionParameter=[
-#     [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/6.0), 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/3.0), 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/2.0), 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, 2.0*(np.pi/3.0), 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, 5.0*(np.pi/6.0), 4.262750825],
-#     [0.22, 0.0053,  17.17547062, 8.959564147, np.pi, 4.262750825]
-# ]
-
-materialFunctionParameter=[
-    [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/12.0), 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/6.0), 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/4.0), 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/3.0), 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, 5.0*(np.pi/12.0), 4.262750825],
-    [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/2.0), 4.262750825]
-]
-
-# ------ Loops through Parameters for Material Function -----------
-for i in range(0,np.shape(materialFunctionParameter)[0]):
+
+dataset_numbers = [0, 1, 2, 3, 4, 5]
+
+for dataset_number in dataset_numbers:
     print("------------------")
-    print("New Loop")
+    print(str(dataset_number) + "th data set")
     print("------------------")
-   # Check output directory
-    outputPath = path + str(i)
-    isExist = os.path.exists(outputPath)
-    if not isExist:
-        # Create a new directory because it does not exist
-        os.makedirs(outputPath)
-        print("The new directory " + outputPath + " is created!")
-
-    # thread = threading.Thread(target=run_CellProblem(executable, pythonModule, pythonPath, LOGFILE))
-    # thread.start()
-
-    #TODO: apperently its not possible to pass a variable via subprocess and "calculate" another input value inside the python file.
-    #      Therefore we use this instead.
-    SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_r",materialFunctionParameter[i][0])
-    SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_h",materialFunctionParameter[i][1])
-    SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_flat",materialFunctionParameter[i][2])
-    SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_target",materialFunctionParameter[i][3])
-    SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_theta",materialFunctionParameter[i][4])    
-
-    LOGFILE = outputPath + "/" + pythonModule + "_output" + "_" + str(i) + ".log"
-
-    processList = []
-    p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
-                                    + " -outputPath " + outputPath
-                                    + " -gamma " + str(gamma) 
-                                    + " | tee " + LOGFILE, shell=True)
-
-    # p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
-    #                                 + " -outputPath " + outputPath
-    #                                 + " -gamma " + str(gamma) 
-    #                                 + " -param_r " + str(materialFunctionParameter[i][0])
-    #                                 + " -param_h " + str(materialFunctionParameter[i][1])
-    #                                 + " -param_omega_flat " + str(materialFunctionParameter[i][2])
-    #                                 + " -param_omega_target " + str(materialFunctionParameter[i][3])
-    #                                 + " -phase2_angle " + str(materialFunctionParameter[i][4])
-    #                                 + " | tee " + LOGFILE, shell=True)
-
-    p.wait() # wait
-    processList.append(p)
-    exit_codes = [p.wait() for p in processList]
-    # ---------------------------------------------------
-    # wait here for the result to be available before continuing
-    # thread.join()
-    f = open(outputPath+"/parameter.txt", "w")
-    f.write("r = "+str(materialFunctionParameter[i][0])+"\n")
-    f.write("h = "+str(materialFunctionParameter[i][1])+"\n")
-    f.write("omega_flat = "+str(materialFunctionParameter[i][2])+"\n")        
-    f.write("omega_target = "+str(materialFunctionParameter[i][3])+"\n")         
-    f.close()   
-    #
+
+    # ----- Setup Paths -----
+    # write_LOG = True   # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+    # path='/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer/results/'  
+    # pythonPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/experiment/wood-bilayer'
+    path = os.getcwd() + '/experiment/wood-bilayer/results_' + str(dataset_number) + '/'
+    pythonPath = os.getcwd() + '/experiment/wood-bilayer'
+    pythonModule = "wood_european_beech"
+    executable = os.getcwd() + '/build-cmake/src/Cell-Problem'
+    # ---------------------------------
+    # Setup Experiment
+    # ---------------------------------
+    gamma = 1.0
+
+    # ----- Define Parameters for Material Function  --------------------
+    # [r, h, omega_flat, omega_target, theta, experimental_kappa]
+    # r = (thickness upper layer)/(thickness)
+    # h = thickness [meter]
+    # omega_flat = moisture content in the flat state before drying [%]
+    # omega_target = moisture content in the target state [%]
+    # theta = rotation angle (not implemented and used)
+    # experimental_kappa = curvature measure in experiment
+
+    #First Experiment:
+
+
+    # Dataset Ratio r = 0.12
+    # materialFunctionParameter=[
+    #    [0.12, 0.0047, 17.32986047, 14.70179844, 0, 1.140351217],
+    #    [0.12, 0.0047, 17.32986047, 13.6246,     0, 1.691038688],
+    #    [0.12, 0.0047, 17.32986047, 12.42994508, 0, 2.243918105],
+    #    [0.12, 0.0047, 17.32986047, 11.69773413, 0, 2.595732726],
+    #    [0.12, 0.0047, 17.32986047, 11.14159987, 0, 2.945361006],
+    #    [0.12, 0.0047, 17.32986047, 9.500670278, 0, 4.001528043],
+    #    [0.12, 0.0047, 17.32986047, 9.005046347, 0, 4.312080261],
+    # ]
+
+    # # Dataset Ratio r = 0.17 
+    # materialFunctionParameter=[
+    #    [0.17, 0.0049, 17.28772791 , 14.75453569, 0, 1.02915975],
+    #    [0.17, 0.0049, 17.28772791 , 13.71227639,  0, 1.573720805],
+    #    [0.17, 0.0049, 17.28772791 , 12.54975012, 0, 2.407706364],
+    #    [0.17, 0.0049, 17.28772791 , 11.83455959, 0, 2.790518802],
+    #    [0.17, 0.0049, 17.28772791 , 11.29089521, 0, 3.173814476],
+    #    [0.17, 0.0049, 17.28772791 , 9.620608917, 0, 4.187433094],
+    #    [0.17, 0.0049, 17.28772791 , 9.101671742, 0, 4.511739072],
+    # ]
+
+    # # Dataset Ratio r = 0.22
+    # materialFunctionParameter=[
+    #    [0.22, 0.0053,  17.17547062, 14.72680026, 0, 1.058078122],
+    #    [0.22, 0.0053,  17.17547062, 13.64338887, 0, 1.544624544],
+    #    [0.22, 0.0053,  17.17547062, 12.41305478, 0, 2.317033799],
+    #    [0.22, 0.0053,  17.17547062, 11.66482931, 0, 2.686043143],
+    #    [0.22, 0.0053,  17.17547062, 11.09781471, 0, 2.967694189],
+    #    [0.22, 0.0053,  17.17547062, 9.435795985, 0, 3.913528418],
+    #    [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825]
+    # ]
+
+    # # Dataset Ratio r = 0.34
+    # materialFunctionParameter=[
+    #    [0.34, 0.0063, 17.14061081 , 14.98380876, 0, 0.789078472],
+    #    [0.34, 0.0063, 17.14061081 , 13.97154915  0, 1.1299263],
+    #    [0.34, 0.0063, 17.14061081 , 12.77309253, 0, 1.738136936],
+    #    [0.34, 0.0063, 17.14061081 , 12.00959929, 0, 2.159520896],
+    #    [0.34, 0.0063, 17.14061081 , 11.42001731, 0, 2.370047499],
+    #    [0.34, 0.0063, 17.14061081 , 9.561447179, 0, 3.088299431],
+    #    [0.34, 0.0063, 17.14061081 , 8.964704969, 0, 3.18097558],
+    # ]
+
+    # # Dataset Ratio r = 0.43
+    # materialFunctionParameter=[
+    #    [0.43, 0.0073, 17.07559686 , 15.11316339, 0, 0.577989364],
+    #    [0.43, 0.0073, 17.07559686 , 14.17997082, 0, 0.829007544],
+    #    [0.43, 0.0073, 17.07559686 , 13.05739844, 0, 1.094211707],
+    #    [0.43, 0.0073, 17.07559686 , 12.32309209, 0, 1.325332511],
+    #    [0.43, 0.0073, 17.07559686 , 11.74608518, 0, 1.400455154],
+    #    [0.43, 0.0073, 17.07559686 , 9.812372466, 0, 1.832325697],
+    #    [0.43, 0.0073, 17.07559686 , 9.10519385 , 0, 2.047483977],
+    # ]
+
+    # # Dataset Ratio r = 0.49
+    # materialFunctionParameter=[
+    #    [0.49, 0.008,  17.01520754, 15.30614414, 0, 0.357615902],
+    #    [0.49, 0.008,  17.01520754, 14.49463867, 0, 0.376287785],
+    #    [0.49, 0.008,  17.01520754, 13.46629742, 0, 0.851008627],
+    #    [0.49, 0.008,  17.01520754, 12.78388234, 0, 0.904475291],
+    #    [0.49, 0.008,  17.01520754, 12.23057715, 0, 1.039744708],
+    #    [0.49, 0.008,  17.01520754, 10.21852839, 0, 1.346405241],
+    #    [0.49, 0.008,  17.01520754, 9.341730605, 0, 1.566568558]
+    # ]
+
+
+
+
+
+    materialFunctionParameter=[
+    [  # Dataset Ratio r = 0.12
+    [0.12, 0.0047, 17.32986047, 14.70179844, 0, 1.140351217],
+    [0.12, 0.0047, 17.32986047, 13.6246,     0, 1.691038688],
+    [0.12, 0.0047, 17.32986047, 12.42994508, 0, 2.243918105],
+    [0.12, 0.0047, 17.32986047, 11.69773413, 0, 2.595732726],
+    [0.12, 0.0047, 17.32986047, 11.14159987, 0, 2.945361006],
+    [0.12, 0.0047, 17.32986047, 9.500670278, 0, 4.001528043],
+    [0.12, 0.0047, 17.32986047, 9.005046347, 0, 4.312080261]
+    ],
+    [  # Dataset Ratio r = 0.17
+    [0.17, 0.0049, 17.28772791 , 14.75453569, 0, 1.02915975],
+    [0.17, 0.0049, 17.28772791 , 13.71227639,  0, 1.573720805],
+    [0.17, 0.0049, 17.28772791 , 12.54975012, 0, 2.407706364],
+    [0.17, 0.0049, 17.28772791 , 11.83455959, 0, 2.790518802],
+    [0.17, 0.0049, 17.28772791 , 11.29089521, 0, 3.173814476],
+    [0.17, 0.0049, 17.28772791 , 9.620608917, 0, 4.187433094],
+    [0.17, 0.0049, 17.28772791 , 9.101671742, 0, 4.511739072]
+    ],
+    [  # Dataset Ratio r = 0.22
+    [0.22, 0.0053,  17.17547062, 14.72680026, 0, 1.058078122],
+    [0.22, 0.0053,  17.17547062, 13.64338887, 0, 1.544624544],
+    [0.22, 0.0053,  17.17547062, 12.41305478, 0, 2.317033799],
+    [0.22, 0.0053,  17.17547062, 11.66482931, 0, 2.686043143],
+    [0.22, 0.0053,  17.17547062, 11.09781471, 0, 2.967694189],
+    [0.22, 0.0053,  17.17547062, 9.435795985, 0, 3.913528418],
+    [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825]
+    ],
+    [  # Dataset Ratio r = 0.34
+    [0.34, 0.0063, 17.14061081 , 14.98380876, 0, 0.789078472],
+    [0.34, 0.0063, 17.14061081 , 13.97154915,  0, 1.1299263],
+    [0.34, 0.0063, 17.14061081 , 12.77309253, 0, 1.738136936],
+    [0.34, 0.0063, 17.14061081 , 12.00959929, 0, 2.159520896],
+    [0.34, 0.0063, 17.14061081 , 11.42001731, 0, 2.370047499],
+    [0.34, 0.0063, 17.14061081 , 9.561447179, 0, 3.088299431],
+    [0.34, 0.0063, 17.14061081 , 8.964704969, 0, 3.18097558]
+    ],
+    [  # Dataset Ratio r = 0.43
+    [0.43, 0.0073, 17.07559686 , 15.11316339, 0, 0.577989364],
+    [0.43, 0.0073, 17.07559686 , 14.17997082, 0, 0.829007544],
+    [0.43, 0.0073, 17.07559686 , 13.05739844, 0, 1.094211707],
+    [0.43, 0.0073, 17.07559686 , 12.32309209, 0, 1.325332511],
+    [0.43, 0.0073, 17.07559686 , 11.74608518, 0, 1.400455154],
+    [0.43, 0.0073, 17.07559686 , 9.812372466, 0, 1.832325697],
+    [0.43, 0.0073, 17.07559686 , 9.10519385 , 0, 2.047483977]
+    ],
+    [  # Dataset Ratio r = 0.49
+    [0.49, 0.008,  17.01520754, 15.30614414, 0, 0.357615902],
+    [0.49, 0.008,  17.01520754, 14.49463867, 0, 0.376287785],
+    [0.49, 0.008,  17.01520754, 13.46629742, 0, 0.851008627],
+    [0.49, 0.008,  17.01520754, 12.78388234, 0, 0.904475291],
+    [0.49, 0.008,  17.01520754, 12.23057715, 0, 1.039744708],
+    [0.49, 0.008,  17.01520754, 10.21852839, 0, 1.346405241],
+    [0.49, 0.008,  17.01520754, 9.341730605, 0, 1.566568558]
+    ]
+    ]
+
+    # --- Second Experiment: Rotate "active" bilayer phase 
+    # materialFunctionParameter=[
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/6.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/3.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/2.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, 2.0*(np.pi/3.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, 5.0*(np.pi/6.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, np.pi, 4.262750825]
+    # ]
+
+    # materialFunctionParameter=[
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, 0, 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/12.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/6.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/4.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/3.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, 5.0*(np.pi/12.0), 4.262750825],
+    #     [0.22, 0.0053,  17.17547062, 8.959564147, (np.pi/2.0), 4.262750825]
+    # ]
+
+    # ------ Loops through Parameters for Material Function -----------
+    for i in range(0,np.shape(materialFunctionParameter)[1]):
+        print("------------------")
+        print("New Loop")
+        print("------------------")
+    # Check output directory
+        outputPath = path + str(i)
+        isExist = os.path.exists(outputPath)
+        if not isExist:
+            # Create a new directory because it does not exist
+            os.makedirs(outputPath)
+            print("The new directory " + outputPath + " is created!")
+
+        # thread = threading.Thread(target=run_CellProblem(executable, pythonModule, pythonPath, LOGFILE))
+        # thread.start()
+
+        #TODO: apperently its not possible to pass a variable via subprocess and "calculate" another input value inside the python file.
+        #      Therefore we use this instead.
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_r",materialFunctionParameter[dataset_number][i][0])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_h",materialFunctionParameter[dataset_number][i][1])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_flat",materialFunctionParameter[dataset_number][i][2])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_omega_target",materialFunctionParameter[dataset_number][i][3])
+        SetParameterMaterialFunction(pythonPath + "/" + pythonModule, "param_theta",materialFunctionParameter[dataset_number][i][4])    
+
+        LOGFILE = outputPath + "/" + pythonModule + "_output" + "_" + str(i) + ".log"
+
+        processList = []
+        p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+                                        + " -outputPath " + outputPath
+                                        + " -gamma " + str(gamma) 
+                                        + " | tee " + LOGFILE, shell=True)
+
+        # p = subprocess.Popen(executable + " " + pythonPath + " " + pythonModule
+        #                                 + " -outputPath " + outputPath
+        #                                 + " -gamma " + str(gamma) 
+        #                                 + " -param_r " + str(materialFunctionParameter[i][0])
+        #                                 + " -param_h " + str(materialFunctionParameter[i][1])
+        #                                 + " -param_omega_flat " + str(materialFunctionParameter[i][2])
+        #                                 + " -param_omega_target " + str(materialFunctionParameter[i][3])
+        #                                 + " -phase2_angle " + str(materialFunctionParameter[i][4])
+        #                                 + " | tee " + LOGFILE, shell=True)
+
+        p.wait() # wait
+        processList.append(p)
+        exit_codes = [p.wait() for p in processList]
+        # ---------------------------------------------------
+        # wait here for the result to be available before continuing
+        # thread.join()
+        f = open(outputPath+"/parameter.txt", "w")
+        f.write("r = "+str(materialFunctionParameter[dataset_number][i][0])+"\n")
+        f.write("h = "+str(materialFunctionParameter[dataset_number][i][1])+"\n")
+        f.write("omega_flat = "+str(materialFunctionParameter[dataset_number][i][2])+"\n")        
+        f.write("omega_target = "+str(materialFunctionParameter[dataset_number][i][3])+"\n")         
+        f.close()   
+        #