clear all
load('2011data.mat')
BCmp=load('2011_bcmap2')
DSTmp=load('2011_dustmap2')
  
newBCmp=BCmp;
newDSTmp=DSTmp;
newEALd=EALd;
newEALb=EALb;
newEALd_TP=EALd_TP;
newEALb_TP=EALb_TP;

newDA_HK=DA_HK;
newDA_WH=DA_WH;
newDA_KK=DA_KK;
newDA_TP=DA_TP;

clear EALd DA_HK DA_WH DA_KK DA_TP EALb EALb_TP EALd_TP BCmp DSTmp EALmap

load('2012data.mat')
BCmp=load('2012_bcmap2')
DSTmp=load('2012_dustmap2')
  
newBCmp=(BCmp.EALmap+newBCmp.EALmap)/2;
newDSTmp=(DSTmp.EALmap+newDSTmp.EALmap)/2;
newEALd=[newEALd EALd];
newEALb=[newEALb EALb];
newEALd_TP=[newEALd_TP EALd_TP];
newEALb_TP=[newEALb_TP EALb_TP];

newDA_HK=[newDA_HK DA_HK];
newDA_WH=[newDA_WH DA_WH];
newDA_KK=[newDA_KK DA_KK];
newDA_TP=[newDA_TP DA_TP];

clear EALd DA_HK DA_WH DA_KK DA_TP EALb EALb_TP EALd_TP BCmp DSTmp EALmap

load('2013data.mat')
BCmp=load('2013_bcmap2')
DSTmp=load('2013_dustmap2')
  
newBCmp=(BCmp.EALmap+newBCmp)/2;
newDSTmp=(DSTmp.EALmap+newDSTmp)/2;
newEALd=[newEALd EALd];
newEALb=[newEALb EALb];
newEALd_TP=[newEALd_TP EALd_TP];
newEALb_TP=[newEALb_TP EALb_TP];

newDA_HK=[newDA_HK DA_HK];
newDA_WH=[newDA_WH DA_WH];
newDA_KK=[newDA_KK DA_KK];
newDA_TP=[newDA_TP DA_TP];
clear EALd DA_HK DA_WH DA_KK DA_TP EALb EALb_TP EALd_TP BCmp DSTmp EALmap

load('2014data.mat')
% BCmp=load('2014_bcmap2')
% DSTmp=load('2014_dustmap2')
%   
% newBCmp=nanmean(BCmp.EALmap+newBCmp.EALmap);
% newDSTmp=nanmean(DSTmp.EALmap+newDSTmp.EALmap);
newEALd=[newEALd EALd];
newEALb=[newEALb EALb];
newEALd_TP=[newEALd_TP EALd_TP];
newEALb_TP=[newEALb_TP EALb_TP];
newDA_HK=[newDA_HK DA_HK];
newDA_WH=[newDA_WH DA_WH];
newDA_KK=[newDA_KK DA_KK];
newDA_TP=[newDA_TP DA_TP];
clear EALd DA_HK DA_WH DA_KK

load('2015data.mat')
BCmp=load('2015_bcmap2')
DSTmp=load('2015_dustmap2')
  
newBCmp=(BCmp.EALmap+newBCmp)/2;
newDSTmp=(DSTmp.EALmap+newDSTmp)/2;
newEALd=[newEALd EALd];
newEALb=[newEALb EALb];
newEALd_TP=[newEALd_TP EALd_TP];
newEALb_TP=[newEALb_TP EALb_TP];
newDA_HK=[newDA_HK DA_HK];
newDA_WH=[newDA_WH DA_WH];
newDA_KK=[newDA_KK DA_KK];
newDA_TP=[newDA_TP DA_TP];
clear EALd DA_HK DA_WH DA_KK DA_TP EALb EALb_TP EALd_TP BCmp DSTmp EALmap

load('2016data.mat')
BCmp=load('2016_bcmap2')
DSTmp=load('2016_dustmap2')
  
newBCmp=(BCmp.EALmap+newBCmp)/2;
newDSTmp=(DSTmp.EALmap+newDSTmp)/2;
newEALd=[newEALd EALd];
newEALb=[newEALb EALb];
newEALd_TP=[newEALd_TP EALd_TP];
newEALb_TP=[newEALb_TP EALb_TP];
newDA_HK=[newDA_HK DA_HK];
newDA_WH=[newDA_WH DA_WH];
newDA_KK=[newDA_KK DA_KK];
newDA_TP=[newDA_TP DA_TP];
clear EALd DA_HK DA_WH DA_KK DA_TP EALb EALb_TP EALd_TP BCmp DSTmp EALmap

