%% Phase Diagram

clear all 
clc


% INPUT Parameters
mu_1 = 1;
rho_1 = 1;
% theta= 0.3;
% theta = 3/4;
theta = 0.05; %interesting
% theta = 0.85; 
% theta = 0.35; 
% theta = 0.068; 
theta = 1/8; %interesting

% alpha = 0.1;  
% beta= 0.05;
% alpha = 15; 
% beta= 25;
alpha = 25; 
beta= 20;

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';      %(elliptic case)
% set_mu_gamma = 'q2';    %(hyperbolic-case)
%set_mu_gamma = 'm';   % mean of q1,q2  

% print_output = true;
print_output = false;


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


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

% x = linspace(-1,1,101);     %~alpha
% x = linspace(10,25,101);     %~alpha
% x = linspace(-40,40,101);     %~alpha
% x = linspace(-10,60,101);     %~alpha
% y = linspace(0.5,42,101);     %~beta
% z = linspace(0.05,0.95,50);  %~theta
% z = linspace(0.01,0.99,100);  %~theta

% x = linspace(-60,60,51);     %~alpha
% y = linspace(0.5,80,51);     %~beta
% z = linspace(0.05,0.95,51);  %~theta


x = linspace(-40,40,101);     %~alpha
y = linspace(0.01,30,101);     %~beta
z = linspace(0.01,0.99,101);  %~theta

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

% Plot Curvature:

[~,~,~,K] = arrayfun(@(a)classifyMIN(mu_1,rho_1,a,beta,theta,set_mu_gamma,print_output),x,'UniformOutput', false);
K = cell2mat(K);

% mark points for a change from alpha = 14 to alpha = 15
[~,~,~,K1] = classifyMIN(mu_1,rho_1,14,25,theta,set_mu_gamma,print_output);
[~,~,~,K2] = classifyMIN(mu_1,rho_1,15,25,theta,set_mu_gamma,print_output);

figure
plot(x,K)
% legend('Type 1', 'Type 2', 'Type 3')
title(['Curvature, theta=', num2str(theta) ])
xlabel('alpha')
ylabel('kappa')
% mark points for a change from alpha = 14 to alpha = 15
% hold on
% scatter(14,K1,'*','r');
% scatter(15,K2,'*','r');

[~,~,~,K] = arrayfun(@(b)classifyMIN(mu_1,rho_1,alpha,b,theta,set_mu_gamma,print_output),y,'UniformOutput', false);
K = cell2mat(K);

figure
plot(y,K)
% legend('Type 1', 'Type 2', 'Type 3')
title(['Curvature, theta=', num2str(theta) ])
xlabel('beta')
ylabel('kappa')



%make a Plot with dependency on theta!
[~,~,~,K] = arrayfun(@(t)classifyMIN(mu_1,rho_1,alpha,beta,t,set_mu_gamma,print_output),z,'UniformOutput', false);
K = cell2mat(K);
figure
plot(z,K)
title(['Curvature' ])
xlabel('theta')
ylabel('kappa')


% Make a surface plot
[X1,Y1] = meshgrid(x,y);


[~,~,~,K]= arrayfun(@(a,b)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X1,Y1,'UniformOutput', false);
K = cell2mat(K);
figure
% surf(X1,Y1,K,'EdgeColor','none')
surf(X1,Y1,K,'EdgeColor','interp')
% surf(X1,Y1,K,'FaceAlpha',0.5)
title(['Curvature, theta=', num2str(theta) ])
xlabel('alpha')
ylabel('beta')
zlabel('kappa')

% hold on 
% scatter3(14,25,K1,'*','w');
% scatter3(15,25,K2, '*','w');

%
% Plot Angle depending on Theta

t = linspace(0.01,0.99,501);  %~theta
alpha = -1;
beta = 30;


[~,Angle,~,~] = arrayfun(@(t)classifyMIN(mu_1,rho_1,alpha,beta,t,set_mu_gamma,print_output),t,'UniformOutput', false);
Angle = cell2mat(Angle);
figure
plot(t,Angle)
hold on 
alpha = 1;
[~,Angle,~,~] = arrayfun(@(t)classifyMIN(mu_1,rho_1,alpha,beta,t,set_mu_gamma,print_output),t,'UniformOutput', false);
Angle = cell2mat(Angle);
plot(t,Angle)
hold on 
alpha = 2;
[~,Angle,~,~] = arrayfun(@(t)classifyMIN(mu_1,rho_1,alpha,beta,t,set_mu_gamma,print_output),t,'UniformOutput', false);
Angle = cell2mat(Angle);
plot(t,Angle)
hold on 
alpha = 10;
[~,Angle,~,~] = arrayfun(@(t)classifyMIN(mu_1,rho_1,alpha,beta,t,set_mu_gamma,print_output),t,'UniformOutput', false);
Angle = cell2mat(Angle);
plot(t,Angle)
legend('alpha = -1', 'alpha = 1', 'alpha = 2', 'alpha = 10')
title(['Angle, alpha=', num2str(alpha), 'beta= ' num2str(beta) ])
xlabel('theta')
ylabel('Angle')
axis([0 1 -2 4])
% mark points for a change from alpha = 14 to alpha = 15
% hold on
% scatter(14,K1,'*','r');
% scatter(15,K2,'*','r');

%
%








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


[A_2D,angle_2D,V_2D,K] = 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);

% fprintf('angles that are not zero: %d', angle_3D(~angle_3D == 0 ) )



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

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);


% Structure 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);

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));

%%% 2D - Plot (fixed Theta)

cm_2D = redblue(3);
% 
% EdgeColor = 'black';
EdgeColor = 'none';



figure
scatter(X1T1,Y1T1,[], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor',EdgeColor);
colormap(cm_2D)
hold on
scatter(X1T2,Y1T2,[], 'MarkerFaceColor',[1 0  0], 'MarkerEdgeColor',EdgeColor);
hold on 
scatter(X1T3,Y1T3,[],'MarkerFaceColor',	[1 1 0], 'MarkerEdgeColor',EdgeColor);

% mark points for a change from alpha = 14 to alpha = 15
% hold on 
% scatter(14,25,'*','w')
% scatter(15,25,'*','w')

% scatter(X1T3,Y1T3,[], angle_2D, 'filled','MarkerEdgeColor','black');
colormap(cm_2D)
legend('Type 1', 'Type 2', 'Type 3')
% legend('Type 1', 'Type 2')
title(['Fixed Theta Plot: theta = ', num2str(theta) ])
xlabel('alpha')
ylabel('beta')
hold off 



%%% 3D - Plot


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);


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');
% scatter3(YT1,ZT1,XT1, [], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','none');
scatter3(YT1,ZT1,XT1, []);
colormap(cm);
hold on
scatter3(YT2,ZT2,XT2, [], 'MarkerFaceColor',[1 0  0],'MarkerEdgeColor', 'none');
hold on 
% scatter3(XT3,YT3,ZT3, [],'MarkerFaceColor',	[1 1 0],'MarkerEdgeColor', 'black');
% scatter3(YT3,ZT3,XT3, [],'MarkerFaceColor',	[1 1 0],'MarkerEdgeColor', 'none');
scatter3(YT3,ZT3,XT3, [], angle_3D, 'filled');
legend('Type 1', 'Type 2', 'Type 3')
title(['Classification of Minimizers, theta:', num2str(theta) ])
% xlabel('alpha')
% ylabel('beta')
% zlabel('theta')
xlabel('beta')
ylabel('theta')
zlabel('alpha')

% Export Data for Paraview:

% vtkwrite('PhaseDiagram.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type', B);