%% Results of a scan to check the critical kT
% We set kTcrit as the value where the heat flux can increase by two order
% of magnitude on 50 t.u. (i.e. from 1e-20 to 1e-18) with an phi backg init
% of amplitude 1e-5 and on the mode ky=0.2. When the last heat flux value
% is larger than 1e-15 we fix it as kTmin. Then we decrease kT by 0.1 and
% rerun. If 1e-15 is not reached, this value is kTcrit. Otherwise iterate.
% Note: this criterium is equivalent to gamma<0.05
% P J kTcrit
kTcritP0 = [...
0 1 3.5;...    
0 2 2.5;...    
0 3 2.3;...    
0 4 2.4;...    
0 5 2.5;...    
0 6 2.7;...    
0 7 2.9;...    
0 8 2.9;...    
];
kTcritP1 = [... %used ky=0.3
1 1 4.4;...    
1 2 3.5;...    
1 3 3.5;...    
1 4 3.8;...    
1 5 4.4;...    
% 1 6 0;...    
% 1 7 0;...    
% 1 8 0;...    
];
kTcritP2 = [...
% 2 1 0;...    
% 2 2 0;...    
% 2 3 0;...    
% 2 4 0;...    
% 2 5 0;...    
% 2 6 0;...    
% 2 7 0;...    
% 2 8 0;...    
];
kTcritP3 = [...
% 3 1 0;...    
% 3 2 0;...    
% 3 3 0;...    
% 3 4 0;...    
% 3 5 0;...    
% 3 6 0;...    
% 3 7 0;...    
% 3 8 0;...    
];
kTcritP4 = [...
4 1 3.5;...    
4 2 2.8;...    
4 3 2.8;...    
4 4 3.0;...    
4 5 4.2;...    
4 6 3.7;...    
4 7 3.4;...    
4 8 3.3;...    
];
kTcritP5 = [...
5 1 3.8;...    
5 2 2.9;...    
5 3 3.0;...    
5 4 3.6;...    
5 5 4.2;...    
5 6 3.9;...    
5 7 3.6;...    
5 8 3.5;...    
];
kTcritP6 = [...
6 1 3.7;...    
6 2 2.9;...    
6 3 2.9;...    
6 4 3.1;...    
6 5 4.2;...    
6 6 3.8;...    
6 7 3.6;...    
6 8 3.4;...    
];
kTcritP7 = [...
7 1 4.0;...    
7 2 3.1;...    
7 3 3.2;...    
7 4 3.7;...    
7 5 4.4;...    
7 6 4.0;...    
7 7 3.7;...    
7 8 3.6;...    
];
kTcritP8 = [...
8 1 4.0;...    
8 2 3.1;...    
8 3 3.0;...    
8 4 3.3;...    
8 5 4.4;...    
8 6 4.0;...    
8 7 3.7;...    
8 8 3.5;...    
];
kTcritP12 = [...
12 2 3.5;...    
12 4 3.9;...    
12 6 4.3;...    
12 8 3.8;...    
12 10 3.7;...    
12 12 4.0;...    
];

kT_tot = [kTcritP0; kTcritP1; kTcritP2; kTcritP3; kTcritP4; kTcritP5; kTcritP6; kTcritP7; kTcritP8; kTcritP12];
% Stability map
Np = 12;
kT_map = zeros(Np);
[XX,YY] = meshgrid(0:Np,0:Np);
sz = size(kT_tot);
for i = 1:sz(1)
  kT_map(kT_tot(i,1)+1,kT_tot(i,2)+1) = kT_tot(i,3)-4;
end
%plots
figure
subplot(121)
plot(kTcritP0(:,2),kTcritP0(:,3),'--o','DisplayName',num2str(kTcritP0(1,1))); hold on;
plot(kTcritP4(:,2),kTcritP4(:,3),'--o','DisplayName',num2str(kTcritP4(1,1))); hold on;
plot(kTcritP5(:,2),kTcritP5(:,3),'--o','DisplayName',num2str(kTcritP5(1,1))); hold on;
plot(kTcritP6(:,2),kTcritP6(:,3),'--o','DisplayName',num2str(kTcritP6(1,1)));
plot(kTcritP7(:,2),kTcritP7(:,3),'--o','DisplayName',num2str(kTcritP7(1,1)));
plot(kTcritP8(:,2),kTcritP8(:,3),'--o','DisplayName',num2str(kTcritP8(1,1)));
plot(kTcritP12(:,2),kTcritP12(:,3),'--o','DisplayName',num2str(kTcritP12(1,1)));
ylim([0 5]);
ylabel('$\kappa_T^{crit}$')
xlabel('$J$');
subplot(122)
imagesc_custom(XX,YY,kT_map);
colormap(bluewhitered)
xlabel('$J$');
ylabel('$P$');
% bubblechart(kTcritP0(:,1),kTcritP0(:,2),kTcritP0(:,3)); hold on
% bubblechart(kTcritP4(:,1),kTcritP4(:,2),kTcritP4(:,3)); hold on
% bubblechart(kTcritP5(:,1),kTcritP5(:,2),kTcritP5(:,3))
% bubblechart(kTcritP6(:,1),kTcritP6(:,2),kTcritP6(:,3))
% bubblechart(kTcritP7(:,1),kTcritP7(:,2),kTcritP7(:,3))
% bubblechart(kTcritP8(:,1),kTcritP8(:,2),kTcritP8(:,3))
% bubblechart(kTcritP12(:,1),kTcritP12(:,2),kTcritP12(:,3))