EALd=newEALd*1000;
EALb=newEALb;
EALb_TP=newEALb_TP;
EALd_TP=newEALd_TP;

DA_HK=newDA_HK;
DA_WH=newDA_WH;
DA_KK=newDA_KK;
DA_TP=newDA_TP;

%EALd=EALb

clear d *bin
%id=EALd(:)>.1
d=prctile(EALd(:),[10:3:100]);
for i=1:length(d)-1
    id=EALd>d(i) & EALd<d(i+1)
     EALdbin(i)=nanmean(EALd(id));
     DA_TPbin(i)=nanmean(DA_TP(id));
     DA_WHbin(i)=nanmean(DA_WH(id));
    
     DA_HKbin(i)=nanmean(DA_HK(id));
     DA_KKbin(i)=nanmean(DA_KK(id));
     
    % DA_CHbin(i)=nanmean(DA_CH(id));
     stdDA_TPbin(i)=nanstd(DA_TP(id));
     stdDA_WHbin(i)=nanstd(DA_WH(id));
   
     stdDA_HKbin(i)=nanstd(DA_HK(id));
     stdDA_KKbin(i)=nanstd(DA_KK(id));
    % stdDA_CHbin(i)=nanstd(DA_CH(id));
end
%C=varycolor(4);

figure(3)
pos3 = [0.1 0.1 0.37 0.4];
pos1 = [0.1 0.57 0.39 0.38];

pos4 = [0.57 0.1 0.37 0.4];
pos2 = [0.57 0.57 0.37 0.4];

%load('mycolormap.mat')
%colormap(cmap)


subplot('position',pos4)

scatter(log10(EALd(:)),(DA_WH(:)),'MarkerFaceColor','r','MarkerEdgeColor','r',...
    'MarkerFaceAlpha',.1,'MarkerEdgeAlpha',.1)
hold on
errorbar(log10(EALdbin),(DA_WHbin),(stdDA_WHbin),'rs','Linewidth',1,'MarkerSize',10,'MarkerFaceColor','red')
hold on
[p2 S]= polyfit(log10(EALdbin),(DA_WHbin),2)
 [y delta] = polyval(p2,log10(EALdbin),S);
