diff --git a/src/Plot_Angle_Lemma1.4.py b/src/Plot_Angle_Lemma1.4.py
index de8e97c2dcebc582786defe2860ab24a5e2045f6..9691882b1b8cc13f7736116b50afac5fe41ed6a8 100644
--- a/src/Plot_Angle_Lemma1.4.py
+++ b/src/Plot_Angle_Lemma1.4.py
@@ -165,7 +165,7 @@ lambda1 = 0.0
 gamma = 1.0/4.0
 
 gamma = 'infinity'  #Elliptic Setting
-gamma = '0'       #Hyperbolic Setting
+# gamma = '0'       #Hyperbolic Setting
 # gamma = 0.5
 
 
diff --git a/src/Plot_CurvatureLemma1.4.py b/src/Plot_CurvatureLemma1.4.py
index f2c9f2de302d95434296b8b634da09653cd2796e..1d0d493daec58d19a44dde3d14c5d50807ae6090 100644
--- a/src/Plot_CurvatureLemma1.4.py
+++ b/src/Plot_CurvatureLemma1.4.py
@@ -14,6 +14,11 @@ 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
@@ -165,7 +170,7 @@ xmax = 0.4
 
 
 
-numPoints = 100
+numPoints = 200
 X_Values = np.linspace(xmin, xmax, num=numPoints)
 print(X_Values)
 
@@ -297,77 +302,90 @@ print('y_rest:', y_rest)
 print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) )
 
 
-# ---------------- Create Plot -------------------
-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]])
-# ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b')
-# ax.plot(x_rest, y_rest, 'b')
-
-
-
-# ax.plot(X_Values, Y_Values)
-# ax.scatter(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])
+# --- Convert to numpy array
+Y_Values = np.array(Y_Values)
+X_Values = np.array(X_Values)
 
-ax.plot(X_Values[X_Values>jump_xValues[0]], Y_Values[X_Values>jump_xValues[0]] ,'b')
-ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'b')
 
-plt.xlabel(xName)
+# ---------------- 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.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
-# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
-
-
-
-
-ax.grid(True)
-
-# # if angle PLOT :
-# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
-# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
-#
-# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
-#
-# # Plot every other line.. not the jumps...
-# 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]] ,'b')
-#         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 = 'g',alpha=0.5, linestyle = 'dashed')
-
-# plt.axvline(x_plotValues[1],ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')
+# plt.ylabel('$\kappa$')
+# ax.grid(True)
+
+# Add transition Points
+if gamma == '0':
+    transition_point1 =  0.13663316582914573
+    transition_point2 =  0.20899497487437185
+    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
+    transition_point2 = 0.1929145728643216
+    transition_point3 = 0.24115577889447234
+    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$')
diff --git a/src/Plot_Prestrain_Lemma1.4.py b/src/Plot_Prestrain_Lemma1.4.py
index b8e051d2572ef7e1cec341458fb3410bdd2c25c3..adef5d8e9a2800f34f8cddba3601947c9f9b74f0 100644
--- a/src/Plot_Prestrain_Lemma1.4.py
+++ b/src/Plot_Prestrain_Lemma1.4.py
@@ -9,40 +9,11 @@ import re
 import matlab.engine
 from HelperFunctions import *
 from ClassifyMin import *
-# from mpl_toolkits.mplot3d import axes3d
-from mpl_toolkits import mplot3d
-
-import matplotlib.ticker as ticker
-# 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)
-
-###########################################################
-
-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)
-
 
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
 
 
 
@@ -65,8 +36,12 @@ def find_nearestIdx(array, value):
 
 
 # ----- Setup Paths -----
-InputFile  = "/inputs/cellsolver.parset"
-OutputFile = "/outputs/output.txt"
+# 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
@@ -84,28 +59,21 @@ 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 = 10.0
 # lambda1 = 10.0
 rho1 = 1.0
-alpha = 5.0
-beta = 2.0
+alpha = 2.0
+beta = 5.0
 theta = 1.0/4.0
-# theta = 1.0/12.0
-
 
 lambda1 = 0.0
-gamma = 1.0/4.0
+# gamma = 1.0/4.0
 
-gamma = 'infinity'
-gamma = '0'
+gamma = 'infinity'  #Elliptic Setting
+# gamma = '0'       #Hyperbolic Setting
+# gamma = 0.5
 
 
 print('mu1: ', mu1)
@@ -117,384 +85,286 @@ print('gamma:', gamma)
 print('----------------------------')
 
 
