%% Phase Diagram

clear all 
clc


%% Data-Storage
% Storage.mesh         = cell(maxlevel-minlevel+1,1);
% Storage.elemCount    = cell(maxlevel-minlevel+1,1);
% Storage.pFes         = cell(maxlevel-minlevel+1,1);


% INPUT Parameters
mu_1 = 1;
rho_1 = 1;
theta= 0.3;
% theta= 0.1;
% alpha = 2; 
% beta= 2;


fprintf('===========================================================================================================================');
fprintf(' Possible cases for global minimizers: (I) , (II) , (III) , (IV) ')
fprintf('============================================== INPUT PARAMETERS ===========================================================');
fprintf('mu_1:       %d              rho:1: %d                 theta: %d',  mu_1,rho_1,theta );
fprintf('===========================================================================================================================');


% choose u_gamma = q3 to be either q1, q2 or something in between 
%set_mu_gamma = 'q1';
set_mu_gamma = 'q2'; %(only hyperbolic-case)
% set_mu_gamma = 'm'



% print_output = true;
print_output = false;

% Test for fixed value
% [A,type] = classifyMIN(mu_1,rho_1,alpha,beta,theta,set_mu_gamma,print_output);

% create function handle for classifyMIN
% f = @(a,b) classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output);



% PLOT 
x = linspace(-20,20,45);     %~alpha
y = linspace(1.5,40,45);     %~beta
z = linspace(0.05,0.95,45);  %~theta

[X1,Y1] = meshgrid(x,y);


[A_2D,angle_2D,V_2D] = arrayfun(@(a,b)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X1,Y1,'UniformOutput', false);

[X,Y,Z] = meshgrid(x,y,z);

[A,angle_3D,V_3D] = arrayfun(@(a,b,theta)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X,Y,Z,'UniformOutput', false);

color_3D = cell2mat(V_3D);
color_2D = cell2mat(V_2D);

angle_2D = cell2mat(angle_2D);
angle_3D = cell2mat(angle_3D);

A = cell2mat(A);
A_2D = cell2mat(A_2D);

% cellfun(@(a) atan2(sqrt(a(2)),sqrt(a(1))), A_2D, 'UniformOutput',false);

X1 = reshape(X1,[],1);
Y1 = reshape(Y1,[],1);

X = reshape(X,[],1);
Y = reshape(Y,[],1);
Z = reshape(Z,[],1);
color_2D = reshape(color_2D,[],1);
color_3D = reshape(color_3D,[],1);


angle_2D = reshape(angle_2D,[],1);
angle_3D = reshape(angle_3D,[],1);


% scatter(X,Y,[],color,'filled'); %2D
% scatter(X1,Y1,[],color_2D,'filled', 'MarkerEdgeColor','black');
% colormap([1 0 0; 0 1 0; 0 0 1;  1 1 0;]);
% title('Fixed Theta Plot')
 
% scatter3(X,Y,Z,[],color_3D,'filled')  %3D
% colormap([1 0 0; 0 1 0; 0 0 1;  1 1 0;]);
% set(gca, 'clim', [0.5 4.5]);
% title('Classification of Minimizers, theta:')
% xlabel('alpha')
% ylabel('beta')
% zlabel('theta')
% view(90,0)
% view(2)



% Divide result depending on Type/Color
% V_2D = reshape(V_2D,[],1);


V_2DT1 = (color_2D == 1);
V_2DT2 = (color_2D == 2);
V_2DT3 = (color_2D == 3);

% [X1,Y1] = meshgrid(x,y);

V_2DT1 = reshape(V_2DT1,[],1);
V_2DT2 = reshape(V_2DT2,[],1);
V_2DT3 = reshape(V_2DT3,[],1);

X1T1 = V_2DT1.*X1;
Y1T1 = V_2DT1.*Y1;
X1T2 = V_2DT2.*X1;
Y1T2 = V_2DT2.*Y1;
X1T3 = V_2DT3.*X1;
Y1T3 = V_2DT3.*Y1;


