diff --git a/src/Plot-Angle-GammaV2.py b/src/Plot-Angle-GammaV2.py
new file mode 100644
index 0000000000000000000000000000000000000000..7982203914e20664c325ff1a58ba07ab397a69bf
--- /dev/null
+++ b/src/Plot-Angle-GammaV2.py
@@ -0,0 +1,368 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+import matplotlib.ticker as ticker
+
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# from matplotlib import rc
+# rc('text', usetex=True) # Use LaTeX font
+#
+# import seaborn as sns
+# sns.set(color_codes=True)
+
+
+def format_func(value, tick_number):
+ # # find number of multiples of pi/2
+ # N = int(np.round(2 * value / np.pi))
+ # if N == 0:
+ # return "0"
+ # elif N == 1:
+ # return r"$\pi/2$"
+ # elif N == 2:
+ # return r"$\pi$"
+ # elif N % 2 > 0:
+ # return r"${0}\pi/2$".format(N)
+ # else:
+ # return r"${0}\pi$".format(N // 2)
+ # find number of multiples of pi/2
+ N = int(np.round(4 * value / np.pi))
+ if N == 0:
+ return "0"
+ elif N == 1:
+ return r"$\pi/4$"
+ elif N == 2:
+ return r"$\pi/2$"
+ elif N % 2 > 0:
+ return r"${0}\pi/2$".format(N)
+ else:
+ return r"${0}\pi$".format(N // 2)
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+print('---- Input parameters: -----')
+alpha = 10
+mu1 = 1.0
+rho1 = 1.0
+beta = 2.0 #5.0
+theta = 1.0/8.0
+#
+
+alpha = -0.5
+beta = 40.0
+theta= 1/8.0
+
+
+
+# # INTERESTING! from pi/2:
+alpha = -0.5
+beta = 40.0
+theta= 1/8.0
+#
+# # # INTERESTING! from pi/2:
+# alpha = -0.2
+# beta = 25.0
+# theta= 1/2
+
+# INTERESTING!:
+# alpha = -0.5
+# beta = 5.0
+# theta= 1/30
+
+
+
+# INTERESTING!:
+# alpha = -0.25
+# beta = 10.0
+# theta= 3/4
+
+
+# # INTERESTING!:
+alpha = -0.25
+beta = 10.0
+theta= 1/8
+
+#
+# INTERESTING!:
+# alpha = -0.25
+# beta = 5.0
+# theta= 1/8
+#
+
+
+# # INTERESTING!:
+alpha = -0.5
+beta = 10.0
+theta= 1/8
+
+
+
+
+
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 1.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=100) # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+# Gamma_Values = '0'
+
+
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+
+print('angles:', angles)
+
+
+
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+determinantVec = np.vectorize(determinant)
+
+detValues = determinantVec(q1,q2,muGammas,q12)
+print('detValues:', detValues)
+
+
+detZeroidx = find_nearestIdx(detValues, 0)
+print('idx where det nearest to zero', idx)
+gammaClose = Gamma_Values[detZeroidx]
+print('gammaClose:', gammaClose)
+
+
+# --- Convert to numpy array
+Gamma_Values = np.array(Gamma_Values)
+angles = np.array(angles)
+
+# ---------------- Create Plot -------------------
+# plt.figure()
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+width = 6.28 *0.5
+height = width / 1.618
+fig = plt.figure()
+# ax = plt.axes((0.15,0.21 ,0.75,0.75))
+ax = plt.axes((0.15,0.21 ,0.8,0.75))
+ax.tick_params(axis='x',which='major', direction='out',pad=5)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax.grid(True,which='major',axis='both',alpha=0.3)
+
+
+# # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# # plt.title(r'angle$-\gamma$-Plot')
+# plt.plot(Gamma_Values, angles)
+# plt.scatter(Gamma_Values, angles)
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+
+# f,ax=plt.subplots(1)
+
+
+# ax.plot(muGammas, angles)
+# ax.scatter(muGammas, angles)
+
+ax.plot(Gamma_Values, angles, 'royalblue', zorder=3, )
+# ax.scatter(Gamma_Values, angles)
+
+# ax.set_xlabel(r"$q_3(\gamma)$")
+ax.set_xlabel(r"$\gamma$")
+# ax.set_ylabel(r"angle $\angle$")
+ax.set_ylabel(r"angle $\alpha$")
+
+# plt.xlabel("$q_3$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# ax.grid(True)
+
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25))
+
+
+# ax.yaxis.set_major_formatter(ticker.FuncFormatter(
+# lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi))
+
+# ---------------------------- show pi values ------------------------------------
+# ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# ax.legend()
+# # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# # ax.set_yticks([0, np.pi/4 ,np.pi/2])
+# # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+# ---------------------------------------------------------------
+
+
+# ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set_yticks([0, np.pi/4 ,np.pi/2])
+# labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+
+
+# OLD :
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+
+ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
+labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']
+ax.set_yticklabels(labels)
+
+
+# Plot Gamma Value that is closest to q3_star
+ax.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', linewidth=1, label='$\gamma^*$')
+
+
+# color elliptic/hyperbolic region
+# ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2)
+# ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2)
+
+
+ax.legend(loc='upper right')
+
+# plt.xlabel("$q_3(\gamma)$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# plt.legend(loc='upper center')
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Angle-Gamma.pdf')
+
+plt.show()
+
+
+
+
+# plt.figure()
+# plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.xlabel("$\mu_\gamma$")
+# plt.ylabel("angle")
+# plt.legend()
+# plt.show()
+#
diff --git a/src/Plot-Curvature-GammaV2.py b/src/Plot-Curvature-GammaV2.py
new file mode 100644
index 0000000000000000000000000000000000000000..4726a9133690a9390b06310685cf4914288a01d9
--- /dev/null
+++ b/src/Plot-Curvature-GammaV2.py
@@ -0,0 +1,337 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+import matplotlib.ticker as ticker
+
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# from matplotlib import rc
+# rc('text', usetex=True) # Use LaTeX font
+#
+# import seaborn as sns
+# sns.set(color_codes=True)
+
+
+def format_func(value, tick_number):
+ # # find number of multiples of pi/2
+ # N = int(np.round(2 * value / np.pi))
+ # if N == 0:
+ # return "0"
+ # elif N == 1:
+ # return r"$\pi/2$"
+ # elif N == 2:
+ # return r"$\pi$"
+ # elif N % 2 > 0:
+ # return r"${0}\pi/2$".format(N)
+ # else:
+ # return r"${0}\pi$".format(N // 2)
+ # find number of multiples of pi/2
+ N = int(np.round(4 * value / np.pi))
+ if N == 0:
+ return "0"
+ elif N == 1:
+ return r"$\pi/4$"
+ elif N == 2:
+ return r"$\pi/2$"
+ elif N % 2 > 0:
+ return r"${0}\pi/2$".format(N)
+ else:
+ return r"${0}\pi$".format(N // 2)
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+print('---- Input parameters: -----')
+
+alpha = 10.0
+mu1 = 1.0
+rho1 = 1.0
+beta = 2.0 #5.0
+theta = 1.0/8.0
+#
+
+
+# # INTERESTING!:
+alpha = -0.25
+beta = 10.0
+theta= 1/8
+
+
+# INTERESTING!:
+# alpha = -0.5
+# beta = 10.0
+# theta= 1/8
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 1.0
+
+gamma_max = 3.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=50) # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+# Gamma_Values = '0'
+
+
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas, mu1, rho1)
+
+print('angles:', angles)
+
+print('curvature:', curvature)
+
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+determinantVec = np.vectorize(determinant)
+
+detValues = determinantVec(q1,q2,muGammas,q12)
+print('detValues:', detValues)
+
+
+detZeroidx = find_nearestIdx(detValues, 0)
+print('idx where det nearest to zero', idx)
+gammaClose = Gamma_Values[detZeroidx]
+print('gammaClose:', gammaClose)
+
+
+# --- Convert to numpy array
+Gamma_Values = np.array(Gamma_Values)
+angles = np.array(angles)
+
+
+curvature = np.array(curvature)
+
+# ---------------- Create Plot -------------------
+# plt.figure()
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+width = 6.28 *0.5
+height = width / 1.618
+fig = plt.figure()
+# ax = plt.axes((0.15,0.21 ,0.75,0.75))
+ax = plt.axes((0.15,0.18 ,0.8,0.75))
+# ax = plt.axes((0.21,0.21 ,0.8,0.75))
+ax.tick_params(axis='x',which='major', direction='out',pad=1)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=1) # changed pad = distance to title to 1 here!
+# ax.xaxis.set_major_locator(MultipleLocator(0.1))
+# ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+
+
+ax.xaxis.set_major_locator(MultipleLocator(0.3))
+ax.xaxis.set_minor_locator(MultipleLocator(0.15))
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax.grid(True,which='major',axis='both',alpha=0.3)
+
+
+# # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# # plt.title(r'angle$-\gamma$-Plot')
+# plt.plot(Gamma_Values, angles)
+# plt.scatter(Gamma_Values, angles)
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+
+# f,ax=plt.subplots(1)
+
+
+# ax.plot(muGammas, angles)
+# ax.scatter(muGammas, angles
+#
+# ax.plot(Gamma_Values, angles, 'royalblue', zorder=3, )
+ax.plot(Gamma_Values, curvature, 'royalblue', zorder=3, )
+
+# ax.scatter(Gamma_Values, angles)
+
+# ax.set_xlabel(r"$q_3(\gamma)$")
+ax.set_xlabel(r"$\gamma$")
+# ax.set_ylabel(r"curvature $\kappa$")
+ax.set_title(r"curvature $\kappa$", fontsize=9, pad = 4)
+
+# plt.xlabel("$q_3$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# ax.grid(True)
+
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25))
+
+
+# ax.yaxis.set_major_formatter(ticker.FuncFormatter(
+# lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi))
+
+# ---------------------------- show pi values ------------------------------------
+# ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# ax.legend()
+# # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# # ax.set_yticks([0, np.pi/4 ,np.pi/2])
+# # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+# ---------------------------------------------------------------
+
+
+# ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set_yticks([0, np.pi/4 ,np.pi/2])
+# labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+
+
+# OLD :
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+
+
+# Plot Gamma Value that is closest to q3_star
+ax.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed',linewidth=1, label='$\gamma^*$')
+
+
+
+# ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2)
+# ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2)
+
+
+ax.legend(loc='lower right')
+
+# plt.xlabel("$q_3(\gamma)$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# plt.legend(loc='upper center')
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-Gamma.pdf')
+
+plt.show()
+
+
+
+
+# plt.figure()
+# plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.xlabel("$\mu_\gamma$")
+# plt.ylabel("angle")
+# plt.legend()
+# plt.show()
+#
diff --git a/src/Plot_Angle_Lemma1.4V2.py b/src/Plot_Angle_Lemma1.4V2.py
new file mode 100644
index 0000000000000000000000000000000000000000..842880f5263a24f5ff7df10f0eb488291fe6d1d1
--- /dev/null
+++ b/src/Plot_Angle_Lemma1.4V2.py
@@ -0,0 +1,666 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+from HelperFunctions import *
+from ClassifyMin import *
+
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# import tikzplotlib
+# # from pylab import *
+# from tikzplotlib import save as tikz_save
+
+
+# Needed ?
+mpl.use('pdf')
+
+# from subprocess import Popen, PIPE
+#import sys
+
+###################### makePlot.py #########################
+# Generalized Plot-Script giving the option to define
+# quantity of interest and the parameter it depends on
+# to create a plot
+#
+# Input: Define y & x for "x-y plot" as Strings
+# - Run the 'Cell-Problem' for the different Parameter-Points
+# (alternatively run 'Compute_MuGamma' if quantity of interest
+# is q3=muGamma for a significant Speedup)
+
+###########################################################
+
+
+
+# figsize argument takes inputs in inches
+# and we have the width of our document in pts.
+# To set the figure size we construct a function
+# to convert from pts to inches and to determine
+# an aesthetic figure height using the golden ratio:
+# def set_size(width, fraction=1):
+# """Set figure dimensions to avoid scaling in LaTeX.
+#
+# Parameters
+# ----------
+# width: float
+# Document textwidth or columnwidth in pts
+# fraction: float, optional
+# Fraction of the width which you wish the figure to occupy
+#
+# Returns
+# -------
+# fig_dim: tuple
+# Dimensions of figure in inches
+# """
+# # Width of figure (in pts)
+# fig_width_pt = width * fraction
+#
+# # Convert from pt to inches
+# inches_per_pt = 1 / 72.27
+#
+# # Golden ratio to set aesthetic figure height
+# # https://disq.us/p/2940ij3
+# golden_ratio = (5**.5 - 1) / 2
+#
+# # Figure width in inches
+# fig_width_in = fig_width_pt * inches_per_pt
+# # Figure height in inches
+# fig_height_in = fig_width_in * golden_ratio
+#
+# fig_dim = (fig_width_in, fig_height_in)
+#
+# return fig_dim
+#
+
+
+
+def format_func(value, tick_number):
+ # # find number of multiples of pi/2
+ # N = int(np.round(2 * value / np.pi))
+ # if N == 0:
+ # return "0"
+ # elif N == 1:
+ # return r"$\pi/2$"
+ # elif N == 2:
+ # return r"$\pi$"
+ # elif N % 2 > 0:
+ # return r"${0}\pi/2$".format(N)
+ # else:
+ # return r"${0}\pi$".format(N // 2)
+ # find number of multiples of pi/2
+ N = int(np.round(4 * value / np.pi))
+ if N == 0:
+ return "0"
+ elif N == 1:
+ return r"$\pi/4$"
+ elif N == 2:
+ return r"$\pi/2$"
+ elif N % 2 > 0:
+ return r"${0}\pi/2$".format(N)
+ else:
+ return r"${0}\pi$".format(N // 2)
+
+
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+
+# TODO
+# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
+# - Also Add option to plot Minimization Output
+
+
+# ----- Setup Paths -----
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+#---------------------------------------------------------------
+
+print('---- Input parameters: -----')
+mu1 = 1.0 #10.0
+# lambda1 = 10.0
+rho1 = 1.0
+alpha = 5.0
+beta = 10.0
+# alpha = 2.0
+# beta = 2.0
+theta = 1.0/8.0 #1.0/4.0
+
+lambda1 = 0.0
+# gamma = 1.0/4.0
+
+# TEST:
+alpha=3.0;
+
+
+
+
+# # INTERESTING!:
+alpha = 5.0
+beta = 10.0
+theta= 1/8
+
+
+
+alpha_Values
+
+
+
+gamma = 'infinity' #Elliptic Setting
+# gamma = '0' #Hyperbolic Setting
+# gamma = 0.5
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+
+# --- define Interval of x-va1ues:
+xmin = 0.01
+xmax = 0.41
+xmax = 0.99
+
+
+Jumps = False
+
+
+numPoints = 2000
+numPoints = 1000
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+
+
+
+
+for theta in X_Values:
+
+ print('Situation of Lemma1.4')
+ q12 = 0.0
+ q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+ q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+ b1 = prestrain_b1(rho1, beta, alpha,theta)
+ b2 = prestrain_b2(rho1, beta, alpha,theta)
+ b3 = 0.0
+
+ q3 = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
+
+
+
+ G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3, mu1, rho1)
+
+ Y_Values.append(angle)
+
+
+
+print("(Output) Values of " + yName + ": ", Y_Values)
+
+
+idx = find_nearestIdx(Y_Values, 0)
+print(' Idx of value closest to 0', idx)
+ValueClose = Y_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', ValueClose)
+
+
+
+# Find Indices where the difference between the next one is larger than epsilon...
+jump_idx = []
+jump_xValues = []
+jump_yValues = []
+tmp = X_Values[0]
+for idx, x in enumerate(X_Values):
+ print(idx, x)
+ if idx > 0:
+ if abs(Y_Values[idx]-Y_Values[idx-1]) > 1:
+ print('jump candidate')
+ jump_idx.append(idx)
+ jump_xValues.append(x)
+ jump_yValues.append(Y_Values[idx])
+
+
+
+
+
+
+
+print("Jump Indices", jump_idx)
+print("Jump X-values:", jump_xValues)
+print("Jump Y-values:", jump_yValues)
+
+y_plotValues = [Y_Values[0]]
+x_plotValues = [X_Values[0]]
+# y_plotValues.extend(jump_yValues)
+for i in jump_idx:
+ y_plotValues.extend([Y_Values[i-1], Y_Values[i]])
+ x_plotValues.extend([X_Values[i-1], X_Values[i]])
+
+
+y_plotValues.append(Y_Values[-1])
+# x_plotValues = [X_Values[0]]
+# x_plotValues.extend(jump_xValues)
+x_plotValues.append(X_Values[-1])
+
+
+print("y_plotValues:", y_plotValues)
+print("x_plotValues:", x_plotValues)
+# Y_Values[np.diff(y) >= 0.5] = np.nan
+
+
+#get values bigger than jump position
+# gamma = infty
+# x_rest = X_Values[X_Values>x_plotValues[1]]
+# Y_Values = np.array(Y_Values) #convert the np array
+# y_rest = Y_Values[X_Values>x_plotValues[1]]
+#
+#
+# # gamma = 0
+# x_rest = X_Values[X_Values>x_plotValues[3]]
+# Y_Values = np.array(Y_Values) #convert the np array
+# y_rest = Y_Values[X_Values>x_plotValues[3]]
+
+# gamma between
+# Y_Values = np.array(Y_Values) #convert the np array
+# X_Values = np.array(X_Values) #convert the np array
+#
+# x_one = X_Values[X_Values>x_plotValues[3]]
+# # ax.scatter(X_Values, Y_Values)
+# y_rest = Y_Values[X_Values>x_plotValues[3]]
+# ax.plot(X_Values[X_Values>0.135], Y_Values[X_Values<0.135])
+#
+#
+#
+
+
+# y_rest = Y_Values[np.nonzero(X_Values>x_plotValues[1]]
+# print('X_Values:', X_Values)
+# print('Y_Values:', Y_Values)
+# print('x_rest:', x_rest)
+# print('y_rest:', y_rest)
+# print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) )
+
+
+
+
+# --- Convert to numpy array
+Y_Values = np.array(Y_Values)
+X_Values = np.array(X_Values)
+
+# ---------------- Create Plot -------------------
+
+#--- change plot style: SEABORN
+# plt.style.use("seaborn-paper")
+
+
+#--- Adjust gobal matplotlib variables
+# mpl.rcParams['pdf.fonttype'] = 42
+# mpl.rcParams['ps.fonttype'] = 42
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+
+
+# plt.rc('font', family='serif', serif='Times')
+# plt.rc('font', family='serif')
+# # plt.rc('text', usetex=True) #also works...
+# plt.rc('xtick', labelsize=8)
+# plt.rc('ytick', labelsize=8)
+# plt.rc('axes', labelsize=8)
+
+
+
+
+
+#---- Scale Figure apropriately to fit tex-File Width
+# width = 452.9679
+
+# width as measured in inkscape
+width = 6.28 *0.5
+height = width / 1.618
+
+#setup canvas first
+fig = plt.figure() #main
+# fig, ax = plt.subplots()
+# fig, (ax, ax2) = plt.subplots(ncols=2)
+# fig,axes = plt.subplots(nrows=1,ncols=2,figsize=(width,height)) # more than one plot
+
+
+# fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97) #TEST
+
+
+# TEST
+# mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1)
+# fig = plt.figure(figsize=(width+0.1,height+0.1))
+
+
+# mpl.rcParams['figure.figsize'] = (width,height)
+# fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=(width,height)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=set_size(width))
+# fig = plt.subplots(1, 1, figsize=set_size(width))
+
+# --- To create a figure half the width of your document:#
+# fig = plt.figure(figsize=set_size(width, fraction=0.5))
+
+
+
+#--- You must select the correct size of the plot in advance
+# fig.set_size_inches(3.54,3.54)
+
+ax = plt.axes((0.15,0.18,0.8,0.8))
+# ax = plt.axes((0.1,0.1,0.5,0.8))
+# ax = plt.axes((0.1,0.1,1,1))
+# ax = plt.axes()
+
+# ax.spines['right'].set_visible(False)
+# ax.spines['left'].set_visible(False)
+# ax.spines['bottom'].set_visible(False)
+# ax.spines['top'].set_visible(False)
+# ax.tick_params(axis='x',which='major',direction='out',length=10,width=5,color='red',pad=15,labelsize=15,labelcolor='green',
+# labelrotation=15)
+# ax.tick_params(axis='x',which='major', direction='out',pad=5,labelsize=10)
+# ax.tick_params(axis='y',which='major', length=5, width=1, direction='out',pad=5,labelsize=10)
+ax.tick_params(axis='x',which='major', direction='out',pad=3)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+ax.xaxis.set_major_locator(MultipleLocator(0.05))
+ax.xaxis.set_minor_locator(MultipleLocator(0.025))
+
+
+#---- print data-types
+print(ax.xaxis.get_major_locator())
+print(ax.xaxis.get_minor_locator())
+print(ax.xaxis.get_major_formatter())
+print(ax.xaxis.get_minor_formatter())
+
+#---- Hide Ticks or Labels
+# ax.yaxis.set_major_locator(plt.NullLocator())
+# ax.xaxis.set_major_formatter(plt.NullFormatter())
+
+#---- Reducing or Increasing the Number of Ticks
+# ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+# ax.yaxis.set_major_locator(plt.MaxNLocator(3))
+
+
+#----- Fancy Tick Formats
+ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+
+
+
+
+
+
+# --- manually change ticks&labels:
+# ax.set_xticks([0.2,1])
+# ax.set_xticklabels(['pos1','pos2'])
+
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+
+a=ax.yaxis.get_major_locator()
+b=ax.yaxis.get_major_formatter()
+c = ax.get_xticks()
+d = ax.get_xticklabels()
+print('xticks:',c)
+print('xticklabels:',d)
+
+ax.grid(True,which='major',axis='both',alpha=0.3)
+
+
+
+
+
+
+# plt.figure()
+
+# f,ax=plt.subplots(1)
+
+# plt.title(r''+ yName + '-Plot')
+# plt.plot(X_Values, Y_Values,linewidth=2, '.k')
+# plt.plot(X_Values, Y_Values,'.k',markersize=1)
+# plt.plot(X_Values, Y_Values,'.',markersize=0.8)
+
+# plt.plot(X_Values, Y_Values)
+
+# ax.plot([[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])
+
+
+
+# Gamma = '0'
+# ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b')
+#
+# ax.plot([x_plotValues[1],x_plotValues[3]], [y_plotValues[2],y_plotValues[3]] , 'b')
+#
+# ax.plot(x_rest, y_rest, 'b')
+
+
+# Gamma between
+
+# x jump values (gamma 0): [0.13606060606060608, 0.21090909090909093]
+
+# ax.plot([[0,jump_xValues[0]], [0, 0]] , 'b')
+# ax.plot([jump_xValues[0],xmin], [y_plotValues[2],y_plotValues[2]] , 'b')
+
+# ax.plot([[0,0.13606060606060608], [0, 0]] , 'b')
+# ax.plot([[0.13606060606060608,xmin], [(math.pi/2),(math.pi/2)]], 'b')
+
+# jump_xValues[0]
+
+
+
+# --- leave out jumps:
+# ax.scatter(X_Values, Y_Values)
+
+ax.set_xlabel(r"volume fraction $\theta$")
+ax.set_ylabel(r"angle $\alpha$")
+
+
+if Jumps:
+
+ # --- leave out jumps:
+ if gamma == 'infinity':
+ ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]] , 'royalblue')
+ ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue')
+
+
+
+ # ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]])
+ # ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]])
+
+
+
+
+ # ax.plot(X_Values[X_Values>0.136], Y_Values[X_Values>0.136])
+ # ax.plot(X_Values[X_Values<0.135], Y_Values[X_Values<0.135])
+ # ax.scatter(X_Values, Y_Values)
+ # ax.plot(X_Values, Y_Values)
+
+ # plt.plot(x_plotValues, y_plotValues,'.')
+ # plt.scatter(X_Values, Y_Values, alpha=0.3)
+ # plt.scatter(X_Values, Y_Values)
+ # plt.plot(X_Values, Y_Values,'.')
+ # plt.plot([X_Values[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])
+ # plt.axis([0, 6, 0, 20])
+
+ # ax.set_xlabel(r"volume fraction $\theta$", size=11)
+ # ax.set_ylabel(r"angle $\angle$", size=11)
+ # ax.set_xlabel(r"volume fraction $\theta$")
+ # # ax.set_ylabel(r"angle $\angle$")
+ # ax.set_ylabel(r"angle $\alpha$")
+ # plt.ylabel('$\kappa$')
+
+ # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+ # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
+
+
+
+
+ # Plot every other line.. not the jumps...
+
+ if gamma == '0':
+ tmp = 1
+ for idx, x in enumerate(x_plotValues):
+ if idx > 0 and tmp == 1:
+ # plt.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]] )
+ ax.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]], 'royalblue', zorder=2)
+ tmp = 0
+ else:
+ tmp = 1
+
+ # plt.plot([x_plotValues[0],x_plotValues[1]] ,[y_plotValues[0],y_plotValues[1]] )
+ # plt.plot([x_plotValues[2],x_plotValues[3]] ,[y_plotValues[2],y_plotValues[3]] )
+ # plt.plot([x_plotValues[4],x_plotValues[5]] ,[y_plotValues[4],y_plotValues[5]] )
+ # plt.plot([x_plotValues[6],x_plotValues[7]] ,[y_plotValues[6],y_plotValues[7]] )
+
+
+ for x in jump_xValues:
+ plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1, zorder=1)
+ # plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', label=r'$\theta_*$')
+
+ # plt.axvline(x_plotValues[1],ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')
+
+ # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+ # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+ # plt.legend()
+
+
+ # -- SETUP LEGEND
+ # ax.legend(prop={'size': 11})
+ # ax.legend()
+
+ # ------------------ SAVE FIGURE
+ # tikzplotlib.save("TesTout.tex")
+ # plt.close()
+ # mpl.rcParams.update(mpl.rcParamsDefault)
+
+ # plt.savefig("graph.pdf",
+ # #This is simple recomendation for publication plots
+ # dpi=1000,
+ # # Plot will be occupy a maximum of available space
+ # bbox_inches='tight',
+ # )
+ # plt.savefig("graph.pdf")
+
+
+
+ # ---- ADD additional scatter:
+ # ax.scatter(X_Values,Y_Values,s=1,c='black',zorder=4)
+
+ # Find transition point
+ lastIdx = len(Y_Values)-1
+
+ for idx, y in enumerate(Y_Values):
+ if idx != lastIdx:
+ if abs(y-0) < 0.01 and abs(Y_Values[idx+1] - 0) > 0.05:
+ transition_point1 = X_Values[idx+1]
+ print('transition point1:', transition_point1 )
+ if abs(y-0.5*np.pi) < 0.01 and abs(Y_Values[idx+1] -0.5*np.pi)>0.01:
+ transition_point2 = X_Values[idx]
+ print('transition point2:', transition_point2 )
+ if abs(y-0) > 0.01 and abs(Y_Values[idx+1] - 0) < 0.01:
+ transition_point3 = X_Values[idx+1]
+ print('transition point3:', transition_point3 )
+
+ # Add transition Points:
+ if gamma == '0':
+ ax.scatter([transition_point1, transition_point2],[np.pi/2,np.pi/2],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2+0.012 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+ else:
+ ax.scatter([transition_point1, transition_point2, transition_point3 ],[np.pi/2,np.pi/2,0 ],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2 +0.011 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point3 +0.009 , 0+0.08, r"$3$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+else:
+ ax.scatter(X_Values,Y_Values,s=1, marker='o', cmap=None, norm=None, facecolor = 'blue',
+ edgecolor = 'none', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+ ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8 ])
+ labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$',r'$5\pi/8$']
+ ax.set_yticklabels(labels)
+ # ax.set_yticks([1.570786327, np.pi/2 ])
+ # labels = [r'$\pi/2-0.0005 $' , r'$\pi/2$']
+ # ax.set_yticklabels(labels)
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Angle-Theta.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+# figureheight = '\\figureheight',
+# figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
+# #---------------------------------------------------------------
diff --git a/src/Plot_Angle_Lemma1.4_ChangeGamma.py b/src/Plot_Angle_Lemma1.4_ChangeGamma.py
new file mode 100644
index 0000000000000000000000000000000000000000..842880f5263a24f5ff7df10f0eb488291fe6d1d1
--- /dev/null
+++ b/src/Plot_Angle_Lemma1.4_ChangeGamma.py
@@ -0,0 +1,666 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+from HelperFunctions import *
+from ClassifyMin import *
+
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# import tikzplotlib
+# # from pylab import *
+# from tikzplotlib import save as tikz_save
+
+
+# Needed ?
+mpl.use('pdf')
+
+# from subprocess import Popen, PIPE
+#import sys
+
+###################### makePlot.py #########################
+# Generalized Plot-Script giving the option to define
+# quantity of interest and the parameter it depends on
+# to create a plot
+#
+# Input: Define y & x for "x-y plot" as Strings
+# - Run the 'Cell-Problem' for the different Parameter-Points
+# (alternatively run 'Compute_MuGamma' if quantity of interest
+# is q3=muGamma for a significant Speedup)
+
+###########################################################
+
+
+
+# figsize argument takes inputs in inches
+# and we have the width of our document in pts.
+# To set the figure size we construct a function
+# to convert from pts to inches and to determine
+# an aesthetic figure height using the golden ratio:
+# def set_size(width, fraction=1):
+# """Set figure dimensions to avoid scaling in LaTeX.
+#
+# Parameters
+# ----------
+# width: float
+# Document textwidth or columnwidth in pts
+# fraction: float, optional
+# Fraction of the width which you wish the figure to occupy
+#
+# Returns
+# -------
+# fig_dim: tuple
+# Dimensions of figure in inches
+# """
+# # Width of figure (in pts)
+# fig_width_pt = width * fraction
+#
+# # Convert from pt to inches
+# inches_per_pt = 1 / 72.27
+#
+# # Golden ratio to set aesthetic figure height
+# # https://disq.us/p/2940ij3
+# golden_ratio = (5**.5 - 1) / 2
+#
+# # Figure width in inches
+# fig_width_in = fig_width_pt * inches_per_pt
+# # Figure height in inches
+# fig_height_in = fig_width_in * golden_ratio
+#
+# fig_dim = (fig_width_in, fig_height_in)
+#
+# return fig_dim
+#
+
+
+
+def format_func(value, tick_number):
+ # # find number of multiples of pi/2
+ # N = int(np.round(2 * value / np.pi))
+ # if N == 0:
+ # return "0"
+ # elif N == 1:
+ # return r"$\pi/2$"
+ # elif N == 2:
+ # return r"$\pi$"
+ # elif N % 2 > 0:
+ # return r"${0}\pi/2$".format(N)
+ # else:
+ # return r"${0}\pi$".format(N // 2)
+ # find number of multiples of pi/2
+ N = int(np.round(4 * value / np.pi))
+ if N == 0:
+ return "0"
+ elif N == 1:
+ return r"$\pi/4$"
+ elif N == 2:
+ return r"$\pi/2$"
+ elif N % 2 > 0:
+ return r"${0}\pi/2$".format(N)
+ else:
+ return r"${0}\pi$".format(N // 2)
+
+
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+
+# TODO
+# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
+# - Also Add option to plot Minimization Output
+
+
+# ----- Setup Paths -----
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+#---------------------------------------------------------------
+
+print('---- Input parameters: -----')
+mu1 = 1.0 #10.0
+# lambda1 = 10.0
+rho1 = 1.0
+alpha = 5.0
+beta = 10.0
+# alpha = 2.0
+# beta = 2.0
+theta = 1.0/8.0 #1.0/4.0
+
+lambda1 = 0.0
+# gamma = 1.0/4.0
+
+# TEST:
+alpha=3.0;
+
+
+
+
+# # INTERESTING!:
+alpha = 5.0
+beta = 10.0
+theta= 1/8
+
+
+
+alpha_Values
+
+
+
+gamma = 'infinity' #Elliptic Setting
+# gamma = '0' #Hyperbolic Setting
+# gamma = 0.5
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+
+# --- define Interval of x-va1ues:
+xmin = 0.01
+xmax = 0.41
+xmax = 0.99
+
+
+Jumps = False
+
+
+numPoints = 2000
+numPoints = 1000
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+
+
+
+
+for theta in X_Values:
+
+ print('Situation of Lemma1.4')
+ q12 = 0.0
+ q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+ q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+ b1 = prestrain_b1(rho1, beta, alpha,theta)
+ b2 = prestrain_b2(rho1, beta, alpha,theta)
+ b3 = 0.0
+
+ q3 = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
+
+
+
+ G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3, mu1, rho1)
+
+ Y_Values.append(angle)
+
+
+
+print("(Output) Values of " + yName + ": ", Y_Values)
+
+
+idx = find_nearestIdx(Y_Values, 0)
+print(' Idx of value closest to 0', idx)
+ValueClose = Y_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', ValueClose)
+
+
+
+# Find Indices where the difference between the next one is larger than epsilon...
+jump_idx = []
+jump_xValues = []
+jump_yValues = []
+tmp = X_Values[0]
+for idx, x in enumerate(X_Values):
+ print(idx, x)
+ if idx > 0:
+ if abs(Y_Values[idx]-Y_Values[idx-1]) > 1:
+ print('jump candidate')
+ jump_idx.append(idx)
+ jump_xValues.append(x)
+ jump_yValues.append(Y_Values[idx])
+
+
+
+
+
+
+
+print("Jump Indices", jump_idx)
+print("Jump X-values:", jump_xValues)
+print("Jump Y-values:", jump_yValues)
+
+y_plotValues = [Y_Values[0]]
+x_plotValues = [X_Values[0]]
+# y_plotValues.extend(jump_yValues)
+for i in jump_idx:
+ y_plotValues.extend([Y_Values[i-1], Y_Values[i]])
+ x_plotValues.extend([X_Values[i-1], X_Values[i]])
+
+
+y_plotValues.append(Y_Values[-1])
+# x_plotValues = [X_Values[0]]
+# x_plotValues.extend(jump_xValues)
+x_plotValues.append(X_Values[-1])
+
+
+print("y_plotValues:", y_plotValues)
+print("x_plotValues:", x_plotValues)
+# Y_Values[np.diff(y) >= 0.5] = np.nan
+
+
+#get values bigger than jump position
+# gamma = infty
+# x_rest = X_Values[X_Values>x_plotValues[1]]
+# Y_Values = np.array(Y_Values) #convert the np array
+# y_rest = Y_Values[X_Values>x_plotValues[1]]
+#
+#
+# # gamma = 0
+# x_rest = X_Values[X_Values>x_plotValues[3]]
+# Y_Values = np.array(Y_Values) #convert the np array
+# y_rest = Y_Values[X_Values>x_plotValues[3]]
+
+# gamma between
+# Y_Values = np.array(Y_Values) #convert the np array
+# X_Values = np.array(X_Values) #convert the np array
+#
+# x_one = X_Values[X_Values>x_plotValues[3]]
+# # ax.scatter(X_Values, Y_Values)
+# y_rest = Y_Values[X_Values>x_plotValues[3]]
+# ax.plot(X_Values[X_Values>0.135], Y_Values[X_Values<0.135])
+#
+#
+#
+
+
+# y_rest = Y_Values[np.nonzero(X_Values>x_plotValues[1]]
+# print('X_Values:', X_Values)
+# print('Y_Values:', Y_Values)
+# print('x_rest:', x_rest)
+# print('y_rest:', y_rest)
+# print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) )
+
+
+
+
+# --- Convert to numpy array
+Y_Values = np.array(Y_Values)
+X_Values = np.array(X_Values)
+
+# ---------------- Create Plot -------------------
+
+#--- change plot style: SEABORN
+# plt.style.use("seaborn-paper")
+
+
+#--- Adjust gobal matplotlib variables
+# mpl.rcParams['pdf.fonttype'] = 42
+# mpl.rcParams['ps.fonttype'] = 42
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+
+
+# plt.rc('font', family='serif', serif='Times')
+# plt.rc('font', family='serif')
+# # plt.rc('text', usetex=True) #also works...
+# plt.rc('xtick', labelsize=8)
+# plt.rc('ytick', labelsize=8)
+# plt.rc('axes', labelsize=8)
+
+
+
+
+
+#---- Scale Figure apropriately to fit tex-File Width
+# width = 452.9679
+
+# width as measured in inkscape
+width = 6.28 *0.5
+height = width / 1.618
+
+#setup canvas first
+fig = plt.figure() #main
+# fig, ax = plt.subplots()
+# fig, (ax, ax2) = plt.subplots(ncols=2)
+# fig,axes = plt.subplots(nrows=1,ncols=2,figsize=(width,height)) # more than one plot
+
+
+# fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97) #TEST
+
+
+# TEST
+# mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1)
+# fig = plt.figure(figsize=(width+0.1,height+0.1))
+
+
+# mpl.rcParams['figure.figsize'] = (width,height)
+# fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=(width,height)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=set_size(width))
+# fig = plt.subplots(1, 1, figsize=set_size(width))
+
+# --- To create a figure half the width of your document:#
+# fig = plt.figure(figsize=set_size(width, fraction=0.5))
+
+
+
+#--- You must select the correct size of the plot in advance
+# fig.set_size_inches(3.54,3.54)
+
+ax = plt.axes((0.15,0.18,0.8,0.8))
+# ax = plt.axes((0.1,0.1,0.5,0.8))
+# ax = plt.axes((0.1,0.1,1,1))
+# ax = plt.axes()
+
+# ax.spines['right'].set_visible(False)
+# ax.spines['left'].set_visible(False)
+# ax.spines['bottom'].set_visible(False)
+# ax.spines['top'].set_visible(False)
+# ax.tick_params(axis='x',which='major',direction='out',length=10,width=5,color='red',pad=15,labelsize=15,labelcolor='green',
+# labelrotation=15)
+# ax.tick_params(axis='x',which='major', direction='out',pad=5,labelsize=10)
+# ax.tick_params(axis='y',which='major', length=5, width=1, direction='out',pad=5,labelsize=10)
+ax.tick_params(axis='x',which='major', direction='out',pad=3)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+ax.xaxis.set_major_locator(MultipleLocator(0.05))
+ax.xaxis.set_minor_locator(MultipleLocator(0.025))
+
+
+#---- print data-types
+print(ax.xaxis.get_major_locator())
+print(ax.xaxis.get_minor_locator())
+print(ax.xaxis.get_major_formatter())
+print(ax.xaxis.get_minor_formatter())
+
+#---- Hide Ticks or Labels
+# ax.yaxis.set_major_locator(plt.NullLocator())
+# ax.xaxis.set_major_formatter(plt.NullFormatter())
+
+#---- Reducing or Increasing the Number of Ticks
+# ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+# ax.yaxis.set_major_locator(plt.MaxNLocator(3))
+
+
+#----- Fancy Tick Formats
+ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+
+
+
+
+
+
+# --- manually change ticks&labels:
+# ax.set_xticks([0.2,1])
+# ax.set_xticklabels(['pos1','pos2'])
+
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+
+a=ax.yaxis.get_major_locator()
+b=ax.yaxis.get_major_formatter()
+c = ax.get_xticks()
+d = ax.get_xticklabels()
+print('xticks:',c)
+print('xticklabels:',d)
+
+ax.grid(True,which='major',axis='both',alpha=0.3)
+
+
+
+
+
+
+# plt.figure()
+
+# f,ax=plt.subplots(1)
+
+# plt.title(r''+ yName + '-Plot')
+# plt.plot(X_Values, Y_Values,linewidth=2, '.k')
+# plt.plot(X_Values, Y_Values,'.k',markersize=1)
+# plt.plot(X_Values, Y_Values,'.',markersize=0.8)
+
+# plt.plot(X_Values, Y_Values)
+
+# ax.plot([[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])
+
+
+
+# Gamma = '0'
+# ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b')
+#
+# ax.plot([x_plotValues[1],x_plotValues[3]], [y_plotValues[2],y_plotValues[3]] , 'b')
+#
+# ax.plot(x_rest, y_rest, 'b')
+
+
+# Gamma between
+
+# x jump values (gamma 0): [0.13606060606060608, 0.21090909090909093]
+
+# ax.plot([[0,jump_xValues[0]], [0, 0]] , 'b')
+# ax.plot([jump_xValues[0],xmin], [y_plotValues[2],y_plotValues[2]] , 'b')
+
+# ax.plot([[0,0.13606060606060608], [0, 0]] , 'b')
+# ax.plot([[0.13606060606060608,xmin], [(math.pi/2),(math.pi/2)]], 'b')
+
+# jump_xValues[0]
+
+
+
+# --- leave out jumps:
+# ax.scatter(X_Values, Y_Values)
+
+ax.set_xlabel(r"volume fraction $\theta$")
+ax.set_ylabel(r"angle $\alpha$")
+
+
+if Jumps:
+
+ # --- leave out jumps:
+ if gamma == 'infinity':
+ ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]] , 'royalblue')
+ ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue')
+
+
+
+ # ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]])
+ # ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]])
+
+
+
+
+ # ax.plot(X_Values[X_Values>0.136], Y_Values[X_Values>0.136])
+ # ax.plot(X_Values[X_Values<0.135], Y_Values[X_Values<0.135])
+ # ax.scatter(X_Values, Y_Values)
+ # ax.plot(X_Values, Y_Values)
+
+ # plt.plot(x_plotValues, y_plotValues,'.')
+ # plt.scatter(X_Values, Y_Values, alpha=0.3)
+ # plt.scatter(X_Values, Y_Values)
+ # plt.plot(X_Values, Y_Values,'.')
+ # plt.plot([X_Values[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])
+ # plt.axis([0, 6, 0, 20])
+
+ # ax.set_xlabel(r"volume fraction $\theta$", size=11)
+ # ax.set_ylabel(r"angle $\angle$", size=11)
+ # ax.set_xlabel(r"volume fraction $\theta$")
+ # # ax.set_ylabel(r"angle $\angle$")
+ # ax.set_ylabel(r"angle $\alpha$")
+ # plt.ylabel('$\kappa$')
+
+ # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+ # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
+
+
+
+
+ # Plot every other line.. not the jumps...
+
+ if gamma == '0':
+ tmp = 1
+ for idx, x in enumerate(x_plotValues):
+ if idx > 0 and tmp == 1:
+ # plt.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]] )
+ ax.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]], 'royalblue', zorder=2)
+ tmp = 0
+ else:
+ tmp = 1
+
+ # plt.plot([x_plotValues[0],x_plotValues[1]] ,[y_plotValues[0],y_plotValues[1]] )
+ # plt.plot([x_plotValues[2],x_plotValues[3]] ,[y_plotValues[2],y_plotValues[3]] )
+ # plt.plot([x_plotValues[4],x_plotValues[5]] ,[y_plotValues[4],y_plotValues[5]] )
+ # plt.plot([x_plotValues[6],x_plotValues[7]] ,[y_plotValues[6],y_plotValues[7]] )
+
+
+ for x in jump_xValues:
+ plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1, zorder=1)
+ # plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', label=r'$\theta_*$')
+
+ # plt.axvline(x_plotValues[1],ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')
+
+ # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+ # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+ # plt.legend()
+
+
+ # -- SETUP LEGEND
+ # ax.legend(prop={'size': 11})
+ # ax.legend()
+
+ # ------------------ SAVE FIGURE
+ # tikzplotlib.save("TesTout.tex")
+ # plt.close()
+ # mpl.rcParams.update(mpl.rcParamsDefault)
+
+ # plt.savefig("graph.pdf",
+ # #This is simple recomendation for publication plots
+ # dpi=1000,
+ # # Plot will be occupy a maximum of available space
+ # bbox_inches='tight',
+ # )
+ # plt.savefig("graph.pdf")
+
+
+
+ # ---- ADD additional scatter:
+ # ax.scatter(X_Values,Y_Values,s=1,c='black',zorder=4)
+
+ # Find transition point
+ lastIdx = len(Y_Values)-1
+
+ for idx, y in enumerate(Y_Values):
+ if idx != lastIdx:
+ if abs(y-0) < 0.01 and abs(Y_Values[idx+1] - 0) > 0.05:
+ transition_point1 = X_Values[idx+1]
+ print('transition point1:', transition_point1 )
+ if abs(y-0.5*np.pi) < 0.01 and abs(Y_Values[idx+1] -0.5*np.pi)>0.01:
+ transition_point2 = X_Values[idx]
+ print('transition point2:', transition_point2 )
+ if abs(y-0) > 0.01 and abs(Y_Values[idx+1] - 0) < 0.01:
+ transition_point3 = X_Values[idx+1]
+ print('transition point3:', transition_point3 )
+
+ # Add transition Points:
+ if gamma == '0':
+ ax.scatter([transition_point1, transition_point2],[np.pi/2,np.pi/2],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2+0.012 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+ else:
+ ax.scatter([transition_point1, transition_point2, transition_point3 ],[np.pi/2,np.pi/2,0 ],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2 +0.011 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point3 +0.009 , 0+0.08, r"$3$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+else:
+ ax.scatter(X_Values,Y_Values,s=1, marker='o', cmap=None, norm=None, facecolor = 'blue',
+ edgecolor = 'none', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+ ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8 ])
+ labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$',r'$5\pi/8$']
+ ax.set_yticklabels(labels)
+ # ax.set_yticks([1.570786327, np.pi/2 ])
+ # labels = [r'$\pi/2-0.0005 $' , r'$\pi/2$']
+ # ax.set_yticklabels(labels)
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Angle-Theta.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+# figureheight = '\\figureheight',
+# figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
+# #---------------------------------------------------------------
diff --git a/src/Plot_CurvatureLemma1.4_alpha3.0.py b/src/Plot_CurvatureLemma1.4_alpha3.0.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f516be57839d72a4cb35d365628eb5433f4659d
--- /dev/null
+++ b/src/Plot_CurvatureLemma1.4_alpha3.0.py
@@ -0,0 +1,408 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+from HelperFunctions import *
+from ClassifyMin import *
+
+import matplotlib.ticker as ticker
+# from subprocess import Popen, PIPE
+#import sys
+
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+###################### makePlot.py #########################
+# Generalized Plot-Script giving the option to define
+# quantity of interest and the parameter it depends on
+# to create a plot
+#
+# Input: Define y & x for "x-y plot" as Strings
+# - Run the 'Cell-Problem' for the different Parameter-Points
+# (alternatively run 'Compute_MuGamma' if quantity of interest
+# is q3=muGamma for a significant Speedup)
+
+###########################################################
+
+def format_func(value, tick_number):
+ # find number of multiples of pi/2
+ N = int(np.round(2 * value / np.pi))
+ if N == 0:
+ return "0"
+ elif N == 1:
+ return r"$\pi/2$"
+ elif N == 2:
+ return r"$\pi$"
+ elif N % 2 > 0:
+ return r"${0}\pi/2$".format(N)
+ else:
+ return r"${0}\pi$".format(N // 2)
+
+
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+
+# TODO
+# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
+# - Also Add option to plot Minimization Output
+
+
+# ----- Setup Paths -----
+InputFile = "/inputs/cellsolver.parset"
+OutputFile = "/outputs/output.txt"
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+#---------------------------------------------------------------
+
+print('---- Input parameters: -----')
+# mu1 = 10.0
+# # lambda1 = 10.0
+rho1 = 1.0
+# alpha = 5.0
+# beta = 10.0
+# theta = 1.0/4.0
+
+
+# mu1 = 10.0
+mu1 = 1.0
+# lambda1 = 10.0
+# rho1 = 10.0
+# alpha = 5.0
+# beta = 2.0
+beta = 10.0
+theta = 1.0/4.0
+
+theta = 1.0/2.0
+# theta = 1.0/12.0
+
+
+# TesT:
+alpha = 3.0
+
+lambda1 = 0.0
+gamma = 1.0/4.0
+
+gamma = 'infinity'
+gamma = '0'
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+# Optional - TODO? :
+# -Ask User for Input ...
+# function = input("Enter value you want to plot (y-value):\n")
+# print(f'You entered {function}')
+# parameter = input("Enter Parameter this value depends on (x-value) :\n")
+# print(f'You entered {parameter}')
+
+# -Add Option to change NumberOfElements used for computation of Cell-Problem
+
+
+# --- Define Quantity of interest:
+# Options: 'q1', 'q2', 'q3', 'q12' ,'q21', 'q31', 'q13' , 'q23', 'q32' , 'b1', 'b2' ,'b3'
+# TODO: EXTRA (MInimization Output) 'Minimizer (norm?)' 'angle', 'type', 'curvature'
+# yName = 'q12'
+# # yName = 'b1'
+# yName = 'q3'
+yName = 'angle'
+yName = 'curvature'
+
+# --- Define Parameter this function/quantity depends on:
+# Options: mu1 ,lambda1, rho1 , alpha, beta, theta, gamma
+# xName = 'theta'
+# xName = 'gamma'
+# xName = 'lambda1'
+xName = 'theta'
+# xName = 'alpha'
+
+
+# --- define Interval of x-values:
+xmin = 0
+xmax = 30
+
+# xmin = 0.245
+# xmax = 0.99
+#
+#
+# xmin = 0.14
+# xmax = 0.19
+
+# xmin = 0.01
+# xmax = 3.0
+
+xmin = 0.01
+xmax = 0.4
+
+
+
+numPoints = 200
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+
+
+
+for theta in X_Values:
+# for alpha in X_Values:
+
+ print('Situation of Lemma1.4')
+ q12 = 0.0
+ q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+ q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+ b1 = prestrain_b1(rho1, beta, alpha,theta)
+ b2 = prestrain_b2(rho1, beta, alpha,theta)
+ b3 = 0.0
+ if gamma == '0':
+ q3 = q2
+ if gamma == 'infinity':
+ q3 = q1
+
+ if yName == 'q1': # TODO: Better use dictionary?...
+ print('q1 used')
+ Y_Values.append(q1)
+ elif yName =='q2':
+ print('q2 used')
+ Y_Values.append(q2)
+ elif yName =='q3':
+ print('q3 used')
+ Y_Values.append(q3)
+ elif yName =='q12':
+ print('q12 used')
+ Y_Values.append(q12)
+ elif yName =='b1':
+ print('b1 used')
+ Y_Values.append(b1)
+ elif yName =='b2':
+ print('b2 used')
+ Y_Values.append(b2)
+ elif yName =='b3':
+ print('b3 used')
+ Y_Values.append(b3)
+ elif yName == 'angle' or yName =='type' or yName =='curvature':
+ G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3, mu1, rho1)
+ if yName =='angle':
+ print('angle used')
+ Y_Values.append(angle)
+ if yName =='type':
+ print('angle used')
+ Y_Values.append(type)
+ if yName =='curvature':
+ print('curvature used')
+ Y_Values.append(curvature)
+
+
+print("(Output) Values of " + yName + ": ", Y_Values)
+
+
+idx = find_nearestIdx(Y_Values, 0)
+print(' Idx of value closest to 0:', idx)
+ValueClose = Y_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*:', ValueClose)
+print('Theta(Idx) with curvature closest to 0:', ValueClose)
+
+
+
+
+
+# Find Indices where the difference between the next one is larger than epsilon...
+jump_idx = []
+jump_xValues = []
+jump_yValues = []
+tmp = X_Values[0]
+for idx, x in enumerate(X_Values):
+ print(idx, x)
+ if idx > 0:
+ if abs(Y_Values[idx]-Y_Values[idx-1]) > 1:
+ print('jump candidate')
+ jump_idx.append(idx)
+ jump_xValues.append(x)
+ jump_yValues.append(Y_Values[idx])
+
+
+
+
+
+
+
+print("Jump Indices", jump_idx)
+print("Jump X-values:", jump_xValues)
+print("Jump Y-values:", jump_yValues)
+
+y_plotValues = [Y_Values[0]]
+x_plotValues = [X_Values[0]]
+# y_plotValues.extend(jump_yValues)
+for i in jump_idx:
+ y_plotValues.extend([Y_Values[i-1], Y_Values[i]])
+ x_plotValues.extend([X_Values[i-1], X_Values[i]])
+
+
+y_plotValues.append(Y_Values[-1])
+# x_plotValues = [X_Values[0]]
+# x_plotValues.extend(jump_xValues)
+x_plotValues.append(X_Values[-1])
+
+
+print("y_plotValues:", y_plotValues)
+print("x_plotValues:", x_plotValues)
+# Y_Values[np.diff(y) >= 0.5] = np.nan
+
+
+#get values bigger than jump position
+x_rest = X_Values[X_Values>x_plotValues[1]]
+
+Y_Values = np.array(Y_Values) #convert the np array
+
+y_rest = Y_Values[X_Values>x_plotValues[1]]
+# y_rest = Y_Values[np.nonzero(X_Values>x_plotValues[1]]
+print('X_Values:', X_Values)
+print('Y_Values:', Y_Values)
+print('x_rest:', x_rest)
+print('y_rest:', y_rest)
+print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) )
+
+
+# --- Convert to numpy array
+Y_Values = np.array(Y_Values)
+X_Values = np.array(X_Values)
+
+
+# ---------------- Create Plot -------------------
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+# width as measured in inkscape
+width = 6.28 *0.5
+height = width / 1.618
+fig = plt.figure()
+ax = plt.axes((0.15,0.18,0.8,0.8))
+ax.tick_params(axis='x',which='major', direction='out',pad=3)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+ax.xaxis.set_major_locator(MultipleLocator(0.05))
+ax.xaxis.set_minor_locator(MultipleLocator(0.025))
+ax.grid(True,which='major',axis='both',alpha=0.3)
+# plt.figure()
+# f,ax=plt.subplots(1)
+
+
+ax.set_xlabel(r"volume fraction $\theta$")
+ax.set_ylabel(r"curvature $\kappa$")
+# plt.xlabel(xName)
+# plt.ylabel(yName)
+# plt.ylabel('$\kappa$')
+# ax.grid(True)
+
+# Add transition Points
+if gamma == '0':
+ # transition_point1 = 0.13663316582914573 #alpha = 5.0
+ # transition_point2 = 0.20899497487437185
+
+ transition_point1 = 0.2133016508254127 # #alpha = 3.0
+ transition_point2 = 0.29554277138569285
+ # transition_point2 = 0.29574287143571787
+
+ plt.axvline(transition_point1,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1)
+ plt.axvline(transition_point2,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1)
+
+ ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue')
+ ax.plot(X_Values[np.where(np.logical_and(X_Values>jump_xValues[0], X_Values<jump_xValues[1])) ], Y_Values[np.where(np.logical_and(X_Values>jump_xValues[0] ,X_Values<jump_xValues[1] ))] ,'royalblue')
+ ax.plot(X_Values[X_Values>jump_xValues[1]], Y_Values[X_Values>jump_xValues[1]], 'royalblue')
+ # ax.plot(x_plotValues,y_plotValues, 'royalblue')
+ ax.scatter([transition_point1, transition_point2],[jump_yValues[0], jump_yValues[1]],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , jump_yValues[0]-0.02, r"$4$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2+0.012 , jump_yValues[1]+0.02, r"$5$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+if gamma == 'infinity':
+ # transition_point1 = 0.13663316582914573 #alpha = 5.0
+ # transition_point2 = 0.1929145728643216
+ # transition_point3 = 0.24115577889447234
+
+
+ transition_point1 = 0.2133016508254127 #alpha = 3.0
+ transition_point2 = 0.2771335667833917
+ transition_point3 = 0.3313606803401701
+ plt.axvline(transition_point1,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1)
+ plt.axvline(transition_point2,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1)
+ plt.axvline(transition_point3,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1)
+ ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue')
+ ax.plot(X_Values[X_Values>jump_xValues[0]], Y_Values[X_Values>jump_xValues[0]], 'royalblue')
+
+ idx1 = find_nearestIdx(X_Values, transition_point1)
+ idx2 = find_nearestIdx(X_Values, transition_point2)
+ print('idx1', idx1)
+ print('idx2', idx2)
+ Y_TP1 = Y_Values[idx1]
+ Y_TP2 = Y_Values[idx2]
+ print('Y_TP1', Y_TP1)
+ print('Y_TP2', Y_TP2)
+
+
+ ax.scatter([transition_point1, transition_point2],[Y_TP1, Y_TP2],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+ edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+ ax.text(transition_point1-0.02 , Y_TP1-0.02, r"$6$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
+ )
+
+ ax.text(transition_point2+0.015 , Y_TP2+0.020, r"$7$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5))
+# for x in jump_xValues:
+# plt.axvline(x,ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-Theta.pdf')
+
+# plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.legend()
+plt.show()
+# #---------------------------------------------------------------
diff --git a/src/Plot_elasticQuantities.py b/src/Plot_elasticQuantities.py
new file mode 100644
index 0000000000000000000000000000000000000000..c9aa66c31a74b8c2c662399a2d9ae04569a9ec9d
--- /dev/null
+++ b/src/Plot_elasticQuantities.py
@@ -0,0 +1,325 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+from HelperFunctions import *
+from ClassifyMin import *
+
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+
+
+def find_nearest(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return array[idx]
+
+
+def find_nearestIdx(array, value):
+ array = np.asarray(array)
+ idx = (np.abs(array - value)).argmin()
+ return idx
+
+
+
+# TODO
+# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
+# - Also Add option to plot Minimization Output
+
+
+# ----- Setup Paths -----
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+#---------------------------------------------------------------
+
+print('---- Input parameters: -----')
+
+mu1 = 1.0
+# mu1 = 10.0
+# lambda1 = 10.0
+rho1 = 1.0
+alpha = 2.0
+beta = 5.0
+theta = 1.0/4.0
+
+lambda1 = 0.0
+# gamma = 1.0/4.0
+
+gamma = 'infinity' #Elliptic Setting
+# gamma = '0' #Hyperbolic Setting
+# gamma = 0.5
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+# --- define Interval of x-va1ues:
+xmin = 0.0
+xmax = 1.0
+
+
+numPoints = 200
+Theta_Values = np.linspace(xmin, xmax, num=numPoints)
+print('Theta_Values:', Theta_Values)
+
+
+
+
+
+
+B1_Values = []
+B2_Values = []
+
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+
+b1_Vec = np.vectorize(prestrain_b1)
+b2_Vec = np.vectorize(prestrain_b2)
+
+harmonicMeanVec = np.vectorize(harmonicMean)
+arithmeticMeanVec = np.vectorize(arithmeticMean)
+
+Theta_Values = np.array(Theta_Values)
+
+# B1_Values_alphaNeg1 = b1_Vec(rho1, beta, -1.0,Theta_Values)
+# B1_Values_alphaNeg10 = b1_Vec(rho1, beta, -10.0,Theta_Values)
+# B1_Values_alpha2= b1_Vec(rho1, beta, 2.0 ,Theta_Values)
+# B1_Values_alpha10= b1_Vec(rho1, beta, 10.0 ,Theta_Values)
+# # B2_Values = b2_Vec(rho1, beta, alpha,Theta_Values)
+# B2_Values_alphaNeg1 = b2_Vec(rho1, beta, -1.0,Theta_Values)
+# B2_Values_alphaNeg10 = b2_Vec(rho1, beta, -10.0,Theta_Values)
+# B2_Values_alpha2= b2_Vec(rho1, beta, 2.0 ,Theta_Values)
+# B2_Values_alpha10= b2_Vec(rho1, beta, 10.0 ,Theta_Values)
+
+Q1_Values_beta1 = (1.0/6.0)*harmonicMeanVec(mu1, 1.0, Theta_Values)
+Q1_Values_beta2 = (1.0/6.0)*harmonicMeanVec(mu1, 2.0, Theta_Values)
+Q1_Values_beta5 = (1.0/6.0)*harmonicMeanVec(mu1, 5.0, Theta_Values)
+Q1_Values_beta10 = (1.0/6.0)*harmonicMeanVec(mu1, 10.0, Theta_Values)
+
+Q2_Values_beta1 = (1.0/6.0)*arithmeticMeanVec(mu1, 1.0, Theta_Values)
+Q2_Values_beta2 = (1.0/6.0)*arithmeticMeanVec(mu1, 2.0, Theta_Values)
+Q2_Values_beta5 = (1.0/6.0)*arithmeticMeanVec(mu1, 5.0, Theta_Values)
+Q2_Values_beta10 = (1.0/6.0)*arithmeticMeanVec(mu1, 10.0, Theta_Values)
+
+print("Q1_Values_beta1 ", Q1_Values_beta1 )
+
+# --- Convert to numpy array
+# B1_Values = np.array(B1_Values)
+# B2_Values = np.array(B2_Values)
+Q1_Values_beta1 = np.array(Q1_Values_beta1 )
+Q1_Values_beta2 = np.array(Q1_Values_beta2 )
+Q1_Values_beta5 = np.array(Q1_Values_beta5 )
+Q1_Values_beta10 = np.array(Q1_Values_beta10 )
+
+Q2_Values_beta1 = np.array(Q2_Values_beta1 )
+Q2_Values_beta2 = np.array(Q2_Values_beta2 )
+Q2_Values_beta5 = np.array(Q2_Values_beta5 )
+Q2_Values_beta10 = np.array(Q2_Values_beta10 )
+
+# ---------------- Create Plot -------------------
+
+#--- change plot style: SEABORN
+# plt.style.use("seaborn-paper")
+
+
+#--- Adjust gobal matplotlib variables
+# mpl.rcParams['pdf.fonttype'] = 42
+# mpl.rcParams['ps.fonttype'] = 42
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+# mpl.rcParams['axes.grid'] = True
+
+# plt.rc('font', family='serif', serif='Times')
+# plt.rc('font', family='serif')
+# # plt.rc('text', usetex=True) #also works...
+# plt.rc('xtick', labelsize=8)
+# plt.rc('ytick', labelsize=8)
+# plt.rc('axes', labelsize=8)
+
+
+
+
+
+#---- Scale Figure apropriately to fit tex-File Width
+# width = 452.9679
+
+# width as measured in inkscape
+# width = 6.28 *0.5
+width = 6.28
+
+height = width / 1.618
+height = width / 2.5
+#setup canvas first
+fig = plt.figure() #main
+# fig, ax = plt.subplots()
+# fig, (ax, ax2) = plt.subplots(ncols=2)
+fig,ax = plt.subplots(nrows=1,ncols=3,figsize=(width,height)) # more than one plot
+
+# --- set overall Title
+# fig.suptitle('Example of a Single Legend Shared Across Multiple Subplots')
+
+# fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97) #TEST
+
+
+# TEST
+# mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1)
+# fig = plt.figure(figsize=(width+0.1,height+0.1))
+
+
+# mpl.rcParams['figure.figsize'] = (width,height)
+# fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=(width,height)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+# fig = plt.figure(figsize=set_size(width))
+# fig = plt.subplots(1, 1, figsize=set_size(width))
+
+# --- To create a figure half the width of your document:#
+# fig = plt.figure(figsize=set_size(width, fraction=0.5))
+
+
+
+#--- You must select the correct size of the plot in advance
+# fig.set_size_inches(3.54,3.54)
+
+
+# ---- TODO ?:
+# ax[0] = plt.axes((0.15,0.18,0.8,0.8))
+
+
+# ax.tick_params(axis='x',which='major', direction='out',pad=3)
+# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+# ax.xaxis.set_major_locator(MultipleLocator(0.1))
+# ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+# a=ax.yaxis.get_major_locator()
+# b=ax.yaxis.get_major_formatter()
+# c = ax.get_xticks()
+# d = ax.get_xticklabels()
+# print('xticks:',c)
+# print('xticklabels:',d)
+
+ax[0].grid(True,which='major',axis='both',alpha=0.3)
+ax[1].grid(True,which='major',axis='both',alpha=0.3)
+ax[2].grid(True,which='major',axis='both',alpha=0.3)
+# ax.plot(Theta_Values,B1_Values , 'royalblue')
+# ax.plot(Theta_Values,B2_Values , 'royalblue')
+
+l1 = ax[0].plot(Theta_Values,Q1_Values_beta1 , label=r"$\theta_\mu = 1.0$")
+l2 = ax[0].plot(Theta_Values,Q1_Values_beta2 , label=r"$\theta_\mu = 2.0$")
+l3 = ax[0].plot(Theta_Values,Q1_Values_beta5 , label=r"$\theta_\mu = 5.0$")
+l4 = ax[0].plot(Theta_Values,Q1_Values_beta10 , label=r"$\theta_\mu = 10.0$")
+
+ax[0].set_xlabel(r"volume fraction $\theta$")
+ax[0].set_ylabel(r" $q_1$")
+ax[0].xaxis.set_major_locator(MultipleLocator(0.25))
+# Labels to use in the legend for each line
+line_labels = [r"$\theta_\mu = 1.0$", r"$\theta_\mu = 2.0$", r"$\theta_\mu = 5.0$", r"$\theta_\mu = 10.0$"]
+
+
+ax[1].plot(Theta_Values,Q2_Values_beta1 , label=r"$\theta_\rho = 1.0$")
+ax[1].plot(Theta_Values,Q2_Values_beta2 , label=r"$\theta_\rho = 2.0$")
+ax[1].plot(Theta_Values,Q2_Values_beta5 , label=r"$\theta_\rho = 5.0$")
+ax[1].plot(Theta_Values,Q2_Values_beta10 , label=r"$\theta_\rho = 10.0$")
+
+ax[1].set_xlabel(r"volume fraction $\theta$")
+ax[1].set_ylabel(r" $q_2$")
+ax[1].xaxis.set_major_locator(MultipleLocator(0.25))
+# ax[1].xaxis.set_minor_locator(MultipleLocator(0.05))
+
+
+ax[2].plot(Theta_Values,Q1_Values_beta1/Q2_Values_beta1 , label=r"$\theta_\rho = 1.0$")
+ax[2].plot(Theta_Values,Q1_Values_beta2/Q2_Values_beta2 , label=r"$\theta_\rho = 2.0$")
+ax[2].plot(Theta_Values,Q1_Values_beta5/Q2_Values_beta5 , label=r"$\theta_\rho = 5.0$")
+ax[2].plot(Theta_Values,Q1_Values_beta10/Q2_Values_beta10 , label=r"$\theta_\rho = 10.0$")
+
+ax[2].set_xlabel(r"volume fraction $\theta$")
+ax[2].set_ylabel(r" $q_1/q_2$")
+ax[2].xaxis.set_major_locator(MultipleLocator(0.25))
+
+
+
+
+plt.subplots_adjust(wspace=0.4, hspace=0)
+plt.tight_layout()
+
+
+
+
+
+
+
+fig.legend([l1, l2, l3, l4], # The line objects
+ labels=line_labels, # The labels for each line
+ loc="center right", # Position of legend
+ borderaxespad=0.15 # Small spacing around legend box
+ # title="Legend Title" # Title for the legend
+ )
+
+# Adjust the scaling factor to fit your legend text completely outside the plot
+# (smaller value results in more space being made for the legend)
+plt.subplots_adjust(right=0.8)
+
+# ------------------ SAVE FIGURE
+# tikzplotlib.save("TesTout.tex")
+# plt.close()
+# mpl.rcParams.update(mpl.rcParamsDefault)
+
+# plt.savefig("graph.pdf",
+# #This is simple recomendation for publication plots
+# dpi=1000,
+# # Plot will be occupy a maximum of available space
+# bbox_inches='tight',
+# )
+# plt.savefig("graph.pdf")
+
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-q1q2-Theta.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+# figureheight = '\\figureheight',
+# figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
+# #---------------------------------------------------------------
diff --git a/src/plot_ElasticRatio.py b/src/plot_ElasticRatio.py
new file mode 100644
index 0000000000000000000000000000000000000000..1da41871db5a2da9f64bb8f628d0c9abc7d92ae6
--- /dev/null
+++ b/src/plot_ElasticRatio.py
@@ -0,0 +1,252 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+
+import time
+
+# ----------- SETUP PATHS
+# InputFile = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+InputFile = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+# -------------------------- Input Parameters --------------------
+# mu1 = 10.0 # TODO : here must be the same values as in the Parset for computeMuGamma
+mu1 = 1.0
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+beta = 5.0
+theta = 1.0/4.0
+
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+# gamma = 1.0
+
+# gamma = 5.0
+
+#added
+# lambda1 = 10.0
+lambda1 = 0.0
+
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+
+print('lambda1: ', lambda1)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+#
+# gamma_min = 0.5
+# gamma_max = 1.0
+#
+# # gamma_min = 1
+# # gamma_max = 1
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=3)
+# # #
+# # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
+# print('(Input) Gamma_Values:', Gamma_Values)
+
+print('type of gamma:', type(gamma))
+# # #
+# Gamma_Values = ['0', 'infinity']
+Gamma_Values = ['infinity']
+Gamma_Values = ['0']
+print('(Input) Gamma_Values:', Gamma_Values)
+
+for gamma in Gamma_Values:
+
+ print('Run for gamma = ', gamma)
+ print('type of gamma:', type(gamma))
+ # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+ # # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+ # print('Test MuGamma:', muGamma)
+
+ # ------- Options --------
+ # print_Cases = True
+ # print_Output = True
+
+ #TODO
+ # generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
+ generalCase = False
+
+ # make_3D_plot = True
+ # make_3D_PhaseDiagram = True
+ make_2D_plot = False
+ make_2D_PhaseDiagram = False
+ make_3D_plot = False
+ make_3D_PhaseDiagram = False
+ # make_2D_plot = True
+ make_2D_PhaseDiagram = True
+ #
+
+ # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
+
+ # SamplePoints_3D = 10 # Number of sample points in each direction
+ # SamplePoints_2D = 10 # Number of sample points in each direction
+ SamplePoints_3D = 300 # Number of sample points in each direction
+ # SamplePoints_3D = 150 # Number of sample points in each direction
+ SamplePoints_3D = 100 # Number of sample points in each direction
+ # SamplePoints_3D = 200 # Number of sample points in each direction
+ # SamplePoints_3D = 400 # Number of sample points in each direction
+ # SamplePoints_2D = 7500 # Number of sample points in each direction
+ SamplePoints_2D = 4000 # 4000 # Number of sample points in each direction
+ SamplePoints_2D = 10 # 4000 # Number of sample points in each direction
+
+ # if make_3D_PhaseDiagram:
+ # alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+ # # alphas_ = np.linspace(-10, 10, SamplePoints_3D)
+ # # betas_ = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
+ # # betas_ = np.linspace(0.01,20.01,SamplePoints_3D) # FULL Range
+ # # betas_ = np.linspace(0.01,0.99,SamplePoints_3D) # weird part
+ # betas_ = np.linspace(1.01,40.01,SamplePoints_3D) #TEST !!!!! For Beta <1 weird tings happen...
+ # thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+ #
+ #
+ # alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+ # classifyMin_anaVec = np.vectorize(classifyMin_ana)
+ #
+ # # Get MuGamma values ...
+ # GetMuGammaVec = np.vectorize(GetMuGamma)
+ # muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
+ # # Classify Minimizers....
+ # G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+ #
+ # # G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1, True, True)
+ # # print('size of G:', G.shape)
+ # # print('G:', G)
+ #
+ # # Option to print angles
+ # # print('angles:', angles)
+ #
+ #
+ # # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # # T = Out[2]
+ # # --- Write to VTK
+ #
+ # GammaString = str(gamma)
+ # VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
+ # gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+ # print('Written to VTK-File:', VTKOutputName )
+
+ if make_2D_PhaseDiagram:
+ # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+ # alphas_ = np.linspace(0, 1, SamplePoints_2D)
+ thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+ alphas_ = np.linspace(-5, 5, SamplePoints_2D)
+ # alphas_ = np.linspace(-5, 15, SamplePoints_2D)
+ # thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
+ betas_ = 10.0
+ # alphas_ = -0.5
+ # betas_ = np.linspace(1.01,10.01,SamplePoints_3D) #TEST !!!!! For Beta <1 weird tings happen...
+
+ alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+
+
+ harmonicMeanVec = np.vectorize(harmonicMean)
+ arithmeticMeanVec = np.vectorize(arithmeticMean)
+ #
+ # q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
+ # q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
+
+ GetMuGammaVec = np.vectorize(GetMuGamma)
+ muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+
+ q1 = harmonicMeanVec(mu1, betas, thetas)
+ q2 = arithmeticMeanVec(mu1, betas, thetas)
+
+ # G, angles, Tq1 = harmonicMeanVec(mu1, betas, thetas)ypes, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
+ # print('size of G:', G.shape)
+ # print('G:', G)
+ # print('Types:', Types)
+ # Out = classifyMin_anaVec(alphas,betas,thetas)
+ # T = Out[2]
+ # --- Write to VTK
+ # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+
+ elasticRatio = q1/q2
+
+ print('type( q1) :', type(q1))
+
+ print('q1:', q1)
+ print('q2:', q2)
+ print('q1/q2:', q1/q2)
+
+ GammaString = str(gamma)
+ VTKOutputName = "outputs/ElasticRatio" #+ "Gamma_" + GammaString
+ gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'elasticRatio': elasticRatio} )
+ print('Written to VTK-File:', VTKOutputName )
+
+
+ # --- Make 3D Scatter plot
+ if(make_3D_plot or make_2D_plot):
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+ colors = cm.plasma(Types)
+ # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+ # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
+ if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
+ if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
+ # cbar=plt.colorbar(pnt3d)
+ # cbar.set_label("Values (units)")
+ plt.axvline(x = 8, color = 'b', linestyle = ':', label='$q_1$')
+ plt.axhline(y = 0.083333333, color = 'b', linestyle = ':', label='$q_1$')
+
+ ax.set_xlabel('alpha')
+ ax.set_ylabel('beta')
+ if make_3D_plot: ax.set_zlabel('theta')
+ plt.show()
+
+
+
+
+
+# ALTERNATIVE
+# colors = ("red", "green", "blue")
+# groups = ("Type 1", "Type2", "Type3")
+#
+# # Create plot
+# fig = plt.figure()
+# ax = fig.add_subplot(1, 1, 1)
+#
+# for data, color, group in zip(Types, colors, groups):
+# # x, y = data
+# ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
+#
+# plt.title('Matplot scatter plot')
+# plt.legend(loc=2)
+# plt.show()