From af804d3bdf1aaadfc455e5e5897aa31ef839184c Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 15 Dec 2021 10:30:24 +0100
Subject: [PATCH] Update plot routines

---
 src/PhaseDiagram.py                     |  23 +-
 src/PhaseDiagram_intermediateGamma.py   | 354 +++++++++++
 src/Plot_Angle_Alpha.py                 | 785 ++++++++++++++++++++++++
 src/Plot_Angle_Lemma1.4V3.py            | 770 +++++++++++++++++++++++
 src/Plot_Angle_Lemma1.4V3_Transition.py | 770 +++++++++++++++++++++++
 src/Plot_Curvature_Alpha.py             | 732 ++++++++++++++++++++++
 src/Plot_Curvature_Lemma1.4V3.py        | 689 +++++++++++++++++++++
 7 files changed, 4114 insertions(+), 9 deletions(-)
 create mode 100644 src/PhaseDiagram_intermediateGamma.py
 create mode 100644 src/Plot_Angle_Alpha.py
 create mode 100644 src/Plot_Angle_Lemma1.4V3.py
 create mode 100644 src/Plot_Angle_Lemma1.4V3_Transition.py
 create mode 100644 src/Plot_Curvature_Alpha.py
 create mode 100644 src/Plot_Curvature_Lemma1.4V3.py

diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py
index 2f269e7c..c7f32e5c 100644
--- a/src/PhaseDiagram.py
+++ b/src/PhaseDiagram.py
@@ -115,11 +115,11 @@ mu1 = 1.0
 rho1 = 1.0
 alpha = 2.0
 beta = 2.0
-beta = 10.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 = 'infinity'
 # gamma = 0.5
 # gamma = 0.25
 # gamma = 1.0
@@ -159,8 +159,9 @@ print('----------------------------')
 
 print('type of gamma:', type(gamma))
 # # #
-Gamma_Values = ['0', 'infinity']
-# Gamma_Values = ['infinity']
+# Gamma_Values = ['0', 'infinity']
+Gamma_Values = ['infinity']
+Gamma_Values = ['0']
 print('(Input) Gamma_Values:', Gamma_Values)
 
 for gamma in Gamma_Values:
@@ -180,13 +181,13 @@ for gamma in Gamma_Values:
     generalCase = False
 
     # make_3D_plot = True
-    make_3D_PhaseDiagram = True
+    # make_3D_PhaseDiagram = True
     make_2D_plot = False
     make_2D_PhaseDiagram = False
     make_3D_plot = False
-    # make_3D_PhaseDiagram = False
+    make_3D_PhaseDiagram = False
     # make_2D_plot = True
-    # make_2D_PhaseDiagram = True
+    make_2D_PhaseDiagram = True
     #
 
     # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
@@ -232,7 +233,7 @@ for gamma in Gamma_Values:
     # 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 = 200 # 4000 # Number of sample points in each direction
 
     if make_3D_PhaseDiagram:
         alphas_ = np.linspace(-20, 20, SamplePoints_3D)
@@ -307,7 +308,11 @@ for gamma in Gamma_Values:
         # betas_  = np.linspace(0.01,40.01,1)
         #fix to one value:
         # betas_ = 2.0;
-        betas_ = 10.0;
+        # betas_ = 10.0;
+        betas_ = 5.0;
+
+
+
 
         # TEST
         # alphas_ = np.linspace(-8, 8, SamplePoints_2D)