%Filter out artificial zeros
% X1T1 = X1T1(Y1T1~=0);
% Y1T1 = Y1T1(Y1T1~=0);
% X1T2 = X1T2(Y1T2~=0);
% Y1T2 = Y1T2(Y1T2~=0);
% X1T3 = X1T3(Y1T3~=0);
% Y1T3 = Y1T3(Y1T3~=0);

% color1_2D = reshape(V_2DT1.*color_2D,[],1);
% color2_2D = reshape(V_2DT2.*color_2D,[],1);
% color3_2D = reshape(V_2DT3.*color_2D,[],1);

%ANGLE BETWEEN VECTOR and x-axis
% atan2(sqrt(a2),sqrt(a1));

% angle_2D = linspace(-3,3,length(X1T1));
cm_2D = redblue(3);

scatter(X1T1,Y1T1,[], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm_2D)
hold on
scatter(X1T2,Y1T2,[], 'MarkerFaceColor',[1 0  0], 'MarkerEdgeColor','black');
hold on 
% scatter(X1T3,Y1T3,[],'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
scatter(X1T3,Y1T3,[], angle_2D, 'filled','MarkerEdgeColor','black');
colormap(cm_2D)
legend('Type 1', 'Type 2', 'Type 3')
title('Fixed Theta Plot')
xlabel('alpha')
ylabel('beta')
hold off 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3D 


V_3DT1 = (color_3D == 1);
V_3DT2 = (color_3D == 2);
V_3DT3 = (color_3D == 3);

% [X1,Y1] = meshgrid(x,y);

V_3DT1 = reshape(V_3DT1,[],1);
V_3DT2 = reshape(V_3DT2,[],1);
V_3DT3 = reshape(V_3DT3,[],1);

XT1 = V_3DT1.*X;
YT1 = V_3DT1.*Y;
ZT1 = V_3DT1.*Z;
XT2 = V_3DT2.*X;
YT2 = V_3DT2.*Y;
ZT2 = V_3DT2.*Z;
XT3 = V_3DT3.*X;
YT3 = V_3DT3.*Y;
ZT3 = V_3DT3.*Z;


%Filter out artificial zeros
% XT1 = XT1(YT1~=0);
% YT1 = YT1(YT1~=0);
% ZT1 = ZT1(YT1~=0);
% XT2 = XT2(YT2~=0);
% YT2 = YT2(YT2~=0);
% ZT2 = ZT2(YT2~=0);
% XT3 = XT3(YT3~=0);
% YT3 = YT3(YT3~=0);
% ZT3 = ZT3(YT3~=0);




% color1 = reshape(V_3DT1.*color_3D,[],1);
% color2 = reshape(V_3DT2.*color_3D,[],1);
% color3 = reshape(V_3DT3.*color_3D,[],1);




%TODO: ANGLE between (e,e1)
% angle = linspace(-3,3,length(XT1)); 

% CosTheta = max(min(dot(u,v)/(norm(u)*norm(v)),1),-1);
% ThetaInDegrees = real(acosd(CosTheta));
%ANGLE BETWEEN VECTOR and x-axis
% atan2(vector.y,vector.x);



cm = redblue(90);
figure
%fixed Color
% scatter3(XT1,YT1,ZT1, [], 'MarkerFaceColor',[0.75 0 0],'MarkerEdgeColor', 'none');
%variing Color
scatter3(YT1,ZT1,XT1, [], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm);
hold on
scatter3(YT2,ZT2,XT2, [], 'MarkerFaceColor',[1 0  0],'MarkerEdgeColor', 'none');
hold on 
% scatter3(XT3,YT3,ZT3, [],'MarkerFaceColor',[0 0 1 ],'MarkerEdgeColor', 'none');
scatter3(YT3,ZT3,XT3, [], angle_3D, 'filled');
legend('Type 1', 'Type 2', 'Type 3')
title('Classification of Minimizers, theta:')
% xlabel('alpha')
% ylabel('beta')
% zlabel('theta')

xlabel('beta')
ylabel('theta')
zlabel('alpha')















%
%
%

%  %% STORE DATA
%  format shortE
%  Storage.mesh{level}           = m;
%  
%  
%  
%  

% %Write Data to LATEX-Code

% writeFile = 'Table1.tex';
% % writeFile = strcat('Table',num2str(i),'.tex');
% C = table2cell(T);
% Head = ["" , "$\lVert u-u_h \rVert_0 $ ","","$\vert u-u_h \lvert_1 $","","$\lVert p-p_h \rVert_0$", "" ,...
%     "$\lVert j_hp-p_h \rVert_0$","", "$\lVert \operatorname{div}u_h \rVert_0$","$\lVert \operatorname{div}\Pi_h u_h \rVert_0$", ...
%     "$\lvert \lvert \lvert u-u_h \rvert \rvert \rvert$",""; ...
%     "Mesh", "error", "ord" "error", "ord","error", "ord" "error", "ord","error", "error", "error", "ord"];
% % C = [T.Properties.VariableNames;C];
% C = [ Head ;C ];
% C = cellfun(@num2str, C, 'UniformOutput',false);
% % Cell2LatexTable(C,  'Compile', false, 'File', 'Table.tex', 'hline', '3lines', 'vline','full', 'Caption', 'Flow-Problem','Alignment','Left')
% 
% input.tableColumnAlignment = 'ccccccccccccc';
% 
% 
% input.data = C;
% % input.tableColumnAlignment = 'c' 
% input.tableBorders = 1;
% % input.tableCaption = 'Flow';
% % generate LaTex code
% latex = latexTable(input);
% 
% 
% % save LaTex code as file
% fid=fopen(writeFile,'w');
% [nrows,ncols] = size(latex);
% 
% fprintf(fid,'\\begin{table}[h] \n ');
% fprintf(fid,'\\resizebox{1\\textwidth}{!}{ \\begin{minipage}{\\textwidth} ');
%  
% fprintf(fid,'\\begin{tabular}{l c l c l c} \n' );
% switch PROBLEM
%     case 0
%       fprintf(fid,' \\text{Problem : } &  No Flow  ');
%     case 1
%       fprintf(fid,' \\text{Problem : } &  Flow  ');
%     case 2
%       fprintf(fid,' \\text{Problem : } &  Potential-Flow ');
% end
% 
% fprintf(fid,'&\\text{Mesh-Distortion: }&  %d  ', mDistortion);
% fprintf(fid,' & \\text{Reconstruction - Force: } & %d\n \\\\ ', modRhs);
% fprintf(fid,'\\text{order: }&  %d ', order);
% if meshType
%     fprintf(fid,' &\\text{Mesh-Type:  }&  $\\mathbf{\\bigtriangleup}$ ' );
% else
%     fprintf(fid,' &\\text{Mesh-Type:  } & $\\mathbf{\\square}$  ' );
% end
% fprintf(fid,'&\\text{Reconstruction - Reaction: } &  %d\n \\\\ ', modLhs);
% fprintf(fid,' \\text{viscosity $\\nu $: } & %s', nu);
% fprintf(fid, ' & \\text{Reaction $ \\alpha $: } & %s ', alpha);
% if meshType
%      fprintf(fid,' & \\text{Elements: } & $ (\\mathbf{P+}_{k},P_{k-1}^{\\text{disc}},\\mathbf{RT}_{[k]})$   \n \\\\ '); 
% else
%      fprintf(fid,' &\\text{Elements: } & $(\\mathbf{Q}_{k}^2,P_{k-1}^{\\text{disc}},\\mathbf{BDM}_{[k]})$  \n \\\\ ');   
% end
% fprintf(fid,'\\hline \\hline ');
% fprintf(fid,' \\end{tabular} \\vspace{3mm}    \\end{minipage}} \n  ');
% 
% for row = 1:nrows
%     fprintf(fid,'%s\n',latex{row,:});
% end
% 
% fclose(fid);
% fprintf(strcat('\n... your LaTex code has been saved as  \n ',  writeFile ,' in your working directory\n'));
 
 