Add Curvature-Gamma Plots

src/PhaseDiagram.py

@@ -172,14 +172,14 @@ print('----------------------------')
 generalCase = False
 
 # make_3D_plot = True
-# make_3D_PhaseDiagram = True
+make_3D_PhaseDiagram = True
 make_2D_plot = False
-# make_2D_PhaseDiagram = 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 ---
 # q1 = harmonicMean(mu1, beta, theta)
@@ -218,14 +218,28 @@ make_2D_PhaseDiagram = True
 # SamplePoints_3D = 10 # Number of sample points in each direction
 # SamplePoints_2D = 10 # Number of sample points in each direction
 SamplePoints_3D = 300 # Number of sample points in each direction
+# SamplePoints_3D = 100 # Number of sample points in each direction
+# SamplePoints_3D = 200 # Number of sample points in each direction
+# SamplePoints_3D = 400 # Number of sample points in each direction
 SamplePoints_2D = 300 # 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)
+    # betas_  = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
+
+
+    # betas_  = np.linspace(0.01,1.01,SamplePoints_3D)  # weird part
+    betas_  = np.linspace(1,40.01,SamplePoints_3D)     #TEST !!!!!  For Beta <1 weird tings happen...
     thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+
+
+
+    # TEST :
+    # alphas_ = np.linspace(-40, 40, SamplePoints_3D)
+    # betas_  = np.linspace(0.01,80.01,SamplePoints_3D) # Full Range
+
     # print('type of alphas', type(alphas_))
     # print('Test:', type(np.array([mu_gamma])) )
     alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
@@ -267,6 +281,7 @@ if make_2D_PhaseDiagram:
     # betas_  = np.linspace(0.01,40.01,1)
     #fix to one value:
     betas_ = 2.0;
+    # betas_ = 0.5;  #TEST!!!
     alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
 
 

src/PhaseDiagram_PlotScript.py

@@ -6,6 +6,8 @@ from paraview.simple import *
 #### disable automatic camera reset on 'Show'
 paraview.simple._DisableFirstRenderCameraReset()
 
+#--- Run in Terminal with: 'pvbatch PhaseDiagram_PlotScript.py'
+
 
 # ----- CREATE A Phase Diagram for the following cases:
 # 1. (Hyperbolic Case) Gamma = '0'
@@ -23,7 +25,9 @@ curvature = 0
 # phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGammainfinity_100SP.vts'])
 # phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGammainfinity300sp.vts'])
 # phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGammainfinityBetaSmallerOne.vts'])
-phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGamma0300sp.vts'])
+# phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGamma0300sp.vts'])
+# phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGammainfinity_greaterOne.vts'])
+phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGamma0_greaterOne.vts'])
 # phaseDiagram3DGammainfinityvts = XMLStructuredGridReader(FileName=['/home/klaus/Desktop/DUNE/dune-microstructure/outputs/PhaseDiagram3DGamma0.vts'])
 
 
@@ -420,7 +424,8 @@ renderView1.AxesGrid.XAxisUseCustomLabels = 1
 renderView1.AxesGrid.ZAxisUseCustomLabels = 1
 renderView1.AxesGrid.YAxisUseCustomLabels = 1
 renderView1.AxesGrid.XAxisLabels = [-20, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0]
-renderView1.AxesGrid.YAxisLabels = [0 , 10.0, 20.0, 30.0, 40.0]
+# renderView1.AxesGrid.YAxisLabels = [0 , 10.0, 20.0, 30.0, 40.0]
+renderView1.AxesGrid.YAxisLabels = [1 ,5, 10.0,15, 20.0,25, 30.0,35, 40.0]
 # renderView1.AxesGrid.YAxisLabels = [0 ,5.0, 10.0, 15.0, 20.0, 25.0, 30.0,35.0, 40.0]
 renderView1.AxesGrid.ZAxisLabels = [ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
 

src/Plot-Angle-Gamma.py

@@ -176,7 +176,7 @@ print('----------------------------')
 
 gamma_min = 0.01
 gamma_max = 1.0
-Gamma_Values = np.linspace(gamma_min, gamma_max, num=200)    # TODO variable Input Parameters...alpha,beta...
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=20)    # TODO variable Input Parameters...alpha,beta...
 print('(Input) Gamma_Values:', Gamma_Values)
 # mu_gamma = []
 

src/Plot-Angle-q3.py

@@ -352,11 +352,15 @@ ax.set_yticklabels(labels)
 # ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2)
 # ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2)
 
