Update

src/CellScript.py

@@ -17,6 +17,11 @@ import sys
 
 # ----------------------------------------------------------------------------------------------------------------------------
 
+
+
+
+
+
 # ----- Setup Paths -----
 InputFile  = "/inputs/cellsolver.parset"
 OutputFile = "/outputs/output.txt"
@@ -50,14 +55,19 @@ print('----------------------------')
 
 
 
+print('RunCellProblem...')
+RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+
+print('Read effective quantities...')
+Q, B = ReadEffectiveQuantities()
+print('Q:', Q)
+print('B:', B)
+
 
-#
-# print('RunCellProblem...')
-# RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
 
-# TEST:
-print('Compare_Classification...')
-Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+# Compare symbolicMinimization with Classification 'ClassifyMin' :
+# print('Compare_Classification...')
+# Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
 
 
 
@@ -73,7 +83,7 @@ Inp = False
 Inp_T = True
 print('Run symbolic Minimization...')
 #Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath)
-G, angle, type, kappa = eng.symMinimization(Inp_T,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
+G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
 # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
 G = np.asarray(G) #cast Matlab Outout to numpy array
 # --- print Output ---

src/ClassifyMin.py

@@ -12,6 +12,24 @@ import sys
 # from subprocess import Popen, PIPE
 
 
+# --------------------------------------------------
+# 'classifyMin' classifies Minimizers by utilizing the result of
+# Lemma1.6
+#
+#
+#
+#
+# 'classifyMin_ana': (Special-Case : Lemma1.4)
+# ..additionally assumes Poisson-ratio=0 => q12==0
+#
+#
+#
+# Output : MinimizingMatrix, Angle, Type, Curvature
+
+
+
+
+
 def harmonicMean(mu_1, beta, theta):
     return mu_1*(beta/(theta+(1-theta)*beta))
 
@@ -55,6 +73,9 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
     b2 = prestrain_b2(rho_1, beta, alpha,theta)
     return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
 
+
+
+
 # --------------------------------------------------------------------
 # Classify Type of minimizer  1 = R1 , 2 = R2 , 3 = R3                          # before : destinction between which axis.. (4Types )
 # where
@@ -67,6 +88,13 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp
 # R3 = E2 U E3 U P1.1 U P2 U H
 # -------------------------------------------------------------------
 def classifyMin(q1, q2, q3, q12, b1, b2,  print_Cases=False, print_Output=False):   #ClassifyMin_hom?
+    # Assumption of Classification-Lemma1.6:
+    #  1. [b3 == 0]
+    #  2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+
+    # TODO: check if Q is orthotropic here - assert()
+
+
     if print_Output: print("Run ClassifyMin...")
     CaseCount = 0
     epsilon = sys.float_info.epsilon #Machine epsilon

src/PhaseDiagram.py

@@ -9,11 +9,13 @@ 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)
 
@@ -34,67 +36,65 @@ import time
 #
 # 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
+#
+
 
 
-# unabhängig von alpha...
-def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
-    # ------------------------------------ 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
-
-
-
-
-
-# --- SETUPS PATHS
+# ----------- SETUP PATHS
 # InputFile  = "/inputs/cellsolver.parset"
 # OutputFile = "/outputs/output.txt"
-
 InputFile  = "/inputs/computeMuGamma.parset"
 OutputFile = "/outputs/outputMuGamma.txt"
 # --------- Run  from src folder:
@@ -114,11 +114,14 @@ mu1 = 10.0               # TODO : here must be the same values as in the Parset
 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 = 'infinity'
 # gamma = 0.5
+# gamma = 0.25
+# gamma = 1.0
 
 print('---- Input parameters: -----')
 print('mu1: ', mu1)
@@ -138,12 +141,14 @@ print('----------------------------')
 # print_Cases = True
 # print_Output = True
 
+                        #TODO
+# generalCase = False #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
 
