diff --git a/Matlab-Programs/symMinimization.m b/Matlab-Programs/symMinimization.m
new file mode 100644
index 0000000000000000000000000000000000000000..74ad4aca2241b62bb84fc7eeaa82983016eb8617
--- /dev/null
+++ b/Matlab-Programs/symMinimization.m
@@ -0,0 +1,367 @@
+
+function [G, angle, Type, kappa] = symMinimization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath) %(Q_hom,B_eff)
+
+
+syms v1 v2 q1 q2 q3 q12 b1 b2 b3
+
+
+% -------- Options ----------
+% % print_Input = false; %effective quantities
+% print_Input = true;
+% % print_statPoint = false;
+% print_statPoint = true;
+% print_Minimizer = true;
+% % make_FunctionPlot = false;
+% make_FunctionPlot = true;
+
+print_Uniqueness = false;
+check_StationaryPoints = false; % could be added to Input-Parameters..
+compare_with_Classification = false; %maybe added later ---
+
+
+%fprintf('Functions to be minimized:')
+f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2-2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
+ + q12*(v1^2*v2^2-b2*v1^2-b1*v2^2+b1*b2);
+f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
+ + q12*(v1^2*v2^2+b2*v1^2+b1*v2^2+b1*b2);
+
+
+% ---- Fix parameters
+
+if ~exist('InputPath','var')
+ % third parameter does not exist, so default it to something
+ absPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs";
+ end
+
+
+% 1. Import effective quantities from CellSolver-Code:
+%read as sparse Matrix...
+try %absolutePath
+ Qmat = spconvert(load(absPath + '' + "/QMatrix.txt"));
+ Bmat = spconvert(load(absPath + '' + "/BMatrix.txt"));
+% fprintf('Use absolute Path')
+catch ME % use relativePath
+ Qmat = spconvert(load('../outputs/QMatrix.txt'));
+ Bmat = spconvert(load('../outputs/BMatrix.txt'));
+% fprintf('Use relative Path')
+end
+%convert to full matrix...
+Qmat = full(Qmat);
+Bmat = full(Bmat);
+
+
+
+
+% --- TODO CHECK: assert if Q is orthotropic ???
+
+
+if print_Input
+ fprintf('effective quadratic form:')
+ Qmat
+ fprintf('effective prestrain')
+ Bmat
+
+ % check if Q is (close to..) symmetric
+ % könnte Anti-symmetric part berechnen und schauen dass dieser klein?
+ % Test: issymmetric(Qmat) does not work for float matrices?
+ % symmetric part 0.5*(Qmat+Qmat')
+ % anti-symmetric part 0.5*(Qmat-Qmat')
+ if norm(0.5*(Qmat-Qmat'),'fro') < 1e-10
+ fprintf('Qmat (close to) symmetric \n')
+ else
+ fprintf('Qmat not symmetric \n')
+ end
+
+
+ % Check if B_eff is diagonal this is equivalent to b3 == 0
+ if abs(Bmat(3)) < 1e-10
+ fprintf('B_eff is diagonal (b3 == 0) \n')
+ else
+ fprintf('B_eff is NOT diagonal (b3 != 0) \n')
+ end
+end
+
+
+
+
+% CAST VALUES TO SYM FIRST? This is done anyway..
+% % Substitute effective quantitites
+f_plus = subs(f_plus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
+ Bmat(1), Bmat(2), Bmat(3)});
+
+f_minus = subs(f_minus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
+ Bmat(1), Bmat(2), Bmat(3)});
+
+
+% Compute the Gradients
+df_plusx = diff(f_plus,v1);
+df_plusy = diff(f_plus,v2);
+df_minusx = diff(f_minus,v1);
+df_minusy = diff(f_minus,v2);
+
+% Setup Equations Grad(f) = 0
+eq1 = df_plusx == 0;
+eq2 = df_plusy == 0;
+eqns_plus = [eq1, eq2];
+eq3 = df_minusx == 0;
+eq4 = df_minusy == 0;
+eqns_minus = [eq3, eq4];
+
+% ------- Symbolically Solve Equations:
+% More robust (works even for values b_3 ~ 1e-08 ):
+S_plus = solve(eqns_plus,v1,v2,'MaxDegree' , 5);
+S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5);
+
+A_plus = S_plus.v1;
+B_plus = S_plus.v2;
+A_minus = S_minus.v1;
+B_minus = S_minus.v2;
+
+if check_StationaryPoints
+ %---------- TEST if Grad(f) = 0 ---------------------
+ fprintf('Testing equation grad(f) = 0 with stationary points')
+ for i = 1:size(A_plus,1)
+ fprintf('Testing %d.point (f_plus): ',i )
+ [ double(subs(subs(df_plusx,v1,A_plus(i)),v2,B_plus(i))), double(subs(subs(df_plusy,v1,A_plus(i)),v2,B_plus(i))) ]
+ end
+ for i = 1:size(A_minus,1)
+ fprintf('Testing %d.point (f_minus): ',i )
+ [double(subs(subs(df_minusx,v1,A_minus(i)),v2,B_minus(i))), double(subs(subs(df_minusy,v1,A_minus(i)),v2,B_minus(i)))]
+ end
+ % ------------------------------------
+end
+
+% --- Extract only Real-Solutions
+% fprintf('real stationary points of f_plus:')
+tmp1 = A_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
+tmp2 = B_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
+A_plus = tmp1;
+B_plus = tmp2;
+SP_plus = [A_plus,B_plus];
+
+% fprintf('real stationary points of f_minus:')
+tmp1 = A_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
+tmp2 = B_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
+A_minus = tmp1;
+B_minus = tmp2;
+SP_minus = [A_minus,B_minus];
+
+% TODO one should use f_plus.subs(A_plus..) to compute function value symbolically?
+% in the end only the stationaryPoints are used.. Ok to compare function values numerically
+
+% Determine global Minimizer from stationary points:
+% fprintf('function values at stationary points (f_plus):')
+T_plus = arrayfun(@(v1,v2) double(f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_plus,B_plus,'UniformOutput', false);
+T_plus = cell2mat(T_plus);
+
+%Test: use Substitution
+% subs(f_plus,{v1, v2}, {A_plus,B_plus})
+
+% fprintf('function values at stationary points (f_minus):')
+T_minus = arrayfun(@(v1,v2) double(f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_minus,B_minus,'UniformOutput', false);
+T_minus = cell2mat(T_minus);
+
+%Test: use Substitution
+% T_minus = subs(f_minus,{v1, v2}, {A_minus,B_minus})
+% double(T_minus)
+
+if print_statPoint
+ fprintf('real stationary points of f_plus: ')
+% SP_Plus %alternative: output as symbolic (can be unwieldy)
+ double(SP_plus)
+ fprintf('real stationary points of f_minus:')
+% SP_Minus %alternative: output as symbolic (can be unwieldy)
+ double(SP_minus)
+ fprintf('function values at stationary points (f_plus):')
+ T_plus
+ fprintf('function values at stationary points (f_minus):')
+ T_minus
+end
+
+% --- Find Stationary Point(s) with minimum Value
+[Min_plus,MinIdx_plus] = min(T_plus, [], 'all', 'linear'); %find one min...
+[Min_minus,MinIdx_minus] = min(T_minus, [], 'all', 'linear');
+% [Min_minus,MinIdx_minus] = min(T_minus) % works with symbolic too
+
+% Compare Minimizers of f_plus & f_minus
+[globalMinimizerValue,GlobalIdx] = min([Min_plus,Min_minus]);
+
+if GlobalIdx == 1 %Min_plus % i.e. Global Minimizer is given by f_plus
+ GlobalMinimizer = SP_plus(MinIdx_plus,:);
+ sign = 1.0;
+elseif GlobalIdx == 2 %Min_minus % i.e. Global Minimizer is given by f_minus
+ GlobalMinimizer = SP_minus(MinIdx_minus,:);
+ sign = -1.0;
+end
+
+% ------ Check if there are more SP with the same value...
+ MinIndices_minus = find(T_minus(:) == globalMinimizerValue); % Find indices of All Minima
+ MinIndices_plus = find(T_plus(:) == globalMinimizerValue); % Find indices of All Minima
+ numMinSP_minus = size(MinIndices_minus,1); % One of these is always >= 2 due to the structure of the roots..
+ numMinSP_plus = size(MinIndices_plus,1);
+
+% AllMinSP_minus = SP_minus(MinIndices_minus,:)
+% AllMinSP_minus = double(SP_minus(MinIndices_minus,:))
+% AllMin = T_minus(MinIndices) %bereits klar dass diese selben funktionswert haben..
+
+ Minimizer = sign*(GlobalMinimizer'*GlobalMinimizer); % global minimizing Matrix G*
+ MinimizerCount = 1;
+ % different Stationary Points might correspond to the same minimizing
+ % Matrix G*... check this:
+
+ % Compare only with other StationaryPoints/Minimizers
+ % remove Index of Minimizer
+ if GlobalIdx == 1
+ MinIndices_plus = MinIndices_plus(MinIndices_plus~=MinIdx_plus);
+ elseif GlobalIdx == 2
+ MinIndices_minus = MinIndices_minus(MinIndices_minus~=MinIdx_minus);
+ end
+
+ MinIndices = cat(1,MinIndices_plus,MinIndices_minus); %[Minimizers-Indices f_plus, Minimizer-Indices f_minus]
+
+ for i = 1:(numMinSP_minus+numMinSP_plus-1) % -1: dont count Minimizer itself..
+ idx = MinIndices(i);
+
+ if i > numMinSP_plus
+ SP = SP_minus(idx,:);
+ else
+ SP = SP_plus(idx,:);
+ end
+
+% SP_value = T_minus(idx) % not needed?
+ Matrix = sign*(SP'*SP);
+
+ if norm(double(Matrix-Minimizer),'fro') < eps %check is this sufficient here?
+ % both StationaryPoints correspond to the same
+ % (Matrix-)Minimizer
+ else
+ % StationaryPoint corresponds to a different (Matrix-)Minimizer
+ MinimizerCount = MinimizerCount + 1;
+ end
+ end
+% ----------------------------------------------------------------------------------------------------------------
+
+
+
+% Output Uniqueness of Minimizers:
+if print_Uniqueness
+ if MinimizerCount == 1
+ fprintf('Unique Minimzier')
+ elseif MinimizerCount == 2
+ fprintf('Two Minimziers')
+ else
+ fprintf('1-Parameter family of Minimziers')
+ end
+end
+
+
+
+
+
+% --- determine the angle of the Minimizer
+% a1 = Minimizer(1,1)
+% a2 = Minimizer(2,2)
+a1 = double(Minimizer(1,1));
+a2 = double(Minimizer(2,2));
+
+% compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
+e = [sqrt((a1/(a1+a2))), sqrt((a2/(a1+a2)))]; % always positive under sqrt here .. basically takes absolute value here
+angle = atan2(e(2), e(1));
+
+% compute curvature kappa
+kappa = (a1 + a2);
+
+% % CHeck off diagonal entries:
+% sqrt(a1*a2);
+% double(Minimizer);
+
+G = double(Minimizer);
+
+
+
+% --- "Classification" / Determine the TYPE of Minimizer by using
+% the number of solutions (Uniqueness?)
+% the angle (axial- or non-axial Minimizer)
+
+% (Alternative compute det[GlobalMinimizer' e1'] where e1 = [1 0] ?)
+
+% Check Uniqueness -- Options: unique/twoMinimizers/1-ParameterFamily
+if MinimizerCount == 1
+% fprintf('Unique Minimzier')
+ % Check if Minimizer is axial or non-axial:
+ if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+ Type = 3;
+ else % non-axial Minimizer
+ Type = 1;
+ end
+elseif MinimizerCount == 2
+% fprintf('Two Minimziers')
+ % Check if Minimizer is axial or non-axial:
+ if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+ Type = 3;
+ else % non-axial Minimizer
+ fprintf('ERROR: Two non-axial Minimizers cannot happen!')
+ end
+else
+% fprintf('1-Parameter family of Minimziers')
+ % Check if Minimizer is axial or non-axial:
+ if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
+% fprintf('ERROR: axial Minimizers cannot happen for 1-Parameter Family!')
+ else % non-axial Minimizer
+ Type = 2;
+ end
+end
+% ------------------------------------------------------------------------------------------------------
+
+
+
+
+if print_Output
+ fprintf(' --------- Output symMinimization --------')
+ fprintf('Global Minimizer v: (%d,%d) \n', GlobalMinimizer(1),GlobalMinimizer(2) )
+ fprintf('Global Minimizer Value f(v): %d \n', sym(globalMinimizerValue) ) %cast to symbolic
+ % fprintf('Global Minimizer Value : %d', globalMinimizerValue )
+
+ fprintf('Global Minimizer G: \n' )
+ G
+ fprintf("Angle = %d \n", angle)
+ fprintf("Curvature = %d \n", kappa)
+ fprintf("Type = %i \n", Type)
+ fprintf(' --------- -------------------- --------')
+end
+
+
+
+if make_FunctionPlot
+ fsurf(@(x,y) f_plus(x,y,q1,q2,q3,q12,b1,b2,b3)) % Plot functions
+ hold on
+ plot3(double(A_plus),double(B_plus),T_plus,'g*')
+ %Plot GlobalMinimizer:
+ hold on
+ plot3(double(GlobalMinimizer(1)),double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
+ % view(90,0)
+ % view(2)
+
+ figure
+ fsurf(@(x,y) f_minus(x,y,q1,q2,q3,q12,b1,b2,b3))
+ hold on
+ plot3(double(A_minus),double(B_minus),T_minus,'g*')
+ hold on
+ plot3(double(GlobalMinimizer(1)), double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
+end
+
+
+return
+
+
+% Write symbolic solution to txt-File in Latex format
+% fileID = fopen('txt.txt','w');
+% fprintf(fileID,'%s' , latex(S_plus.v1));
+% fclose(fileID);
+
+
+
+
+
+
+