diff --git a/src/PhaseDiagram_intermediateGamma.py b/src/PhaseDiagram_intermediateGamma.py
new file mode 100644
index 00000000..a142890b
--- /dev/null
+++ b/src/PhaseDiagram_intermediateGamma.py
@@ -0,0 +1,354 @@
+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
+# print(sys.executable)
+
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters):   alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+#   compute q1,q2,b1,b2 from Formula
+#       Option 1.1 :
+#           set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+#       Option 1.2 :
+#           compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+#   compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+#
+#
+# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+#                 OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+#     # ------------------------------------ get mu_gamma ------------------------------
+#     # ---Scenario 1.1: extreme regimes
+#     if gamma == '0':
+#         print('extreme regime: gamma = 0')
+#         mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+#         print("mu_gamma:", mu_gamma)
+#     elif gamma == 'infinity':
+#         print('extreme regime: gamma = infinity')
+#         mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta)   # = q1
+#         print("mu_gamma:", mu_gamma)
+#     else:
+#         # --- Scenario 1.2:  compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+#         # print("Run computeMuGamma for Gamma = ", gamma)
+#         with open(InputFilePath, 'r') as file:
+#             filedata = file.read()
+#         filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+#         # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+#         filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+#         filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+#         filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+#         filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+#         f = open(InputFilePath,'w')
+#         f.write(filedata)
+#         f.close()
+#         # --- Run Cell-Problem
+#
+#         # Check Time
+#         # t = time.time()
+#         # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+#         #                                      capture_output=True, text=True)
+#         # --- Run Cell-Problem_muGama   -> faster
+#         # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+#         #                                              capture_output=True, text=True)
+#         # --- Run Compute_muGamma (2D Problem much much faster)
+#
+#         subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+#                                                              capture_output=True, text=True)
+#         # print('elapsed time:', time.time() - t)
+#
+#         #Extract mu_gamma from Output-File                                           TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+#         with open(OutputFilePath, 'r') as file:
+#             output = file.read()
+#         tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+#         s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+#         mu_gamma = float(s[0])
+#         # print("mu_gamma:", mu_gammaValue)
+#     # --------------------------------------------------------------------------------------
+#     return mu_gamma
+#
+
+
+
+# ----------- 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
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+# beta = 10.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...
+
+Gamma_Values = [0.05, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5]
+# Gamma_Values = ['infinity']
+print('(Input) Gamma_Values:', Gamma_Values)
+# #
+for gamma in Gamma_Values:
+
+
+        # 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
+
+
+    # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+    # q1 = harmonicMean(mu1, beta, theta)
+    # q2 = arithmeticMean(mu1, beta, theta)
+    # --- Set q12
+    # q12 = 0.0  # (analytical example)              # TEST / TODO read from Cell-Output
+
+
+
+
+
+    # b1 = prestrain_b1(rho1, beta, alpha, theta)
+    # b2 = prestrain_b2(rho1, beta, alpha, theta)
+    #
+    # print('---- Input parameters: -----')
+    # print('mu1: ', mu1)
+    # print('rho1: ', rho1)
+    # print('alpha: ', alpha)
+    # print('beta: ', beta)
+    # print('theta: ', theta)
+    # print("q1: ", q1)
+    # print("q2: ", q2)
+    # print("mu_gamma: ", mu_gamma)
+    # print("q12: ", q12)
+    # print("b1: ", b1)
+    # print("b2: ", b2)
+    # print('----------------------------')
+    # print("machine epsilon", sys.float_info.epsilon)
+
+    # G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
+    # Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+    # print("Test", Test)
+
+    # ---------------------- 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_2D = 200 # Number of sample points in each direction
+
+
+
+    if make_3D_PhaseDiagram:
+        alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+        betas_  = np.linspace(0.01,40.01,SamplePoints_3D)
+        thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+        # print('type of alphas', type(alphas_))
+        # print('Test:', type(np.array([mu_gamma])) )
+        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!!!
+        # print('size of G:', G.shape)
+        # print('G:', G)
+        # 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)
+        # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+        # alphas_ = np.linspace(-5, 15, SamplePoints_2D)
+        # thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
+
+
+        # good range:
+        # alphas_ = np.linspace(9, 10, SamplePoints_2D)
+        # thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
+
+        #bEFORE PrestrainChange:
+        # alphas_ = np.linspace(8, 10, SamplePoints_2D)
+        # thetas_ = np.linspace(0.05,0.16,SamplePoints_2D)
+
+        alphas_ = np.linspace(-2, 1, SamplePoints_2D)
+        thetas_ = np.linspace(0.4,0.6,SamplePoints_2D)
+
+            # alphas_ = np.linspace(8, 12, SamplePoints_2D)
+            # thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
+        # betas_  = np.linspace(0.01,40.01,1)
+        #fix to one value:
+        # betas_ = 2.0;
+        betas_ = 10.0;
+        alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+
+
+        if generalCase:
+            classifyMin_matVec = np.vectorize(classifyMin_mat)
+            GetCellOutputVec = np.vectorize(GetCellOutput, otypes=[np.ndarray, np.ndarray])
+            Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,lambda1, InputFilePath ,OutputFilePath )
+
+
+            # print('type of Q:', type(Q))
+            # print('Q:', Q)
+            G, angles, Types, curvature = classifyMin_matVec(Q,B)
+
+        else:
+            classifyMin_anaVec = np.vectorize(classifyMin_ana)
+            GetMuGammaVec = np.vectorize(GetMuGamma)
+            muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+            G, angles, Types, 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"
+
+        GammaString = str(gamma)
+        VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
+        gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+        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()
+        # plt.savefig('common_labels.png', dpi=300)
+        # print('T:', T)
+        # print('Type 1 occured here:', np.where(T == 1))
+        # print('Type 2 occured here:', np.where(T == 2))
+
+
+        # print(alphas_)
+        # print(betas_)
+
+
+
+
+
+# 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()
diff --git a/src/Plot_Angle_Alpha.py b/src/Plot_Angle_Alpha.py
new file mode 100644
index 00000000..d23e7b7e
--- /dev/null
+++ b/src/Plot_Angle_Alpha.py
@@ -0,0 +1,785 @@
+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 = 3
+beta = 10.0
+theta= 1/8
+
+
+
+
+
+
+gamma = 'infinity'  #Elliptic Setting
+# gamma = '0'       #Hyperbolic Settings
+# 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
+
+
+xmin = -2.0
+xmax = 1.0
+
+
+Jumps = True
+
+
+numPoints = 10000
+# numPoints = 100
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+Angle_alpha0 = []
+Angle_alphaNeg0125 = []
+Angle_alphaNeg025 = []
+Angle_alphaNeg05 = []
+Angle_alphaNeg075 = []
+Angle_alphaNeg1 = []
+Angle_alpha3 = []
+
+Angle_alphaNeg0625 = []
+Angle_alphaNeg0875 = []
+
+theta = 0.5
+
+
+
+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
+    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)
+    # G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg1.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg05 .append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg025.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
+    # Angle_alpha3.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(-0.75,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg075.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
+    # Angle_alpha0.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg0125.append(angle)
+    #
+    # G, angle, Type, curvature = classifyMin_ana(-0.625,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg0625.append(angle)
+    # G, angle, Type, curvature = classifyMin_ana(-0.875,beta,theta, q3,  mu1, rho1)
+    # Angle_alphaNeg0875.append(angle)
+
+
+
+print("(Output) Values of angle: ", 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)
+
+
+Angle_alphaNeg1 = np.array(Angle_alphaNeg1)
+Angle_alphaNeg05 = np.array(Angle_alphaNeg05)
+Angle_alphaNeg025 = np.array(Angle_alphaNeg025)
+Angle_alphaNeg075 = np.array(Angle_alphaNeg075)
+Angle_alpha3 = np.array(Angle_alpha3)
+Angle_alphaNeg0 = np.array(Angle_alpha0)
+Angle_alphaNeg0125 = np.array(Angle_alphaNeg0125)
+
+Angle_alphaNeg0625 = np.array(Angle_alphaNeg0625)
+Angle_alphaNeg0875 = np.array(Angle_alphaNeg0875)
+# ---------------- 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
+width = 6.28
+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.15,0.18,0.6,0.6))
+# 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))
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+ax.xaxis.set_major_locator(MultipleLocator(0.5))
+ax.xaxis.set_minor_locator(MultipleLocator(0.25))
+
+#---- 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.11 , 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_point2+0.06 , np.pi/2-0.02, r"$1$", 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.11 , 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_point2 + 0.06 , 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_point3 +0.04 , 0+0.04, 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)
+         # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1, label=r"$\theta_\rho = -1.0$")
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        #
+        # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
+        # 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.legend(handles=[l1,l2,l3,l4, l5, l6, l7],
+        #           labels= [ r"$\theta_\rho = 0$", r"$\theta_\rho = -0.125$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.75$",  r"$\theta_\rho = -1.0$",  r"$\theta_\rho = 3.0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+       # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0875,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg075,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg0625,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg05,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg025,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alphaNeg0125,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l8 = ax.scatter(X_Values,Angle_alpha0,s=2, marker='s', edgecolor = 'black', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+
+        # l1 = ax.plot(X_Values,Angle_alphaNeg1, color='red', linewidth=1.5, zorder=3, label = 'test')
+        # l2 = ax.plot(X_Values,Angle_alphaNeg0875,color='orangered', linewidth=1.5 ,linestyle = '--' ,zorder=3)
+        # l3 = ax.plot(X_Values,Angle_alphaNeg075,color='orange', linewidth=1.5,linestyle = '--' ,zorder=3)
+        # l4 = ax.plot(X_Values,Angle_alphaNeg0625, linewidth=1.5,linestyle = '--' ,  zorder=3)
+        # l5 = ax.plot(X_Values,Angle_alphaNeg05, color='teal',linestyle = '--', linewidth=1.5 ,  zorder=3)
+        # l6 = ax.plot(X_Values,Angle_alphaNeg025,color='lightskyblue', linewidth=1.5,linestyle = '--',  zorder=3)
+        # l7 = ax.plot(X_Values,Angle_alphaNeg0125,color='dodgerblue', linewidth=1.5,linestyle = '--', zorder=3)
+        # l8 = ax.plot(X_Values,Angle_alpha0, color='blue', linewidth=1.5 ,zorder=3)
+        #
+        #
+        #
+        #
+        # ax.legend(handles=[l1[0],l2[0],l3[0],l4[0], l5[0], l6[0], l7[0], l8[0]],
+        #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+
+
+
+        l1 = ax.plot(X_Values,Y_Values, color='red', linewidth=1.5, zorder=3, label = 'test')
+        # l1 = ax.scatter(X_Values,Y_Values,s=1, marker='o', edgecolor = 'red', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+
+
+# ax.plot(X_Values, Y_Values,   marker='o',  markerfacecolor='orange', markeredgecolor='black', markeredgewidth=1,  linewidth=1, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+
+        # line_labels = [r"$\theta_\rho = -1.0$",r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -  \frac{3}{4}$", r"$\theta_\rho = -  \frac{5}{8}$",r"$\theta_\rho = - 0.5 $" ,r"$\theta_\rho = -  0.25", r"$\theta_\rho = -  \frac{1}{8}" , r"$\theta_\rho = 0$"]
+        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.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
+        #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+
+
+
+
+
+
+
+        # fig.legend([l1, l2, l3, l4],     # The line objects
+        #            labels=line_labels,   # The labels for each line
+        #            # loc="upper center",   # Position of legend
+        #            loc='upperleft', bbox_to_anchor=(1,1),
+        #            borderaxespad=0.15    # Small spacing around legend box
+        #            # title="Legend Title"  # Title for the legend
+        #            )
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Angle-Alpha.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.4V3.py b/src/Plot_Angle_Lemma1.4V3.py
new file mode 100644
index 00000000..bdbcef34
--- /dev/null
+++ b/src/Plot_Angle_Lemma1.4V3.py
@@ -0,0 +1,770 @@
+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 = 3
+beta = 10.0
+theta= 1/8
+
+
+
+
+
+
+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 = 100
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+Angle_alpha0 = []
+Angle_alphaNeg0125 = []
+Angle_alphaNeg025 = []
+Angle_alphaNeg05 = []
+Angle_alphaNeg075 = []
+Angle_alphaNeg1 = []
+Angle_alpha3 = []
+
+Angle_alphaNeg0625 = []
+Angle_alphaNeg0875 = []
+
+
+
+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)
+    G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg1.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg05 .append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg025.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
+    Angle_alpha3.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.75,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg075.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
+    Angle_alpha0.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0125.append(angle)
+
+    G, angle, Type, curvature = classifyMin_ana(-0.625,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0625.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.875,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0875.append(angle)
+
+
+
+print("(Output) Values of angle: ", 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)
+
+
+Angle_alphaNeg1 = np.array(Angle_alphaNeg1)
+Angle_alphaNeg05 = np.array(Angle_alphaNeg05)
+Angle_alphaNeg025 = np.array(Angle_alphaNeg025)
+Angle_alphaNeg075 = np.array(Angle_alphaNeg075)
+Angle_alpha3 = np.array(Angle_alpha3)
+Angle_alphaNeg0 = np.array(Angle_alpha0)
+Angle_alphaNeg0125 = np.array(Angle_alphaNeg0125)
+
+Angle_alphaNeg0625 = np.array(Angle_alphaNeg0625)
+Angle_alphaNeg0875 = np.array(Angle_alphaNeg0875)
+# ---------------- 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
+width = 6.28
+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.15,0.18,0.6,0.6))
+# 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))
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+
+#---- 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)
+         # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1, label=r"$\theta_\rho = -1.0$")
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        #
+        # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
+        # 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.legend(handles=[l1,l2,l3,l4, l5, l6, l7],
+        #           labels= [ r"$\theta_\rho = 0$", r"$\theta_\rho = -0.125$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.75$",  r"$\theta_\rho = -1.0$",  r"$\theta_\rho = 3.0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+       # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0875,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg075,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg0625,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg05,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg025,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alphaNeg0125,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l8 = ax.scatter(X_Values,Angle_alpha0,s=2, marker='s', edgecolor = 'black', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+
+        l1 = ax.plot(X_Values,Angle_alphaNeg1, color='red', linewidth=1.5, zorder=3, label = 'test')
+        l2 = ax.plot(X_Values,Angle_alphaNeg0875,color='orangered', linewidth=1.5 ,linestyle = '--' ,zorder=3)
+        l3 = ax.plot(X_Values,Angle_alphaNeg075,color='orange', linewidth=1.5,linestyle = '--' ,zorder=3)
+        l4 = ax.plot(X_Values,Angle_alphaNeg0625, linewidth=1.5,linestyle = '--' ,  zorder=3)
+        l5 = ax.plot(X_Values,Angle_alphaNeg05, color='teal',linestyle = '--', linewidth=1.5 ,  zorder=3)
+        l6 = ax.plot(X_Values,Angle_alphaNeg025,color='lightskyblue', linewidth=1.5,linestyle = '--',  zorder=3)
+        l7 = ax.plot(X_Values,Angle_alphaNeg0125,color='dodgerblue', linewidth=1.5,linestyle = '--', zorder=3)
+        l8 = ax.plot(X_Values,Angle_alpha0, color='blue', linewidth=1.5 ,zorder=3)
+
+        ax.legend(handles=[l1[0],l2[0],l3[0],l4[0], l5[0], l6[0], l7[0], l8[0]],
+                  labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+                  loc='upper left',
+                  bbox_to_anchor=(1,1))
+
+
+
+# ax.plot(X_Values, Y_Values,   marker='o',  markerfacecolor='orange', markeredgecolor='black', markeredgewidth=1,  linewidth=1, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+
+        # line_labels = [r"$\theta_\rho = -1.0$",r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -  \frac{3}{4}$", r"$\theta_\rho = -  \frac{5}{8}$",r"$\theta_\rho = - 0.5 $" ,r"$\theta_\rho = -  0.25", r"$\theta_\rho = -  \frac{1}{8}" , r"$\theta_\rho = 0$"]
+        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.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
+        #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+
+
+
+
+
+
+
+        # fig.legend([l1, l2, l3, l4],     # The line objects
+        #            labels=line_labels,   # The labels for each line
+        #            # loc="upper center",   # Position of legend
+        #            loc='upperleft', bbox_to_anchor=(1,1),
+        #            borderaxespad=0.15    # Small spacing around legend box
+        #            # title="Legend Title"  # Title for the legend
+        #            )
+
+
+
+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.4V3_Transition.py b/src/Plot_Angle_Lemma1.4V3_Transition.py
new file mode 100644
index 00000000..bdbcef34
--- /dev/null
+++ b/src/Plot_Angle_Lemma1.4V3_Transition.py
@@ -0,0 +1,770 @@
+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 = 3
+beta = 10.0
+theta= 1/8
+
+
+
+
+
+
+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 = 100
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+Angle_alpha0 = []
+Angle_alphaNeg0125 = []
+Angle_alphaNeg025 = []
+Angle_alphaNeg05 = []
+Angle_alphaNeg075 = []
+Angle_alphaNeg1 = []
+Angle_alpha3 = []
+
+Angle_alphaNeg0625 = []
+Angle_alphaNeg0875 = []
+
+
+
+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)
+    G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg1.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg05 .append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg025.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
+    Angle_alpha3.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.75,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg075.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
+    Angle_alpha0.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0125.append(angle)
+
+    G, angle, Type, curvature = classifyMin_ana(-0.625,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0625.append(angle)
+    G, angle, Type, curvature = classifyMin_ana(-0.875,beta,theta, q3,  mu1, rho1)
+    Angle_alphaNeg0875.append(angle)
+
+
+
+print("(Output) Values of angle: ", 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)
+
+
+Angle_alphaNeg1 = np.array(Angle_alphaNeg1)
+Angle_alphaNeg05 = np.array(Angle_alphaNeg05)
+Angle_alphaNeg025 = np.array(Angle_alphaNeg025)
+Angle_alphaNeg075 = np.array(Angle_alphaNeg075)
+Angle_alpha3 = np.array(Angle_alpha3)
+Angle_alphaNeg0 = np.array(Angle_alpha0)
+Angle_alphaNeg0125 = np.array(Angle_alphaNeg0125)
+
+Angle_alphaNeg0625 = np.array(Angle_alphaNeg0625)
+Angle_alphaNeg0875 = np.array(Angle_alphaNeg0875)
+# ---------------- 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
+width = 6.28
+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.15,0.18,0.6,0.6))
+# 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))
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+
+#---- 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)
+         # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1, label=r"$\theta_\rho = -1.0$")
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        #
+        # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
+        # 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.legend(handles=[l1,l2,l3,l4, l5, l6, l7],
+        #           labels= [ r"$\theta_\rho = 0$", r"$\theta_\rho = -0.125$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.75$",  r"$\theta_\rho = -1.0$",  r"$\theta_\rho = 3.0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+       # ---------------------------------------------------------------
+        # l1 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1)
+        # l2 = ax.scatter(X_Values,Angle_alphaNeg0875,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l3 = ax.scatter(X_Values,Angle_alphaNeg075,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l4 = ax.scatter(X_Values,Angle_alphaNeg0625,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg05,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l6 = ax.scatter(X_Values,Angle_alphaNeg025,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alphaNeg0125,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        # l8 = ax.scatter(X_Values,Angle_alpha0,s=2, marker='s', edgecolor = 'black', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+
+        l1 = ax.plot(X_Values,Angle_alphaNeg1, color='red', linewidth=1.5, zorder=3, label = 'test')
+        l2 = ax.plot(X_Values,Angle_alphaNeg0875,color='orangered', linewidth=1.5 ,linestyle = '--' ,zorder=3)
+        l3 = ax.plot(X_Values,Angle_alphaNeg075,color='orange', linewidth=1.5,linestyle = '--' ,zorder=3)
+        l4 = ax.plot(X_Values,Angle_alphaNeg0625, linewidth=1.5,linestyle = '--' ,  zorder=3)
+        l5 = ax.plot(X_Values,Angle_alphaNeg05, color='teal',linestyle = '--', linewidth=1.5 ,  zorder=3)
+        l6 = ax.plot(X_Values,Angle_alphaNeg025,color='lightskyblue', linewidth=1.5,linestyle = '--',  zorder=3)
+        l7 = ax.plot(X_Values,Angle_alphaNeg0125,color='dodgerblue', linewidth=1.5,linestyle = '--', zorder=3)
+        l8 = ax.plot(X_Values,Angle_alpha0, color='blue', linewidth=1.5 ,zorder=3)
+
+        ax.legend(handles=[l1[0],l2[0],l3[0],l4[0], l5[0], l6[0], l7[0], l8[0]],
+                  labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+                  loc='upper left',
+                  bbox_to_anchor=(1,1))
+
+
+
+# ax.plot(X_Values, Y_Values,   marker='o',  markerfacecolor='orange', markeredgecolor='black', markeredgewidth=1,  linewidth=1, zorder=3)
+        # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
+        # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+        # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+
+
+        # line_labels = [r"$\theta_\rho = -1.0$",r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -  \frac{3}{4}$", r"$\theta_\rho = -  \frac{5}{8}$",r"$\theta_\rho = - 0.5 $" ,r"$\theta_\rho = -  0.25", r"$\theta_\rho = -  \frac{1}{8}" , r"$\theta_\rho = 0$"]
+        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.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
+        #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
+        #           loc='upper left',
+        #           bbox_to_anchor=(1,1))
+
+
+
+
+
+
+
+        # fig.legend([l1, l2, l3, l4],     # The line objects
+        #            labels=line_labels,   # The labels for each line
+        #            # loc="upper center",   # Position of legend
+        #            loc='upperleft', bbox_to_anchor=(1,1),
+        #            borderaxespad=0.15    # Small spacing around legend box
+        #            # title="Legend Title"  # Title for the legend
+        #            )
+
+
+
+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_Curvature_Alpha.py b/src/Plot_Curvature_Alpha.py
new file mode 100644
index 00000000..dde04df8
--- /dev/null
+++ b/src/Plot_Curvature_Alpha.py
@@ -0,0 +1,732 @@
+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 = 3
+beta = 10.0
+theta= 1/8
+
+
+
+theta = 0.5
+
+
+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
+
+
+xmin = -2.0
+xmax = 1.0
+
+Jumps = False
+
+
+numPoints = 2000
+# numPoints = 100
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+Curvature_alpha0 = []
+Curvature_alphaNeg0125 = []
+Curvature_alphaNeg025 = []
+Curvature_alphaNeg05 = []
+Curvature_alphaNeg075 = []
+Curvature_alphaNeg1 = []
+Curvature_alpha3 = []
+Curvature_alphaNeg5 = []
+
+
+
+
+
+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
+    q3 = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
+
+    G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3,  mu1, rho1)
+    Y_Values.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg1.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg05 .append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg025.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
+    # Curvature_alpha3.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-0.75,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg075.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
+    # Curvature_alpha0.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg0125.append(curvature)
+    # G, angle, Type, curvature = classifyMin_ana(-5.0,beta,theta, q3,  mu1, rho1)
+    # Curvature_alphaNeg5.append(curvature)
+
+print("(Output) Values of Curvature: ", 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]) > 0.5:
+            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)
+
+
+Curvature_alphaNeg1 = np.array(Curvature_alphaNeg1)
+Curvature_alphaNeg05 = np.array(Curvature_alphaNeg05)
+Curvature_alphaNeg025 = np.array(Curvature_alphaNeg025)
+Curvature_alphaNeg075 = np.array(Curvature_alphaNeg075)
+Curvature_alpha3 = np.array(Curvature_alpha3)
+Curvature_alphaNeg0 = np.array(Curvature_alpha0)
+Curvature_alphaNeg0125 = np.array(Curvature_alphaNeg0125)
+Curvature_alphaNeg5 = np.array(Curvature_alphaNeg5)
+# ---------------- 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
+width = 6.28
+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.15,0.18,0.6,0.6))
+# 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))
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+
+ax.xaxis.set_major_locator(MultipleLocator(0.5))
+ax.xaxis.set_minor_locator(MultipleLocator(0.25))
+
+#---- 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"Curvature $\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')
+
+    # 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':
+        # 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)
+
+
+        plt.axvline(jump_xValues[0],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')
+
+        # plt.axvline(jump_xValues[0],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)
+        #                )
+
+    else :
+        plt.axvline(jump_xValues[0],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))
+        ax.scatter(jump_xValues,jump_yValues,s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
+                                 edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        ax.text(jump_xValues[0]+0.05 , jump_yValues[0]+0.02, r"$6$", 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)
+    # ---------------------------------------------------------------
+    # l1 = ax.scatter(X_Values,Curvature_alphaNeg5,s=1, marker='o',  cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=3)
+    # l2 = ax.scatter(X_Values,Curvature_alphaNeg1,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3, label=r"$\theta_\rho = -1.0$")
+    # l3 = ax.scatter(X_Values,Curvature_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+    # l4 = ax.scatter(X_Values,Curvature_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+    # l5 = ax.scatter(X_Values,Curvature_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+    # l6 = ax.scatter(X_Values,Curvature_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+    # l7 = ax.scatter(X_Values,Curvature_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
+    # l8 = ax.scatter(X_Values,Curvature_alpha3,s=1, marker='o',  cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=4)
+    # # l4 = ax.scatter(X_Values,Curvature_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+    #
+    # ax.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
+    #           labels= [r"$\theta_\rho = -5.0$",  r"$\theta_\rho = -1.0$",r"$\theta_\rho = -0.75$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.125$",  r"$\theta_\rho = 0$",    r"$\theta_\rho = 3.0$"  ],
+    #           loc='upper left',
+    #           bbox_to_anchor=(1,1))
+    # ---------------------------------------------------------------
+
+
+    # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
+    # 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.legend([l1, l2, l3, l4],     # The line objects
+    #            labels=line_labels,   # The labels for each line
+    #            # loc="upper center",   # Position of legend
+    #            loc='upperleft', bbox_to_anchor=(1,1),
+    #            borderaxespad=0.15    # Small spacing around legend box
+    #            # title="Legend Title"  # Title for the legend
+    #            )
+
+
+
+
+
+    # l1 = ax.plot(X_Values,Y_Values, color='red', linewidth=1.5, zorder=3, label = 'test')
+    l1 = ax.scatter(X_Values,Y_Values,s=1, marker='o', edgecolor = 'red', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-Alpha.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+#            figureheight = '\\figureheight',
+#            figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
+# #---------------------------------------------------------------
diff --git a/src/Plot_Curvature_Lemma1.4V3.py b/src/Plot_Curvature_Lemma1.4V3.py
new file mode 100644
index 00000000..171664ac
--- /dev/null
+++ b/src/Plot_Curvature_Lemma1.4V3.py
@@ -0,0 +1,689 @@
+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 = 3
+beta = 10.0
+theta= 1/8
+
+
+
+
+
+
+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 = 100
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+
+Curvature_alpha0 = []
+Curvature_alphaNeg0125 = []
+Curvature_alphaNeg025 = []
+Curvature_alphaNeg05 = []
+Curvature_alphaNeg075 = []
+Curvature_alphaNeg1 = []
+Curvature_alpha3 = []
+Curvature_alphaNeg5 = []
+
+
+
+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(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg1.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg05 .append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg025.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
+    Curvature_alpha3.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-0.75,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg075.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
+    Curvature_alpha0.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg0125.append(curvature)
+    G, angle, Type, curvature = classifyMin_ana(-5.0,beta,theta, q3,  mu1, rho1)
+    Curvature_alphaNeg5.append(curvature)
+
+print("(Output) Values of Curvature: ", 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)
+
+
+Curvature_alphaNeg1 = np.array(Curvature_alphaNeg1)
+Curvature_alphaNeg05 = np.array(Curvature_alphaNeg05)
+Curvature_alphaNeg025 = np.array(Curvature_alphaNeg025)
+Curvature_alphaNeg075 = np.array(Curvature_alphaNeg075)
+Curvature_alpha3 = np.array(Curvature_alpha3)
+Curvature_alphaNeg0 = np.array(Curvature_alpha0)
+Curvature_alphaNeg0125 = np.array(Curvature_alphaNeg0125)
+Curvature_alphaNeg5 = np.array(Curvature_alphaNeg5)
+# ---------------- 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
+width = 6.28
+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.15,0.18,0.6,0.6))
+# 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))
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+
+#---- 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"Curvature $\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')
+
+
+
+
+
+    # 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)
+        l1 = ax.scatter(X_Values,Curvature_alphaNeg5,s=1, marker='o',  cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=3)
+        l2 = ax.scatter(X_Values,Curvature_alphaNeg1,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3, label=r"$\theta_\rho = -1.0$")
+        l3 = ax.scatter(X_Values,Curvature_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        l4 = ax.scatter(X_Values,Curvature_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        l5 = ax.scatter(X_Values,Curvature_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        l6 = ax.scatter(X_Values,Curvature_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
+        l7 = ax.scatter(X_Values,Curvature_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
+        l8 = ax.scatter(X_Values,Curvature_alpha3,s=1, marker='o',  cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=4)
+        # l4 = ax.scatter(X_Values,Curvature_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
+
+
+
+        # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
+        # 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)
+        ax.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
+                  labels= [r"$\theta_\rho = -5.0$",  r"$\theta_\rho = -1.0$",r"$\theta_\rho = -0.75$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.125$",  r"$\theta_\rho = 0$",    r"$\theta_\rho = 3.0$"  ],
+                  loc='upper left',
+                  bbox_to_anchor=(1,1))
+
+        # fig.legend([l1, l2, l3, l4],     # The line objects
+        #            labels=line_labels,   # The labels for each line
+        #            # loc="upper center",   # Position of legend
+        #            loc='upperleft', bbox_to_anchor=(1,1),
+        #            borderaxespad=0.15    # Small spacing around legend box
+        #            # title="Legend Title"  # Title for the legend
+        #            )
+
+
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-Theta.pdf')
+
+
+
+
+# tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')
+
+# tikz_save('fig.tikz',
+#            figureheight = '\\figureheight',
+#            figurewidth = '\\figurewidth')
+
+# ----------------------------------------
+
+
+plt.show()
+# #---------------------------------------------------------------
-- 
GitLab