-make_3D_plot = True
+# make_3D_plot = True
 make_3D_PhaseDiagram = True
 make_2D_plot = False
 make_2D_PhaseDiagram = False
-# make_3D_plot = False
+make_3D_plot = False
 # make_3D_PhaseDiagram = False
 # make_2D_plot = True
 # make_2D_PhaseDiagram = True
@@ -183,9 +188,9 @@ make_2D_PhaseDiagram = False
 
 # ---------------------- 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 = 20 # Number of sample points in each direction
+# 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 = 10 # Number of sample points in each direction
 
 
@@ -216,23 +221,32 @@ if make_3D_PhaseDiagram:
     print('Written to VTK-File:', VTKOutputName )
 
 if make_2D_PhaseDiagram:
-    alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+    # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+    # thetas_ = np.linspace(0.01,0.99,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;
-    thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
-    # 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)
 
 
+    # if generalCase:              #TODO
+    #     classifyMinVec = np.vectorize(classifyMin)
+    #     GetCellOutputVec = np.vectorize(GetCellOutput)
+    #     Q, B = GetCellOutputVec(alpha,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+    #
+    #     print('type of Q:', type(Q))
+    #     print('Q:', Q)
+    #
+    # else:
+    classifyMin_anaVec = np.vectorize(classifyMin_ana)
     GetMuGammaVec = np.vectorize(GetMuGamma)
-    muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1)
+    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)
+    # print('Types:', Types)
     # Out = classifyMin_anaVec(alphas,betas,thetas)
     # T = Out[2]
     # --- Write to VTK
@@ -265,8 +279,13 @@ if(make_3D_plot or make_2D_plot):
     # print('Type 2 occured here:', np.where(T == 2))
 
 
-print(alphas_)
-print(betas_)
+# print(alphas_)
+# print(betas_)
+
+
+
+
+
 # ALTERNATIVE
 # colors = ("red", "green", "blue")
 # groups = ("Type 1", "Type2", "Type3")

src/Plotq3-Angle.py

@@ -9,74 +9,29 @@ 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
-# from subprocess import Popen, PIPE
-#import sys
-
-# unabhängig von alpha...
-def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
-    # ------------------------------------ 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
 
 
+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"
-# 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)
@@ -88,42 +43,29 @@ print("InputFilepath: ", InputFilePath)
 print("OutputFilepath: ", OutputFilePath)
 print("Path: ", path)
 
-
-#1. Define Inputs Gamma-Array..
-#2. for(i=0; i<length(array)) ..compute Q_hom, B_eff from Cell-Problem
-#3
-
-# matrix = np.loadtxt(path + 'Qmatrix.txt', usecols=range(3))
-# print(matrix)
-
-# Use Shell Commands:
-# subprocess.run('ls', shell=True)
-
-
-#---------------------------------------------------------------
-
-# -------------------------- Input Parameters --------------------
-mu1 = 10.0
-rho1 = 1.0
-alpha = 10   #1.05263158
-beta = 40.0  #5.0
-# theta = 1.0/4.0
-theta = 1.0/8.0  # 0.5
-
-# InterestingParameterSet :
+print('---- Input parameters: -----')
 # mu1 = 10.0
 # rho1 = 1.0
 # alpha = 10
 # beta = 40.0
-# theta = 1.0/8.0
+# theta = 0.125
+#
 
 
+mu1 = 10.0
+rho1 = 1.0
+# alpha = 10.02021333
+alpha = 10.0
+beta = 2.0
+theta = 0.125
+# 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('---- Input parameters: -----')
 print('mu1: ', mu1)
 print('rho1: ', rho1)
 print('alpha: ', alpha)
@@ -131,116 +73,90 @@ print('beta: ', beta)
 print('theta: ', theta)
 # print('gamma:', gamma)
 print('----------------------------')
-# ----------------------------------------------------------------
-
-
 
