classifyMIN.m

function [A, angle, type, kappa] = classifyMIN (mu_1,rho_1,a,b,t,set_mu_gamma,print_output)

% returns
%   A : Matrix of basis coefficients [a1,a2,a3]
%
%   type : 
%   Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
%



type = 0; % either 1,2,3,4 


mu_h = @(b,t) mu_1.*(b./(t+(1-t).*b));  % harmonic mean
mu_bar = @(b,t) mu_1.*((1-t)+t.*b);     % mu_bar

if (set_mu_gamma == 'q1')
  mu_gamma = @(b,t) mu_h(b,t);
end
if (set_mu_gamma == 'q2')
  mu_gamma = @(b,t) mu_bar(b,t);
end
if (set_mu_gamma == 'm')
  mu_gamma = @(b,t) 0.5*(mu_h(b,t) + mu_bar(b,t));
end

%  q1 q2 q3..
q1 = mu_h(b,t);
q2 = mu_bar(b,t); 
q3 = mu_gamma(b,t);




% values for q1,q2,q3 should be positiv
% assert((q1 > 0 ) & (q2 > 0 ) & (q3 > 0), 'At least one of q1,q2 or q3 is not positive' )


% Compute components of B_eff (old)
% b1 = (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
% b2 =  mu_1.*(rho_1/8).*(1-t.*(1+b.*a));

%TEST (new)
b1 =  (3*rho_1/2).*b.*(1-t.*(1+a));
b2 =  (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a));


% H = [q1 q3; q3 q2];         
%check condition of H first
% fprintf('condition number of Matrix H: %d \n', cond(H));


CaseCount = 0; %check if more than one case ever happens

epsilon = eps;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if abs(det(A)) < epsilon * min(abs(det(A)),0)  
if abs(q1*q2-q3^2) < epsilon     
    
    fprintf('determinant equal zero (parabolic case)')
    fprintf('SHOULD NOT HAPPEN')
    
   
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ELLIPTIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (q1*q2-q3^2 > epsilon)
%     fprintf('determinant greater than zero (elliptic case)');
    
%     a1_star =  (q2(b,t).*q1(b,t).*b1(a,b,t) - q3(b,t).*q2(b,t).*b2(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
%     a2_star =  (q1(b,t).*q2(b,t).*b2(a,b,t) - q3(b,t).*q1(b,t).*b1(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
    a1_star =  (q2.*q1.*b1 - q3.*q2.*b2)./(q1.*q2-q3.^2); 
    a2_star =  (q1.*q2.*b2 - q3.*q1.*b1)./(q1.*q2-q3.^2);
    
    prod = a1_star*a2_star;
    
    
    if(prod > epsilon) % (E1)  inside Lamba             %
        % (a1_star,a2_star) is unique minimizer lies inside Lambda 
        % therefore Minimizer not aligned with axes
        
%         fprintf('\n elliptic-case: (E1)');
        a1 = a1_star;
        a2 = a2_star;
        type = 3;
        CaseCount = CaseCount + 1;
    end
    % Make distinction between boundary & outside (prod < 0 ) 
    if(abs(prod) < epsilon)  % (E2) on boundary of Lambda 
%         fprintf('\n elliptic-case: (E2)');
          
        % Uniqueness of gloal minimizer if lies on boundary  if prod = 0 
        % ----- % 
        

        % global minimizer lies on the boundary of Lambda depending on
        % condition:
        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
            a1 = b1;
            a2 = 0*b1;
            type = 1;  % Minimizer aligned with x1-axis
            CaseCount = CaseCount + 1;
        end
        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
            a1 = 0*b1;
            a2 = b2;
            type = 2; % Minimizer aligned with x2-axis
            CaseCount = CaseCount + 1;
        end
        if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
            % two Minimizers ..pick one
            a1 = b1;
            a2 = 0*b1;
            type = 4;
            CaseCount = CaseCount + 1;
        end
    end
    if((prod) < -1*epsilon) %Outside of Lambda 
        
        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
            a1 = b1;
            a2 = 0*b1;
            type = 1;  % Minimizer aligned with x1-axis
            CaseCount = CaseCount + 1;
        end
        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
            a1 = 0*b1;
            a2 = b2;
            type = 2; % Minimizer aligned with x2-axis
            CaseCount = CaseCount + 1;
        end
    end
    
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HYPERBOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (q1*q2-q3^2 < -1*epsilon) 
%     fprintf('determinant less than zero (hyperbolic case)');

    if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
        a1 = b1;
        a2 = 0*b1;
        type = 1;  % Minimizer aligned with x1-axis
        CaseCount = CaseCount + 1;
    end
    if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
        a1 = 0*b1;
        a2 = b2;
        type = 2; % Minimizer aligned with x2-axis
        CaseCount = CaseCount + 1;
    end
    if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
        % two Minimizers ..pick one
        a1 = b1;
        a2 = 0*b1;
        type = 4;
        CaseCount = CaseCount + 1;
    end

    
    % CAN NOT BE TYPE 3!!
    
end



% Compute a3 from a1 % a2
a3 = sqrt(2*a1*a2);                    % WRONG ? 




% compute angle between [sqrt(a1) , sqrt(a2)] and e1:
% angle = atan2(sqrt(a2),sqrt(a1));

v = [sqrt(a1), sqrt(a2)];
e1 = [1 0];

if (type == 3 )
%    angle = atan2(a2,a1);
%    angle = acos(v*e1/(sqrt(a1+a2)));
else
   angle = 0;
end
% angle = atan2(norm(cross(a,b)), dot(a,b))

% Try to compute angle of Matrices

E1 = [1 0;0 0 ];
% V = [a1 sqrt(a1*a2)/sqrt(2); sqrt(a1*a2)/sqrt(2) a2]; % This ?? NO

V = [a1 sqrt(a1*a2); sqrt(a1*a2) a2];         % This!

% CHECK
% U = trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1)))   
% if abs(U) > 1
%    fprintf('value greater 1') 
% end


% angle = atan2(sqrt(abs(a2)), sqrt(abs(a1)));   % does this make sense? 



% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1))));    %angle in radians
% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1))))/pi * 180; % angle in degrees


