diff --git a/microstructure_testsuite/microstructure_testsuite.py b/microstructure_testsuite/microstructure_testsuite.py
new file mode 100644
index 0000000000000000000000000000000000000000..4cf31a637d831f249c534c4936bf413f9b1449df
--- /dev/null
+++ b/microstructure_testsuite/microstructure_testsuite.py
@@ -0,0 +1,537 @@
+import subprocess
+import re
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import fileinput
+import time
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import codecs
+import sys
+# import sympy as sym
+# from prettytable import PrettyTable
+# import latextable
+# from texttable import Texttable
+# from tabulate import tabulate
+#-------------------------------------------------------------------------------------------------------
+def run_CellProblem(executable, parset, gridLevel, gamma, mu1,lambda1, rho1, alpha, beta, theta, material_prestrain_imp, outputPath, write_materialFunctions, write_prestrainFunctions, write_VTK,write_L2Error ,write_IntegralMean, write_LOG ):
+ print('----- RUN Cell-Problem ----')
+ processList = []
+ LOGFILE = "Cell-Problem_output.log"
+ print('LOGFILE:',LOGFILE)
+ print('executable:',executable)
+ print('parset:',parset)
+
+ # start_time = time.time()
+ if write_LOG:
+ p = subprocess.Popen(executable + parset
+ # + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
+ + " -gamma " + str(gamma)
+ + " -mu1 " + str(mu1)
+ + " -lambda1 " + str(lambda1)
+ + " -alpha " + str(alpha)
+ + " -beta " + str(beta)
+ + " -theta " + str(theta)
+ + " -material_prestrain_imp " + str(material_prestrain_imp)
+ + " -write_materialFunctions " + str(write_materialFunctions)
+ + " -outputPath " + str(outputPath)
+ + " | tee " + LOGFILE, shell=True)
+
+ else:
+ p = subprocess.Popen(executable + parset
+ # + " -numLevels " + str(gridLevel) + " " + str(gridLevel)
+ + " -gamma " + str(gamma)
+ + " -mu1 " + str(mu1)
+ + " -lambda1 " + str(lambda1)
+ + " -alpha " + str(alpha)
+ + " -beta " + str(beta)
+ + " -theta " + str(theta)
+ + " -material_prestrain_imp " + str(material_prestrain_imp)
+ + " -write_materialFunctions " + str(write_materialFunctions)
+ + " -outputPath " + str(outputPath), stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, shell=True ) # surpress Logging
+
+ p.wait() # wait
+
+ # print("--- %s seconds ---" % (time.time() - start_time))
+ # print('------FINISHED PROGRAM on level:' + str(gridLevel))
+ processList.append(p)
+ ###Wait for all simulation subprocesses before proceeding
+ exit_codes = [p.wait() for p in processList]
+
+ return
+
+#-------------------------------------------------------------------------------------------------------
+def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
+ # Read Output Matrices (effective quantities)
+ # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+ # -- Read Matrix Qhom
+ X = []
+ # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ Q = np.array([[X[0][2], X[1][2], X[2][2]],
+ [X[3][2], X[4][2], X[5][2]],
+ [X[6][2], X[7][2], X[8][2]] ])
+
+ # -- Read Beff (as Vector)
+ X = []
+ # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+ with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+ for line in f:
+ s = line.split()
+ X.append([float(s[i]) for i in range(len(s))])
+ B = np.array([X[0][2], X[1][2], X[2][2]])
+ return Q, B
+
+#-------------------------------------------------------------------------------------------------------
+#-------------------------------------------------------------------------------------------------------
+########################
+#### SET PARAMETERS ####
+########################
+# ----- Setup Paths -----
+outputPath = " ../outputs"
+QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt'
+BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'
+
+print('---- Input parameters: -----')
+mu1 = 10.0
+rho1 = 1.0
+alpha = 2.8
+beta = 2.0
+theta = 1.0/4.0
+gamma = 0.75
+lambda1= 10.0
+
+gridLevel = 4
+#############################################
+# Choose preferred Geometry/Prestrain/Material
+#############################################
+#material_prestrain_imp= "analytical_Example"
+material_prestrain_imp= "parametrized_Laminate"
+
+
+#############################################
+# Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER
+#############################################
+Solvertype = 1 #(default = 1)
+Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+set_IntegralZero = True
+
+#############################################
+# Output-Options
+#############################################
+write_materialFunctions = False
+write_prestrainFunctions = False
+write_VTK = False
+# write_L2Error = True
+# write_IntegralMean = True
+write_L2Error = False
+write_IntegralMean = False
+
+write_LOG = False # writes Cell-Problem output-LOG in "Cell-Problem_output.log"
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('material_prestrain_imp: ', material_prestrain_imp)
+print('gridLevel: ', gridLevel)
+
+parset = ' ../inputs/cellsolver.parset'
+executable = ' ../build-cmake/src/Cell-Problem'
+print('---------------------------------------------------------')
+###########################################################################################################
+
+
+run_CellProblem(executable,
+ parset,
+ gridLevel,
+ gamma,
+ mu1,
+ lambda1,
+ rho1,
+ alpha,
+ beta,
+ theta,
+ material_prestrain_imp,
+ outputPath,
+ write_materialFunctions,
+ write_prestrainFunctions,
+ write_VTK,
+ write_L2Error,
+ write_IntegralMean,
+ write_LOG
+ )
+
+print('Read effective quantities...')
+Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+# Q, B = ReadEffectiveQuantities()
+print('Q:', Q)
+print('B:', B)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#
+#
+# for kappa in [4]:
+#
+# print('--- RUN COMPUTATION FOR KAPPA=',kappa)
+#
+# # --- map R2 into R2
+# if domainDim == 2 and targetDim ==2:
+# print('domainDim == 2 and targetDim == 2')
+#
+# if kappa == 0 or kappa == 1:
+# initialIterate = "initialIterateR2_R2"
+# else:
+# initialIterate = "initialIterateR2_R2_deg"+str(kappa)
+# # initialIterate = "initialIterateR2_R2"
+# # initialIterate = "initialIterateR2_R2_deg2"
+# # initialIterate = "initialIterateR2_R2_deg3"
+# # initialIterate = "initialIterateR2_R2_deg4"
+# # initialIterate = "initialIterateR2_R2_deg5"
+# # initialIterate = "initialIterateR2_R2_deg8"
+#
+# executable = '../build-cmake/src/harmonicmaps_2d'
+# parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
+#
+# # --- map R2 into R3
+# if domainDim == 2 and targetDim ==3:
+# print('domainDim == 2 and targetDim == 3')
+#
+# # initialIterate = "initialIterateR2_R3"
+# # initialIterate = "initialIterateR2_R3_smoothed"
+# if kappa == 0 or kappa == 1:
+# # initialIterate = "initialIterateR2_R3"
+# initialIterate = "initialIterateR2_R3_smoothed"
+# else:
+# initialIterate = "initialIterateR2_R2_deg"+str(kappa)
+#
+# # initialIterate = "initialIterateR2_R3_deg2"
+# # initialIterate = "initialIterateR2_R3_deg3"
+# # initialIterate = "initialIterateR2_R3_deg4"
+# # initialIterate = "initialIterateR2_R3_deg5"
+# executable = '../build-cmake/src/harmonicmaps_2d'
+# parset = ' ../parsets/harmonicmaps_EOC_2d.parset'
+#
+# # --- map R3 into R3
+# if domainDim == 3 and targetDim ==3:
+# print('domainDim == 3 and targetDim == 3')
+#
+# if kappa == 0 or kappa == 1:
+# initialIterate = "initialIterateR3_R3"
+# else:
+# initialIterate = "initialIterateR3_R3_deg"+str(kappa)
+#
+#
+# # initialIterate = "initialIterateR3_R3_pertubed"
+# # initialIterate = "initialIterateR2_R2"
+# # initialIterate = "initialIterateR3_R3_deg2"
+# # initialIterate = "initialIterateR3_R3_deg3"
+# # initialIterate = "initialIterateR3_R3_deg4"
+# # initialIterate = "initialIterateR3_R3_deg5"
+# # initialIterate = "initialIterateR3_R3_deg8"
+# # initialIterate = "initialIterateR3_R3_deg12"
+# # initialIterate = "initialIterateR3_R3_deg3_Pert"
+# # initialIterate = "initialIterateR3_R3_deg3_newOrigin"
+#
+#
+# executable = '../build-cmake/src/harmonicmaps_3d'
+# parset = ' ../parsets/harmonicmaps_EOC_3d.parset'
+#
+#
+#
+#
+#
+#
+#
+# print(' InitialIterate used: ', initialIterate)
+#
+# path = os.getcwd()
+# print("Path: ", path)
+#
+# for order in [1]:
+
+ # minLevel = 5 + 1 - order;
+ # maxLevel = 10 + 1 - order;
+ # minLevel = 2 + 1 - order;
+ # maxLevel = 8 + 1 - order;
+
+ #--- Setup LatexTable
+ # table_1 = Texttable()
+ # # table_1.maketitle = "Test"
+ # table_1.set_cols_align(["l", "r", "c"])
+ # table_1.set_cols_valign(["t", "m", "b"])
+ # table_1.add_rows([["Name", "Age", "Nickname"],
+ # ["Mr\nXavier\nHuon", 32, "Xav'"],
+ # ["Mr\nBaptiste\nClement", 1, "Baby"],
+ # ["Mme\nLouise\nBourgeau", 28, "Lou\n \nLoue"]])
+ #
+ # # table_1.add_rows("\multicolumn{3}{c}{ Projection-based-FE of order: 1, $\kappa = 2$}")
+ # print('-- Example 1: Basic --')
+ # print('Texttable Output:')
+ # print(table_1.draw())
+ # print('\nLatextable Output:')
+ # print(latextable.draw_latex(table_1, caption="An example table.", label="table:example_table"))
+ #
+
+
+
+ # rows = [['Rocket', 'Organisation', 'LEO Payload (Tonnes)', 'Maiden Flight'],
+ # ['Saturn V', 'NASA', '140', '1967'],
+ # ['Space Shuttle', 'NASA', '24.4', '1981'],
+ # ['Falcon 9 FT-Expended', 'SpaceX', '22.8', '2017'],
+ # ['Ariane 5 ECA', 'ESA', '21', '2002']]
+ #
+ # table = Texttable()
+ # table.set_cols_align(["c"] * 4)
+ # table.set_deco(Texttable.HEADER | Texttable.VLINES)
+ # table.add_rows(rows)
+ #
+ # print('Tabulate Table:')
+ # print(tabulate(rows, headers='firstrow'))
+ #
+ # # print('\nTexttable Table:')
+ # # print(table.draw())
+ #
+ # print('\nTabulate Latex:')
+ # print(tabulate(rows, headers='firstrow', tablefmt='latex'))
+
+ # print('\nTexttable Latex:')
+ # print(latextable.draw_latex(table, caption="A comparison of rocket features."))
+
+
+
+#
+# #--- Setup Output-Table:
+# x = PrettyTable()
+# x.title = interpolationMethod + "-FE" + ' of order ' + str(order) + " into R" + str(targetDim)
+# # x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" ]
+# x.field_names = ["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$", "$\lnorm{u_h-u_{2h}$", "$\lnorm{u_{2h}-u_{4h}$", "$\seminorm{u_h-u_{2h}}$", "$\seminorm{u_{2h}-u_{4h}}$" , "wall-time" ]
+# x.align["k"] = "l" # Left align city names
+# x.padding_width = 1 # One space between column edges and contents (default)
+#
+# rows = []
+# rows.append(["k", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" , "wall-time" ])
+# # rows.append(["k", "#Triang/#DOF", "# TR-Steps","wall-time","$E^{Di}(0)$", "$E^{Di}$" , "$EOC_h^{L^2}$", "$EOC_h^{H^1}$", "$EOC_h^{E}$" ])
+# print('rows:', rows)
+#
+# if compute_HarmonicMaps:
+#
+# computeHarmonicMap(minLevel,maxLevel,domainDim,targetDim,order,interpolationMethod,maxTrustRegionSteps,initialIterate,randomInitialIterate,readConfiguration,configurationFile,pertubation,pertubationRadius,epsilon,tolerance,executable,parset )
+#
+# # --------------------
+#
+#
+# if compute_Error:
+# print("-------------------------------------------------------")
+# print("Now measuring errors with discretizationErrorMode:" + str(discretizationErrorMode))
+# # subprocess.call(["echo", "Now measuring errors"])
+#
+# # EOC_l2_list = ["-","-"]
+#
+# energyList = []
+# initial_energyList = []
+# stepList = []
+# ndofList = []
+# timeList = []
+# # EOC_l2_list = []
+#
+# #
+# for numLevels in range(minLevel,maxLevel+1):
+#
+# print("NUMBER OF LEVELS:", numLevels)
+# # Measure the discretization errors against the solution on the finest grid
+# # LOGFILE = "./compute-disc-error_" + str(order) + "_" + str(numLevels) + ".log"
+#
+#
+# #subprocess.Popen(["../build-cmake/src/compute-disc-error", "compute-disc-error-skyrmions-hexagon.parset",
+# #"-order", str(order),
+# #"-numLevels", str(numLevels),
+# #"-numReferenceLevels", str(maxLevel),
+# #"-simulationData", "harmonicmaps-result-" + str(order) + "-" + str(numLevels) + ".data",
+# #"-referenceData", "harmonicmaps-result-" + str(order) + "-" + str(maxLevel) + ".data"])
+#
+#
+# # ------------ Get Energy / nDofs/Steps etc. -----------------------------------------------------------------------------
+# LOGFILE_comp = "./harmonicmapsR"+ str(domainDim) +"_intoR"+ str(targetDim) + "_deg" + str(kappa) + "_"+ interpolationMethod + "_" + str(order) + "_" + str(numLevels) + ".log"
+# # Read Energy Values:
+# with open(LOGFILE_comp, 'r') as file:
+# output = file.read()
+#
+# try:
+# # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1) # muss nun nichtmehr am Zeilenanfang stehen! :)
+# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1)
+# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
+# # tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
+# tmp_energy = re.findall(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
+# tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
+# tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
+# # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
+# # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
+# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
+#
+# except AttributeError:
+# # tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output) # muss nun nichtmehr am Zeilenanfang stehen! :)
+# tmp_energy = re.search(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
+# tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
+# tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
+# # tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
+# # tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
+# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+\d?)',output)
+#
+# if maxTrustRegionSteps>0 :
+# print('tmp_step:',tmp_step)
+# print('tmp_step last:', int(tmp_step[-1]))
+#
+# print('tmp_dofs:',tmp_dofs)
+# print('tmp_dofs last:', int(tmp_dofs[-1]))
+#
+# print('type tmp_energy:', type(tmp_energy))
+# print('tmp_energy:', tmp_energy)
+# print('tmp_energy last:', tmp_energy[-1])
+#
+# print('tmp_time:', tmp_time)
+#
+# energy = float(tmp_energy[-1]) #[1] since otherwise it recognizes "2" from L2error...
+# energyList.append(energy)
+#
+# initialEnergy = float(tmp_energy[0])
+# initial_energyList.append(initialEnergy)
+#
+# steps = int(tmp_step[-1])+1 #starts with zero therefore +1
+# stepList.append(steps)
+#
+# ndofs = int(tmp_dofs[-1])
+# ndofList.append(ndofs)
+#
+# time = float(tmp_time[-1])
+# timeList.append(time)
+# # -----------------------------------------------------------------------------------------------
+#
+# if numLevels >= (minLevel+2) and discretizationErrorMode=='discrete':
+#
+# [EOC_l2, EOC_h1, EOC_energy, L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
+# = computeError_discrete(numLevels,minLevel,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList)
+#
+# # x.add_row([numLevels, currentDofs , currentSteps,currentInitialEnergy, currentEnergy, EOC_l2, EOC_h1, EOC_energy , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
+# # x.add_row([numLevels, "-" , "-", "energy", EOC_l2, EOC_h1, "EOC_h^{E}" , format(l2_fine, "5.3e"), format(l2_coarse, "5.3e"), format(h1_fine, "5.3e"), format(h1_coarse, "5.3e") ])
+#
+# elif numLevels >= (minLevel+1) and discretizationErrorMode=='analytical':
+#
+# [EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse, constraintError] \
+# = computeError_analytical(numLevels,domainDim,targetDim,order,interpolationMethod,discretizationErrorMode,maxTrustRegionSteps,energyList,ndofList,stepList,initial_energyList,referenceSolution,numReferenceLevels)
+#
+# else:
+# EOC_l2 = "-"
+# EOC_h1 = "-"
+# EOC_energy = "-"
+# L2_fine = "-"
+# L2_coarse = "-"
+# H1_fine = "-"
+# H1_coarse = "-"
+# constraintError = "-"
+#
+# # if maxTrustRegionSteps > 0 :
+# # print('energyList', energyList)
+# # # currentEnergy = energyList[numLevels-1]
+# # currentEnergy = energyList[-1] # better like this?
+# # print('currentEnergy', currentEnergy)
+# #
+# # # currentDofs = ndofList[numLevels-1]
+# # currentDofs = ndofList[-1]
+# # print('currentDofs', currentDofs)
+# #
+# # # currentSteps = stepList[numLevels-1]
+# # currentSteps = stepList[-1]
+# # print('currentSteps', currentSteps)
+# #
+# # print('initial_energyList', initial_energyList)
+# # currentInitialEnergy = initial_energyList[-1]
+# # print('currentInitialEnergy', currentInitialEnergy)
+# #
+# # else:
+# # currentEnergy = "-"
+# # currentDofs = "-"
+# # currentSteps = "-"
+# # currentInitialEnergy = "-"
+#
+# # x.add_row([numLevels, currentDofs , currentSteps, currentInitialEnergy, currentEnergy, "-", "-" , "-" , "-", "-", "-", "-"])
+# # print('Energy List:', energyList)
+#
+#
+# if maxTrustRegionSteps > 0 :
+# print('energyList', energyList)
+# # currentEnergy = energyList[numLevels-1]
+# currentEnergy = energyList[-1] # better like this?
+# print('currentEnergy', currentEnergy)
+#
+# # currentDofs = ndofList[numLevels-1]
+# currentDofs = ndofList[-1]
+# print('currentDofs', currentDofs)
+#
+# # currentSteps = stepList[numLevels-1]
+# currentSteps = stepList[-1]
+# print('currentSteps', currentSteps)
+#
+# print('initial_energyList', initial_energyList)
+# currentInitialEnergy = initial_energyList[-1]
+# print('currentInitialEnergy', currentInitialEnergy)
+#
+#
+#
+#
+# else:
+# currentEnergy = "-"
+# currentDofs = "-"
+# currentSteps = "-"
+# currentInitialEnergy = "-"
+#
+# print('timeList:', timeList)
+# currentTime = timeList[-1]
+# print('currentSteps:', currentSteps)
+#
+# x.add_row([numLevels, currentDofs, currentSteps, currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy,L2_fine , L2_coarse , H1_fine , H1_coarse , currentTime])
+# rows.append([numLevels, currentDofs, str(currentSteps), currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy , currentTime])
+#
+# ## Add extra column:
+# # print(*EOC_l2_list)
+# # x.add_column('EOC', EOC_l2_list)
+#
+# # Write Table to Text-File:
+# tablefile = open("EOC-table_R"+str(domainDim)+"intoR"+ str(targetDim)+ "_deg" + str(kappa) + "_" + interpolationMethod+ "_order" +str(order) + "tol" + str(tolerance) + ".txt", "w")
+# tablefile.write(str(x))
+# tablefile.write('\n')
+# tablefile.write(str(tabulate(rows, headers='firstrow', tablefmt='latex')))
+#
+# tablefile.close()
+# # print Table
+# print(x)
+#
+#
+#
+# #--- Try Tabulate:
+#
+# print(tabulate(rows, headers='firstrow', tablefmt='latex'))
+# ########## end of kappa-loop ##########################