%figure;scatter(mdsi_na(:,100),mdsi_sa(:,100),'jitter','on','jitterAmount',0.1);xlim([0 100]);ylim([0 100])



%===============
%   JOINT PDF 
%===============

mna1=mdsi_na(:,1:100);
msa1=mdsi_sa(:,1:100);

X1=cat(2,mna1(:),msa1(:));

%bwn=[0.0422,0.0097];
%[f1,x1] = ksdensity(X1,'bandwidth',bwn);
[f1,x1] = ksdensity(X1);

xl=[0,100];yl=[0,100];
%[xq,yq] = meshgrid(linspace(-4.75,4.75),linspace(-1,1));
[xq,yq] = meshgrid(linspace(xl(1),xl(2),500),linspace(yl(1),yl(2),500));
z1 = griddata(x1(:,1),x1(:,2),f1,xq,yq);
%z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);

%clb.Label.HorizontalAlignment='right';
%[c2,h2]=contour(xq,yq,z2,'showtext','on','linecolor',clrs.red,'linestyle','-.');
%[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.red,'linestyle','-.');
%clabel(c2,h2,'Color',clrs.red)



figure
hold on
[c1,h1]=contourf(xq,yq,z1);%,0.1:0.2:1.5,'showtext','on');
colormap(flip(gray))
%clb=colorbar;
%clb.Label.String='Probability Density';
%clb.Label.FontSize=16;


% PLOT EACH MEGADROUGHT PDF INDIVIDUALLY
X2=cat(2,mdagr_na(:),mdagr_sa(:));
[f2,x2] = ksdensity(X2);
z2 = griddata(x2(:,1),x2(:,2),f2,xq,yq);

%[c2,h2]=contour(xq,yq,z2,0.1:0.2:1.5,'showtext','on','linecolor',clrs.red,'linestyle','-','linewidth',1.25);
[c2,h2]=contour(xq,yq,z2,'linecolor',clrs.blue,'linestyle','-','linewidth',1.25);
clabel(c2,h2,'Color',clrs.blue)


scatter(mdagr_na(:,1),mdagr_sa(:,1),'jitter','on','jitterAmount',0.25,'MarkerEdgeColor',clrs.navy)

set(gca,'fontsize',16)
xlabel('NASW drought prob.','fontsize',18)
ylabel('Pampas drought prob.','fontsize',18)
box on

print(['/home/nsteiger/Dropbox/sa-phyda/sa_na_droughtprob_null.png'],'-dpng','-r300');

return

% Find the max and plot it
[~,m1]=max(max(z2,[],1));[~,m2]=max(max(z2,[],2));
plot(xq(m2,m1),yq(m2,m1),'*','linewidth',2,'color',clrs.red,'markersize',10,'markerfacecolor',clrs.red)

line([xl(1),xl(2)],[0,0],'color','k','linestyle','--','linewidth',1.5)
line([0,0],[yl(1),yl(2)],'color','k','linestyle','--','linewidth',1.5)
%set(gca,'tickdir','both')
%xlim([-3.75 3.75])
xlim([-2.5 2.5])
ylim([-0.7 0.7])
set(gca,'fontsize',16)
xlabel('Nino3.4 (annual temp. anomaly ^{\circ}C)','fontsize',18)
ylabel('AMO (temp. anomaly ^{\circ}C)','fontsize',18)
box on