+
+ax.axvline(x = q3_star, color = 'midnightblue', linestyle = 'dashed', label='$q_3^*$')
 ax.axvline(x = q1, color = 'forestgreen', linestyle = 'dashed', label='$q_1$')
 ax.axvline(x = q2, color = 'darkorange', linestyle = 'dashed', label='$q_2$')
 
 
-ax.legend(loc='upper center')
+# ax.legend(loc='upper center')
+ax.legend(loc='upper left', bbox_to_anchor=(0.55, 0.85))
+# ax.legend(loc='upper left')
 
 # plt.xlabel("$q_3(\gamma)$")
 # plt.xlabel("$\gamma$")

src/Plot-Curvature-Gamma.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+import matplotlib.ticker as ticker
+
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# from matplotlib import rc
+# rc('text', usetex=True) # Use LaTeX font
+#
+# import seaborn as sns
+# sns.set(color_codes=True)
+
+
+def format_func(value, tick_number):
+    # # find number of multiples of pi/2
+    # N = int(np.round(2 * value / np.pi))
+    # if N == 0:
+    #     return "0"
+    # elif N == 1:
+    #     return r"$\pi/2$"
+    # elif N == 2:
+    #     return r"$\pi$"
+    # elif N % 2 > 0:
+    #     return r"${0}\pi/2$".format(N)
+    # else:
+    #     return r"${0}\pi$".format(N // 2)
+    # find number of multiples of pi/2
+    N = int(np.round(4 * value / np.pi))
+    if N == 0:
+        return "0"
+    elif N == 1:
+        return r"$\pi/4$"
+    elif N == 2:
+        return r"$\pi/2$"
+    elif N % 2 > 0:
+        return r"${0}\pi/2$".format(N)
+    else:
+        return r"${0}\pi$".format(N // 2)
+
+
+
+def find_nearest(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return array[idx]
+
+
+def find_nearestIdx(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return idx
+
+
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+print('---- Input parameters: -----')
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 2.56140350877193 #2.56140350877193, 4.0852130325814535
+# beta = 2.0  #5.0
+# theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+# theta = 0.075  # 0.5
+# mu1 = 10.0
+# rho1 = 1.0
+alpha = 10.0
+# beta = 40.0
+# theta = 0.125
+mu1 = 10.0
+rho1 = 1.0
+# alpha = 2.0
+beta = 2.0  #5.0
+# theta = 1.0/4.0
+theta = 1.0/8.0
+#
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10.0
+# beta = 2.0  #5.0
+# theta = 1.0/8.0
+# #
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 10.0
+# beta = 2.0
+# theta = 0.125
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 9.0
+# beta = 2.0
+# theta = 0.075
+
+#
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 4.0
+# beta = 10.0
+# theta = 1.0/4.0
+
+# alpha = 10   #1.05263158
+# beta = 40.0  #5.0
+# # theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+#
+#
+# alpha = 2.0
+# beta = 2.0 #5.0
+# theta = 1/4.0
+# rho1 = 1.0
+# mu1 = 10.0
+
+# InterestingParameterSet :
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+
+
+# theta = 0.124242
+# gamma = 0.75
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 1.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=200)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+# Gamma_Values = '0'
+
+
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+print('angles:', angles)
+
+print('curvature:', curvature)
+
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+determinantVec = np.vectorize(determinant)
+
+detValues = determinantVec(q1,q2,muGammas,q12)
+print('detValues:', detValues)
+
+
+detZeroidx = find_nearestIdx(detValues, 0)
+print('idx where det nearest to zero', idx)
+gammaClose = Gamma_Values[detZeroidx]
+print('gammaClose:', gammaClose)
+
+
+# --- Convert to numpy array
+Gamma_Values = np.array(Gamma_Values)
+angles = np.array(angles)
+
+
+curvature = np.array(curvature)
+
+# ---------------- Create Plot -------------------
+# plt.figure()
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+width = 6.28 *0.5
+height = width / 1.618
+fig = plt.figure()
+# ax = plt.axes((0.15,0.21 ,0.75,0.75))
+ax = plt.axes((0.15,0.18 ,0.8,0.75))
+# ax = plt.axes((0.21,0.21 ,0.8,0.75))
+ax.tick_params(axis='x',which='major', direction='out',pad=1)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=1) # changed pad = distance to title to 1 here!
+ax.xaxis.set_major_locator(MultipleLocator(0.1))
+ax.xaxis.set_minor_locator(MultipleLocator(0.05))
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax.grid(True,which='major',axis='both',alpha=0.3)
+
+
+# # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# # plt.title(r'angle$-\gamma$-Plot')
+# plt.plot(Gamma_Values, angles)
+# plt.scatter(Gamma_Values, angles)
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+
+# f,ax=plt.subplots(1)
+
+
+# ax.plot(muGammas, angles)
+# ax.scatter(muGammas, angles
+#
+# ax.plot(Gamma_Values, angles, 'royalblue', zorder=3, )
+ax.plot(Gamma_Values, curvature, 'royalblue', zorder=3, )
+
+# ax.scatter(Gamma_Values, angles)
+
+# ax.set_xlabel(r"$q_3(\gamma)$")
+ax.set_xlabel(r"$\gamma$")
+# ax.set_ylabel(r"curvature $\kappa$")
+ax.set_title(r"curvature $\kappa$", fontsize=9, pad = 4)
+
+# plt.xlabel("$q_3$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# ax.grid(True)
+
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25))
+
+
+# ax.yaxis.set_major_formatter(ticker.FuncFormatter(
+# lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi))
+
+# ---------------------------- show pi values ------------------------------------
+# ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# ax.legend()
+# # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# # ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+# # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+# ---------------------------------------------------------------
+
+
+# ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+# labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+
+
+# OLD :
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+
+
+# Plot Gamma Value that is closest to q3_star
+ax.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', label='$\gamma^*$')
+ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2)
+ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2)
+
+
+ax.legend(loc='upper right')
+
+# plt.xlabel("$q_3(\gamma)$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# plt.legend(loc='upper center')
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-Gamma.pdf')
+
+plt.show()
+
+
+
+
+# plt.figure()
+# plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.xlabel("$\mu_\gamma$")
+# plt.ylabel("angle")
+# plt.legend()
+# plt.show()
+#