+# --- define Interval of x-va1ues:
+xmin = 0.0
+xmax = 1.0
 
-# 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
+numPoints = 20
+Theta_Values = np.linspace(xmin, xmax, num=numPoints)
+print('Theta_Values:', Theta_Values)
 
 
-# --- 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'
+B1_Values = []
+B2_Values = []
 
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
 
-xName = 'theta'
-xName = 'alpha'
+b1_Vec = np.vectorize(prestrain_b1)
+b2_Vec = np.vectorize(prestrain_b2)
 
+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)
 
-# --- define Interval of x-values:
-xmin = 2
-xmax = 6
+# print('B1_Values:', B1_Values)
+# print('B2_Values:', B2_Values)
 
 
-numPoints = 400
-X_Values = np.linspace(xmin, xmax, num=numPoints)
-print(X_Values)
 
+# --- Convert to numpy array
+# B1_Values = np.array(B1_Values)
+# B2_Values  = np.array(B2_Values)
+B1_Values_alphaNeg1  = np.array(B1_Values_alphaNeg1)
+B1_Values_alphaNeg10  = np.array(B1_Values_alphaNeg10)
+B1_Values_alpha2  = np.array(B1_Values_alpha2)
+B1_Values_alpha10  = np.array(B1_Values_alpha10)
+B2_Values_alphaNeg1  = np.array(B2_Values_alphaNeg1)
+B2_Values_alphaNeg10  = np.array(B2_Values_alphaNeg10)
+B2_Values_alpha2  = np.array(B2_Values_alpha2)
+B2_Values_alpha10  = np.array(B2_Values_alpha10)
+# ---------------- Create Plot -------------------
 
-Y_Values = []
+#--- change plot style:  SEABORN
+# plt.style.use("seaborn-paper")
 
-SamplePoints_2D = 20
 
-prestrain_b1Vec = np.vectorize(prestrain_b1)
-prestrain_b2Vec = np.vectorize(prestrain_b2)
+#--- 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)
 
-alphas_ = np.linspace(-20, 20, SamplePoints_2D)
-thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
-#fix to one value:
-betas_ = 2.0;
-alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
 
-B1 = prestrain_b1Vec(rho1, betas, alphas, thetas)
-B2 = prestrain_b2Vec(rho1, betas, alphas, thetas)
-print('B1:', B1)
 
-B1 = np.array(B1)  #convert the np array
-B2 = np.array(B2)
 
 
+#---- Scale Figure apropriately to fit tex-File Width
+# width = 452.9679
 
-f,ax=plt.subplots(1)
+# width as measured in inkscape
+# width = 6.28 *0.5
+width = 6.28
 
-ax.plot(thetas, B1)
-ax.scatter(thetas, B1)
+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=2,figsize=(width,height)) # more than one plot
 
+# --- set overall Title
+# fig.suptitle('Example of a Single Legend Shared Across Multiple Subplots')
 
-plt.xlabel("$q_3$")
-plt.xlabel("$\gamma$")
-plt.ylabel("angle")
-ax.grid(True)
+# fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97)  #TEST
 
 
-# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
-# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+# TEST
+# mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1)
+# fig = plt.figure(figsize=(width+0.1,height+0.1))
 
-# 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))
+# 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))
 
-# 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.legend()
-# ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+#--- You must select the correct size of the plot in advance
+# fig.set_size_inches(3.54,3.54)
 
-# 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$']
-ax.set_yticks([0, np.pi/8, np.pi/4 ])
-labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
-ax.set_yticklabels(labels)
 
+# ---- TODO ?:
+# ax[0] = plt.axes((0.15,0.18,0.8,0.8))
 
-# Plot Gamma Value that is closest to q3_star
-ax.axvline(x = gammaClose, color = 'b', linestyle = 'dashed', label='$\gamma^*$')
-ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.5)
-ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.5)
 
+# 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.plot(Theta_Values,B1_Values , 'royalblue')
+# ax.plot(Theta_Values,B2_Values , 'royalblue')
 
+l1 = ax[0].plot(Theta_Values,B1_Values_alphaNeg1 , label=r"$\theta_\rho = -1.0$")
+l2 = ax[0].plot(Theta_Values,B1_Values_alphaNeg10 , label=r"$\theta_\rho = -10.0$")
+l3 = ax[0].plot(Theta_Values,B1_Values_alpha2 , label=r"$\theta_\rho = 2.0$")
+l4 = ax[0].plot(Theta_Values,B1_Values_alpha10 , label=r"$\theta_\rho = 10.0$")
+ax[0].set_xlabel(r"volume fraction $\theta$")
+ax[0].set_ylabel(r"prestrain $b_1$")
 
