x = 1:0.05:10;
y = 1:0.05:10;
z = 1:0.05:10;



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


s = size(X,1);


% A1 = repmat(1,size,3);
% A2 = repmat(2,size,3);
% A3 = rand(size,4)+ones(size,4);
% % A3 = repmat(2,10,4);


A1 = repmat(1,s,50);
A2 = repmat(2,s,50);
A3 = rand(s,81)+ones(s,81);
% A3 = repmat(2,10,4);


A  = [A1 A3 A2];
B  = repmat(A,1,1,s);


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



% INPUT Parameters
mu_1 = 1;
rho_1 = 1;
theta = 0.05; 
alpha = 25; 
beta= 20;

set_mu_gamma = 'q1';
% set_mu_gamma = 'q2';
print_output= false;

x = linspace(0,10,30);     %~alpha   % most interesting area for alpha >= 0 
% y = linspace(0.05,100,121);     %~beta
y = linspace(2,10,10);     %~beta %TEST
y = 2;   % 2D Phase Diagram
z = linspace(0.01,0.99,30);  %~theta              %theta = 0 , 1  -> parabolic case
[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);


angle_3D = cell2mat(angle_3D);
V_3D = cell2mat(V_3D);

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

vtkwrite('PhaseDiagramType.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',V_3D);



% 2D Phase Diagram:
% fix parameter (beta here) 
% beta= 5;
% 
% x = linspace(-50,50,121);     %~alpha
% z = linspace(0.01,0.99,121);  %~theta              %theta = 0 , 1  -> parabolic case
% [X,Z] = meshgrid(x,z);
% 
% [A,angle_2D,V_2D,~] = arrayfun(@(a,theta)classifyMIN(mu_1,rho_1,a,beta,theta,set_mu_gamma,print_output),X,Z,'UniformOutput', false);
% 
% angle_2D = cell2mat(angle_2D);
% V_2D = cell2mat(V_2D);
% 

% Y = ones(size(X))*beta;
% vtkwrite('PhaseDiagram2D.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',angle_2D);
% 
% vtkwrite('PhaseDiagramType2D.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',V_2D);