hold on
plot((log10(EALdbin)),y,'--r','Linewidth',3)
box on
[px2 rmse] = corrcoef(log10(EALdbin),(DA_WHbin)); 
caption = sprintf('Y = %.2f * X^2 +%.2f * X + %.2f', p2(1), p2(2),p2(3));
yut=1-(S.normr/norm(y - mean(y)))^2
text(.1,.95, caption, 'Fontsize',15,'Units','normalized','color','r');
%yut=1-(S.normr/norm(y - mean(y)))^2
 text(.5,.12,[ 'WH: R^2 =', num2str(yut,'%.2f') ],'Fontsize',15,'Units','normalized','color','r')

  %text(.5,.12,[ 'WH: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','r')
xlim([1.7 3.5])
ylim([2 8])
  xlabel('log [EAL_{dust} (\mug/m^3)]')
    ylabel('\Delta Albedo (%)')
  set(gca,'fontsize',15)
 %  text(.2,.95,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','r') 




subplot('position',pos3)
scatter(log10(EALd(:)),(DA_HK(:)),'MarkerFaceColor','b','MarkerEdgeColor','b',...
    'MarkerFaceAlpha',.1,'MarkerEdgeAlpha',.1)
hold on
errorbar(log10(EALdbin),(DA_HKbin),(stdDA_HKbin),'bs','Linewidth',1,'MarkerSize',10,'MarkerFaceColor','b')
hold on
[p2 S]= polyfit(log10(EALdbin),(DA_HKbin),2)
 [y delta] = polyval(p2,log10(EALdbin),S);
hold on
plot(log10(EALdbin),y,'--b','Linewidth',3)
box on
[px2 rmse] = corrcoef(log10(EALdbin),(DA_HKbin)); 
caption = sprintf('Y = %.2f * X^2 + %.2f * X + %.2f', p2(1), p2(2),p2(3));
yut=1-(S.normr/norm(y - mean(y)))^2
 text(.5,.12,[ 'HK: R^2 =', num2str(yut,'%.2f') ],'Fontsize',15,'Units','normalized','color','b')

text(.1,.95, caption, 'Fontsize',15,'Units','normalized','color','b');

% text(.5,.12,[ 'HK: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','b')
xlim([1.7 3.5])
  xlabel('log [EAL_{dust} (\mug/m^3)]')
    ylabel('\Delta Albedo (%)')
  set(gca,'fontsize',15)
  % text(.2,.9,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','b') 
ylim([2 10])



subplot('position',pos2)
scatter(log10(EALd(:)),(DA_KK(:)),'MarkerFaceColor','k','MarkerEdgeColor','k',...
    'MarkerFaceAlpha',.1,'MarkerEdgeAlpha',.1)
hold on
errorbar(log10(EALdbin),(DA_KKbin),(stdDA_KKbin),'ks','Linewidth',1,'MarkerSize',10,'MarkerFaceColor','k')
hold on
[p2 S]= polyfit(log10(EALdbin),(DA_KKbin),2)
 [y delta] = polyval(p2,log10(EALdbin),S);
hold on
plot(log10(EALdbin),y,'--k','Linewidth',3)
box on
[px2 rmse] = corrcoef(log10(EALdbin),(DA_KKbin)); 
caption = sprintf('Y = %.2f * X^2 +%.2f * X + %.2f', p2(1), p2(2),p2(3));
yut=1-(S.normr/norm(y - mean(y)))^2
text(.1,.95, caption, 'Fontsize',15,'Units','normalized','color','k');
 text(.5,.12,[ 'KK: R^2 = ',num2str(yut,'%.2f') ],'Fontsize',15,'Units','normalized','color','k')
xlim([1.7 3.5])
 % xlabel('log [EAL_{dust} (\mug/m^3)]')
    ylabel('\Delta Albedo (%)')
  set(gca,'fontsize',15)
 ylim([2 7])
 
 
 load('Dustmap.mat')
 newdustmap=dustmap'

varid= netcdf.inqVarID(ncid,'lat'); 
LAT= netcdf.getVar(ncid,varid);
XLAT=(mean(LAT(:,:,1),3)) ;




varid= netcdf.inqVarID(ncid,'lon'); 
LONG= netcdf.getVar(ncid,varid);
XLONG=(mean(LONG(:,:,1),3)) ;


varid= netcdf.inqVarID(ncid,'lev'); 
lev= netcdf.getVar(ncid,varid);
level=(mean(lev(:,:,1),3)) ;

varid= netcdf.inqVarID(ncid,'time'); 
lev= netcdf.getVar(ncid,varid);
time=lev ;
load('/Volumes/Untitled/HMA/paper2_analysis/mycolormap.mat')

 subplot('position',pos1)
 jet=colormap('parula');
 clear inferno
inferno=inferno(100);

%magma=magma();
%inferno=inferno();
%plasma=plasma();
%viridis=viridis();

%colormap(cmap)
colormap('inferno')
% axes1 = axes('Parent',figure1)
% hold(axes1,'all');
[c,h]=contourf(LONG, LAT, newDSTmp',100)
set(h,'edgecolor','none')
%colorbar('northoutside')
c=colorbar('Position',[0.2 0.96 0.3 0.01],'Location','north');
c.Label.String = 'EAL_{dust} (mg/m^3)';
c.Label.FontSize = 15;

ylabel('Latitude (^o N)')
xlabel('Longitude (^o E)')
set(gca,'fontsize',15)
hold on

S=shaperead('/Volumes/Untitled/HMA/paper1_analysis/00_rgi60_regions/00_rgi60_O2Regions.shp');
plot([S.X],[S.Y],'LineWidth',2,'color','k','Linestyle','-')
xlim([67 90])
ylim([22 38])
%caxis([.001 1])
caxis([.3 .8])
%rectangle('Position',[67 25 10 6],'LineWidth',2)
  % text(.2,.85,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','m') 
  
  annotation('textbox',[0.01 .89 0.25 .05],'String',{'a)'},'LineStyle','none','FontSize',20);
annotation('textbox',[0.52 .89 0.25 .05],'String',{'b)'},'LineStyle','none','FontSize',20);
annotation('textbox',[0.01 .43 0.25 .05],'String',{'c)'},'LineStyle','none','FontSize',20);
annotation('textbox',[0.52 .43 0.25 .05],'String',{'d)'},'LineStyle','none','FontSize',20);
  
 
% text(83,32.5,'Tibetan Plateau (TP)','Rotation',1,'fontsize',13)
 text(70.25,36,'HK','color','b','Rotation',5,'fontsize',23)
 text(74.5,37,'KK','color','k','Rotation',-35,'fontsize',23)
 text(75,34.7,'WH','color','r','Rotation',-35,'fontsize',23)
 %text(78,30.7,'C Himalaya (CH)','Rotation',-30,'fontsize',13)
% text(86,28.5,'E Him (EH)','Rotation',-25,'fontsize',13)

%%