% CHECK does case (0,-b) ever happen?  Yes 
% if  (q2*b2^2 > q1*b1^2 && b2 < 0)
%     fprintf('point lies on negative half of y-axis'); 
% end



% CHECK if Minimizer ever lies inside lambda ... YES 
% if  (a1 ~= 0 && a2 ~= 0)
%     fprintf('minimizer lies inside lambda'); 
% end


% Alternative atan2d(y,x): returns angle in degrees 


% CHECK
e = [sqrt(a1), sqrt(a2)];        % Might be complex...
norm = sqrt(a1+a2);              % Might be complex...
e = e./norm;
angle = atan2(e(2), e(1));   %TEST

if (imag(e(1))~= 0 || imag(e(2))~=0)
   fprintf('complex vector e'); 
end

% if (imag(norm)~=0)
%    fprintf('complex norm'); % happens a lot .. 
% end


% do it this way to avoid complex numbers: 
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)); 


if (imag(e(1))~= 0 || imag(e(2))~=0)
   fprintf('complex vector e'); 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST  -------------------
% the magnitude kappa of the vector does not matter for the angle ..
% In order to also get negative angles just use:
% atan2(a2,a1) 

% angle = atan2(a2,a1);         %% FEHLER : muss [sqrt(a1) sqrt(a2)] betrachten..
% angle = atan2(sqrt(a2),sqrt(a1));  % Inputs must be real... 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% angle2 = acos(e(1)/(sqrt(e(1)^2 + e(2)^2) ));
% if angle ~= angle2
%    interesting = 1; 
% end

% if (angle ~= 0 | angle ~= 3.1416)
%    fprintf('angle not zero or pi.. plot type')
%    type  
%    angle
% end

%compute Kappa? 
% fprintf('Output Kappa:')
% k = sqrt(abs(a1) + abs(a2));  % ? 
kappa = (a1 + a2); 



% e = [sqrt(a1), sqrt(a2)];
% norm = sqrt((a1+a2));
% e = e./norm;
% K.*(e'*e);


if (CaseCount > 1)
  fprintf('=============================== ERROR========================================= \n')
  fprintf('Error: More than one Case happened!')
end



% Coefficients of minimizer 

if(print_output)
   fprintf(' \n ') 
   fprintf('=============================== OUTPUT ========================================= \n')
   fprintf('Minimizer is of Type: %d  \n' , type);
   fprintf('Coefficients a1,a2,a3 given by : %d, %d, %d   \n', a1, a2, a3);
   fprintf('================================================================================ \n')
end


A = [a1, a2, a3];

end