diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py
new file mode 100644
index 0000000000000000000000000000000000000000..87e9f935fb8773b2261ece1826eb6c727f4b6ce4
--- /dev/null
+++ b/src/ClassifyMin.py
@@ -0,0 +1,199 @@
+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 subprocess import Popen, PIPE
+
+
+
+def harmonicMean(mu_1,beta,theta):
+ return mu_1*(beta/(theta+(1-theta)*beta))
+
+
+def arithmeticMean(mu_1,beta,theta):
+ return mu_1*((1-theta)+theta*beta)
+
+
+def prestrain_b1(rho_1,beta,alpha):
+ return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)));
+
+def prestrain_b2(rho_1,beta,alpha):
+ return (3.0*rho_1/(4.0*((1.0-theta)+ theta*beta)))*(1-theta*(1+beta*alpha));
+
+
+# define function to be minimized
+def f(a1,a2,q1,q2,q3,q12,b1,b2):
+ A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+ print(A)
+ B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+ print(B)
+ a = np.array([a1, a2])
+ print(a)
+ tmp = np.dot(A,a)
+ print(tmp)
+ tmp2 = np.dot(a,tmp)
+ print(tmp2)
+ tmpB = np.dot(B,a)
+ print(tmpB)
+ # print(q1*(b1**2))
+ # print(q2*(b2**2))
+ # print(q12*b1*b2)
+ return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
+
+
+# Classify Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
+def classifyMin(q1, q2, q3, q12, b1, b2, print_output = "True"):
+ print("Run ClassifyMin...")
+ CaseCount = 0
+ epsilon = sys.float_info.epsilon
+
+ #TEST
+ # A = np.array([[1, 2], [3, 4]])
+ # kappa = 1.0
+ # angle = 1.57
+ # type = 2
+
+ determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
+ print("determinant:", determinant)
+
+
+ if abs(determinant) < epsilon:
+ print('parabolic case (determinant equal zero)')
+
+
+ if determinant >= epsilon:
+ print('elliptic case (determinant greater zero)')
+ a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+ a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+ prod = a1_star*a2_star
+
+ b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
+ b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
+
+ if prod >= epsilon:
+ print('(E1) - inside Lambda ')
+ a1 = a1_star
+ a2 = a2_star
+ type = 3 # non-axial Minimizer
+ CaseCount += 1
+ if abs(prod) < epsilon:
+ print('(E2) - on the boundary of Lambda ')
+ # if q2*b2**2 < q1*b1**2: # Needs to be updated to Mixed Term!!! just define f as a function and check value?!
+ #check
+ if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ a1 = b1_star
+ a2 = 0.0
+ type = 1
+ CaseCount += 1
+ if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ a1 = 0
+ a2 = b2_star
+ type = 2
+ CaseCount += 1
+ if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) = f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ # Two Minimizers pick one
+ a1 = b1_star
+ a2 = 0.0
+ type = 4
+ CaseCount += 1
+
+ if prod < -1.0*epsilon:
+ print('(E3) - Outside Lambda ')
+ if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ a1 = b1_star
+ a2 = 0.0
+ type = 1
+ CaseCount += 1
+ if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+ a1 = 0
+ a2 = b2_star
+ type = 2
+ CaseCount += 1
+ #TODO ...does type4 happen here?
+
+
+
+
+ return A, angle ,type ,kappa
+
+
+
+
+
+# ---------------------------------------------- Main ---------------------
+print('Running Python Code')
+
+
+# --- Input Parameters ----
+
+mu_1 = 1.0
+rho_1 = 1.0
+alpha = 14.0
+beta = 20.0
+theta = 1.0/4.0
+
+
+# define q1, q2 , mu_gamma, q12
+# 1. read from Cell-Output
+# 2. define values from analytic formulas (expect for mu_gamma)
+q1 = harmonicMean(mu_1,beta,theta)
+q2 = arithmeticMean(mu_1,beta,theta)
+q12 = 0.0 # (analytical example)
+
+# set mu_gamma to value or read from Cell-Output
+mu_gamma = q1 # TODO read from Cell-Output
+b1 = prestrain_b1(rho_1,beta,alpha)
+b2 = prestrain_b2(rho_1,beta,alpha)
+
+
+print("q1: ", q1)
+print("q2: ", q2)
+print("mu_gamma : ", mu_gamma)
+print("b1: ", b1)
+print("b2: ", b2)
+
+print("machine epsilon", sys.float_info.epsilon)
+
+## ------- options --------
+# print_output = false;
+
+
+
+
+
+
+
+A, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2)
+
+print("Matrix A")
+print(A)
+print("angle",angle)
+print("type",type)
+print("kappa", kappa)
+
+
+Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+print("Test", Test)
+
+
+
+# Run several times / make Array :
+# Gamma_Values = np.linspace(0.01, 2.5, num=12)
+# print(Gamma_Values)
+
+
+# Make Plot
+# plt.figure()
+# plt.title(r'$\mu_\gamma(\gamma)$-Plot') # USE MATHEMATICAL EXPRESSIONS IN TITLE
+# plt.plot(Gamma_Values, mu_gamma)
+# plt.scatter(Gamma_Values, mu_gamma)
+# plt.xlabel("$\gamma$")
+# plt.ylabel("$\mu_\gamma$")
+# plt.legend()
+# plt.show()