src/Plot-Curvature-q3.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+import matplotlib.ticker as ticker
+
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# from matplotlib import rc
+# rc('text', usetex=True) # Use LaTeX font
+#
+# import seaborn as sns
+# sns.set(color_codes=True)
+
+
+def format_func(value, tick_number):
+    # # find number of multiples of pi/2
+    # N = int(np.round(2 * value / np.pi))
+    # if N == 0:
+    #     return "0"
+    # elif N == 1:
+    #     return r"$\pi/2$"
+    # elif N == 2:
+    #     return r"$\pi$"
+    # elif N % 2 > 0:
+    #     return r"${0}\pi/2$".format(N)
+    # else:
+    #     return r"${0}\pi$".format(N // 2)
+    # find number of multiples of pi/2
+    N = int(np.round(4 * value / np.pi))
+    if N == 0:
+        return "0"
+    elif N == 1:
+        return r"$\pi/4$"
+    elif N == 2:
+        return r"$\pi/2$"
+    elif N % 2 > 0:
+        return r"${0}\pi/2$".format(N)
+    else:
+        return r"${0}\pi$".format(N // 2)
+
+
+
+def find_nearest(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return array[idx]
+
+
+def find_nearestIdx(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return idx
+
+
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+print('---- Input parameters: -----')
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 2.56140350877193 #2.56140350877193, 4.0852130325814535
+# beta = 2.0  #5.0
+# theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+# theta = 0.075  # 0.5
+# mu1 = 10.0
+# rho1 = 1.0
+alpha = 10.0
+# beta = 40.0
+# theta = 0.125
+mu1 = 10.0
+rho1 = 1.0
+# alpha = 2.0
+beta = 2.0  #5.0
+# theta = 1.0/4.0
+theta = 1.0/8.0
+#
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10.0
+# beta = 2.0  #5.0
+# theta = 1.0/8.0
+# #
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 10.0
+# beta = 2.0
+# theta = 0.125
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 9.0
+# beta = 2.0
+# theta = 0.075
+
+#
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 4.0
+# beta = 10.0
+# theta = 1.0/4.0
+
+# alpha = 10   #1.05263158
+# beta = 40.0  #5.0
+# # theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+#
+#
+# alpha = 2.0
+# beta = 2.0 #5.0
+# theta = 1/4.0
+# rho1 = 1.0
+# mu1 = 10.0
+
+# InterestingParameterSet :
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+
+
+# theta = 0.124242
+# gamma = 0.75
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 3.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=200)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+# Gamma_Values = '0'
+
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+
+print('ratio(left) =', b2/b1)
+print('ratio(right) = ', q1/q3_star)
+
+# return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+print('angles:', angles)
+
+print('curvature:', curvature)
+
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+determinantVec = np.vectorize(determinant)
+
+detValues = determinantVec(q1,q2,muGammas,q12)
+print('detValues:', detValues)
+
+
+detZeroidx = find_nearestIdx(detValues, 0)
+print('idx where det nearest to zero', idx)
+gammaClose = Gamma_Values[detZeroidx]
+print('gammaClose:', gammaClose)
+
+
+# --- Convert to numpy array
+Gamma_Values = np.array(Gamma_Values)
+angles = np.array(angles)
+
+curvature = np.array(curvature)
+
+# ---------------- Create Plot -------------------
+# plt.figure()
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams["font.family"] = "serif"
+mpl.rcParams["font.size"] = "9"
+width = 6.28 *0.5
+height = width / 1.618
+fig = plt.figure()
+# ax = plt.axes((0.15,0.21 ,0.75,0.75))
+# ax = plt.axes((0.15,0.21 ,0.8,0.75))
+# ax = plt.axes((0.19,0.21 ,0.8,0.75))
+ax = plt.axes((0.15,0.18 ,0.8,0.75))
+ax.tick_params(axis='x',which='major', direction='out',pad=1)
+# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
+ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=1) # changed pad = distance to title to 1 here!
+ax.xaxis.set_major_locator(MultipleLocator(0.02))
+ax.xaxis.set_minor_locator(MultipleLocator(0.01))
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+ax.grid(True,which='major',axis='both',alpha=0.3)
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+ax.set(xlim=(1.760, 1.880))
+
+
+# # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# # plt.title(r'angle$-\gamma$-Plot')
+# plt.plot(Gamma_Values, angles)
+# plt.scatter(Gamma_Values, angles)
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+
+# f,ax=plt.subplots(1)
+
+
+# ax.plot(muGammas, angles)
+# ax.scatter(muGammas, angles)
+
+# ax.plot(muGammas, angles, 'royalblue', zorder=3, )
+ax.plot(muGammas, curvature, 'royalblue', zorder=3, )
+# ax.scatter(Gamma_Values, angles)
+ax.set_xlabel(r"$q_3(\gamma)$")
+# ax.set_xlabel(r"$\gamma$")
+# ax.set_ylabel(r"curvature $\kappa$")
+ax.set_title(r"curvature $\kappa$", fontsize=9, pad = 4)
+
+# plt.xlabel("$q_3$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# ax.grid(True)
+
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+
+# ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
+# ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+# ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25))
+
+
+# ax.yaxis.set_major_formatter(ticker.FuncFormatter(
+# lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0'))
+# ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi))
+
+# ---------------------------- show pi values ------------------------------------
+# ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+# ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+# ax.legend()
+# # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# # ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+# # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+# ---------------------------------------------------------------
+
+
+# ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+# ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+# labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+
+
+# OLD :
+# ax.set_yticks([0, np.pi/8, np.pi/4 ])
+# labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+# ax.set_yticklabels(labels)
+
+
+# Plot Gamma Value that is closest to q3_star
+# ax.axvline(x = gammaClose, color = 'midnightblue', linestyle = 'dashed', label='$\gamma^*$')
+# ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.2)
+# ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.2)
+
+
+ax.axvline(x = q3_star, color = 'midnightblue', linestyle = 'dashed', label='$q_3^*$')
+ax.axvline(x = q1, color = 'forestgreen', linestyle = 'dashed', label='$q_1$')
+ax.axvline(x = q2, color = 'darkorange', linestyle = 'dashed', label='$q_2$')
+
+
+# ax.legend(loc='upper center')
+ax.legend(loc='upper left', bbox_to_anchor=(0.55, 0.85))
+# ax.legend(loc='upper left')
+
+# plt.xlabel("$q_3(\gamma)$")
+# plt.xlabel("$\gamma$")
+# plt.ylabel("angle")
+# plt.legend(loc='upper center')
+
+fig.set_size_inches(width, height)
+fig.savefig('Plot-Curvature-q3.pdf')
+
+plt.show()
+
+
+
+
+# plt.figure()
+# plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+# plt.plot(muGammas, angles)
+# plt.scatter(muGammas, angles)
+# # plt.axis([0, 6, 0, 20])
+# # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.xlabel("$\mu_\gamma$")
+# plt.ylabel("angle")
+# plt.legend()
+# plt.show()
+#