-# plt.xlabel("$q_3(\gamma)$")
-plt.xlabel("$\gamma$")
-plt.ylabel("angle")
-plt.legend(loc='upper center')
+# Labels to use in the legend for each line
+line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -10.0$", r"$\theta_\rho = 2.0$", r"$\theta_\rho = 10.0$"]
 
-plt.show()
 
+ax[1].plot(Theta_Values,B2_Values_alphaNeg1 , label=r"$\theta_\rho = -1.0$")
+ax[1].plot(Theta_Values,B2_Values_alphaNeg10 , label=r"$\theta_\rho = -10.0$")
+ax[1].plot(Theta_Values,B2_Values_alpha2 , label=r"$\theta_\rho = 2.0$")
+ax[1].plot(Theta_Values,B2_Values_alpha10 , label=r"$\theta_\rho = 10.0$")
+ax[1].set_xlabel(r"volume fraction $\theta$")
+ax[1].set_ylabel(r"prestrain $b_2$")
 
+plt.subplots_adjust(wspace=0.4, hspace=0)
+plt.tight_layout()
+# ax.plot(Theta_Values,B2_Values_alphaNeg1 )
+# ax.plot(Theta_Values,B2_Values_alphaNeg10 )
+# ax.plot(Theta_Values,B2_Values_alpha2 )
+# ax.plot(Theta_Values,B2_Values_alpha10 )
 
 
+# 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')
 
-#
-# fig = plt.figure(figsize=(14,6))
-#
-# # `ax` is a 3D-aware axis instance because of the projection='3d' keyword argument to add_subplot
-# ax = fig.add_subplot(1, 2, 1, projection='3d')
-#
-# p = ax.plot_surface(alphas, thetas, B1, rstride=4, cstride=4, linewidth=0)
-#
-# # surface_plot with color grading and color bar
-# ax = fig.add_subplot(1, 2, 2, projection='3d')
-# p = ax.plot_surface(alphas, thetas, B1, rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
-# cb = fig.colorbar(p, shrink=0.5)
-#
-#
-#
-# fig = plt.figure()
-# ax = plt.axes(projection='3d')
-#
-# # fig = plt.figure(figsize=(6,6))
-# # ax = fig.add_subplot(111, projection='3d')
-# ax = fig.add_subplot(111)
-#
-# # f,ax=plt.subplots(1)
-#
-# # Plot a 3D surface
-# # ax.plot_surface(alphas, thetas, B1,cmap='viridis', edgecolor='none')
-# # ax.plot_surface(alphas, thetas, B1)
-# # ax.plot_surface(alphas, thetas, B1)
-#
-# ax.plot3D(alphas, thetas, B1, 'gray')
-# # ax = fig.add_subplot(211, projection='3d')
-# #
-# #
-# # # Plot a 3D surface
-# # ax.plot_surface(alphas, thetas, B1)
-#
-# plt.show()
-
-# 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()
+# jump_xValues[0]
 
 
-#
-#
-# # 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('angle 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)
-#
-#
-#
-# # 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]) )
-#
-#
-# # ---------------- Create Plot -------------------
-# 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]])
-# # ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b')
-# # ax.plot(x_rest, y_rest, 'b')
-#
-#
-#
+
+# --- leave out jumps:
+# ax.scatter(X_Values, Y_Values)
+
+
+
+
+
+
+
+# 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])
-#
-# plt.xlabel(xName)
-# plt.ylabel(yName)
-#
-# # plt.ylabel('$\kappa$')
-#
-# # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
-# # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
-#
-# 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_formatter(plt.FuncFormatter(format_func))
-#
-# # Plot every other line.. not the jumps...
-# 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]] )
-# #         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 = 'g',alpha=0.5, linestyle = 'dashed')
-#
-# # 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()
-# plt.show()
+
+# 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)
+
+
+# plt.ylabel('$\kappa$')
+
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
+
+
+
+
+
+# -- SETUP LEGEND
+# ax.legend(prop={'size': 11})
+# ax[0].legend()
+# ax[1].legend(
+# loc='upper left',
+# bbox_to_anchor=(1,1))
+
+# Create the legend
+# handles, labels = ax.get_legend_handles_labels()
+
+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-Prestrain-Theta.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+#            figureheight = '\\figureheight',
+#            figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
 # #---------------------------------------------------------------