+# ----------------------------------------------------------------
 
 
-Gamma_Values = np.linspace(0.01, 5, num=15)    # TODO variable Input Parameters...alpha,beta...
+gamma_min = 0.01
+gamma_max = 1.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=100)    # TODO variable Input Parameters...alpha,beta...
 print('(Input) Gamma_Values:', Gamma_Values)
 # mu_gamma = []
 
 
-#
 # # --- Options
-# RUN = True
-# # RUN = False
 # # make_Plot = False
-make_Plot = True      # vll besser : Plot_muGamma
-#
-# if RUN:
-#     for gamma in Gamma_Values:
-#         print("Run Cell-Problem for Gamma = ", gamma)
-#         # print('gamma='+str(gamma))
-#         with open(InputFilePath, 'r') as file:
-#             filedata = file.read()
-#         filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
-#         f = open(InputFilePath,'w')
-#         f.write(filedata)
-#         f.close()
-#         # --- Run Cell-Problem
-#         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_gammaValue = float(s[0])
-#         print("mu_gamma:", mu_gammaValue)
-#         mu_gamma.append(mu_gammaValue)
-#     # ------------end of for-loop -----------------
-#     print("(Output) Values of mu_gamma: ", mu_gamma)
-# # ----------------end of if-statement -------------
-#
-# # mu_gamma=[2.06099, 1.90567, 1.905]
-# # mu_gamma=[2.08306, 1.905, 1.90482, 1.90479, 1.90478, 1.90477]
-#
-# ##Gamma_Values = np.linspace(0.01, 20, num=12) :
-# #mu_gamma= [2.08306, 1.91108, 1.90648, 1.90554, 1.90521, 1.90505, 1.90496, 1.90491, 1.90487, 1.90485, 1.90483, 1.90482]
-#
-# ##Gamma_Values = np.linspace(0.01, 2.5, num=12)
-# # mu_gamma=[2.08306, 2.01137, 1.96113, 1.93772, 1.92592, 1.91937, 1.91541, 1.91286, 1.91112, 1.90988, 1.90897, 1.90828]
-#
-# # Gamma_Values = np.linspace(0.01, 2.5, num=6)
-# # mu_gamma=[2.08306, 1.95497, 1.92287, 1.91375, 1.9101, 1.90828]
-
+make_Plot = True
 
+# Get values for mu_Gamma
 GetMuGammaVec = np.vectorize(GetMuGamma)
-muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1)
+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)
 
-classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
 
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
 G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
 # _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
 
 print('angles:', angles)
 
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+
+
+
 # Make Plot
 if make_Plot:
     plt.figure()
-    plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
-    plt.plot(muGammas, angles)
-    plt.scatter(muGammas, angles)
+    # 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 = 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.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^*$')
+
+
+    # Plot Gamma Value that is closest to q3_star
+    plt.axvline(x = gammaClose, color = 'b', linestyle = 'dashed', label='$\gamma^*$')
+
+    plt.axvspan(gamma_min, gammaClose, color='red', alpha=0.5)
+    plt.axvspan(gammaClose, gamma_max, color='green', alpha=0.5)
+
+    plt.xlabel("$\gamma$")
     plt.ylabel("angle")
     plt.legend()
     plt.show()
-
-#
-# # ------------- RUN Matlab symbolic minimization program
-# eng = matlab.engine.start_matlab()
-# # s = eng.genpath(path + '/Matlab-Programs')
-# s = eng.genpath(path)
-# eng.addpath(s, nargout=0)
-# # print('current Matlab folder:', eng.pwd(nargout=1))
-# eng.cd('Matlab-Programs', nargout=0)  #switch to Matlab-Programs folder
-# # print('current Matlab folder:', eng.pwd(nargout=1))
-# Inp = False
-# print('Run symbolic Minimization...')
-# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
-# # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
-# G = np.asarray(G) #cast Matlab Outout to numpy array
-#
-# # --- print Output ---
-# print('Minimizer G:')
-# print(G)
-# print('angle:', angle)
-# print('type:', type )
-# print('curvature:', kappa)