%    This Code is for L&OL article "Transport and environmental impact of ash 
%induced by the Hunga Tonga-Hunga Ha'apai volcanic eruption"
%    Authors:Yuntao Wang, Huaguo Zhang, Chang Liu, Mei Zheng, Huan-Huan Chen, 
%Wenting Cao, Shuangling Chen, Chen Wang, Xuehua Fan, Zhihua Mao, Fei Chai
%--------------------------------------------------------------------------
% Fig.2----CALIPSO Backscatter
clear;
load("BackScatter.mat")
xlim_Reg = [70 200];ylim_Reg = [-28 -8];
xtam=diff(xlim_Reg)*cos(mean(ylim_Reg)*pi/180);
ytam=diff(ylim_Reg);
c0 = [0 64 183;1 65 184;11 147 255;12 148 255;13 149 255;
      14 150 255;15 151 255;13 251 184;1 144 147;2 181 104];
c1 = [254 252 19;255 253 20;255 218 3;253 184 0;251 148 16;...
      255 106 1;234 51 35;252 71 102;254 112 145;255 150 183;...
      88 86 87;119 119 119;148 148 148;171 171 171;];
for i = 1:14
    c11(i*5-4:i*5,1:3) = repmat(c1(i,:),[5,1]);
end
c2 = [193 193 193];
c22(1:20,1:3) = repmat(c2(1,:),[20,1]);
c3 = [211 211 211;231 231 231;240 240 240;244 244 244;...
      245 245 245;248 248 248;250 250 250;252 252 252;255 255 255];
for i = 1:9
    c33(i*100-99:i*100,1:3) = repmat(c3(i,:),[100,1]);
end
color = [c0;c11;c22;c33]/256;
yt = 0:5000:40000;
yl = ["0","5","10","15","20","25","30","35","40"];

ax1 = subplot('position',[0.05,0.55,0.4,0.4]);
xxx = cell2mat(lon_BS(1));
xxx = flipud(xxx);
xt = 150:floor((length(xxx)-150-150)/4):length(xxx)-150;
xi = xxx(xt);
ll = cell2mat(lat_BS(1));
ll = flipud(ll);
lt = 150:floor((length(ll)-150-150)/4):length(ll)-150;
xlon = {'177.44°W','176.90°W','176.35°W','175.79°W','175.21°W'};
xlat = {'15.76°S','18.14°S','20.51°S','22.88°S','25.42°S'};
temp = flipud(cell2mat(BackScatter(1)));
t = wiener2(temp,[10 10]);
pcolor(xxx,Height,t'),shading flat
colormap(color);
caxis([1*10^-4 0.1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
xticks(xi);
xticklabels(xlon);
text(xi-0.17,[-5000,-5000,-5000,-5000,-5000],xlat,'fontsize',11.5,'FontWeight','bold');
yticks(yt);
yticklabels(yl);
text(xxx(100),0,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(xxx(100),0,'A','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(280),0,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(xxx(280),0,'15 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(1)-0.2,-8000,'Altitude (km)','fontsize',11.5,'FontWeight','bold','Rotation',90);
ax5 = axes('Position',[0.2,0.9,0.25 0.05],'Box','on');
axes(ax5)
contour(lon_coast,lat_coast,Coast,[1,1],'-k')
hold on
plot(xxx,ll,'or','MarkerSize',2,'MarkerFaceColor','r');
xlim([90 200]);ylim([-28,-4]);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);

ax2 = subplot('position',[0.48,0.55,0.4,0.4]);
xxx = cell2mat(lon_BS(2));
xxx = flipud(xxx);
xt = 150:floor((length(xxx)-150-150)/4):length(xxx)-150;
xi = xxx(xt);
ll = cell2mat(lat_BS(2));
ll = flipud(ll);
lt = 150:floor((length(ll)-150-150)/4):length(ll)-150;
xlon = {'159.22°E','159.87°E','160.49°E','161.11°E','161.70°E'};
xlat = {'25.96°S','23.34°S','20.71°S','18.08°S','15.45°S'};
temp = flipud(cell2mat(BackScatter(2)));
pcolor(xxx,Height,temp'),shading flat
colormap(color);
caxis([1*10^-4 0.1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
xticks(xi);
xticklabels(xlon);
text(xi-0.155,[-5000,-5000,-5000,-5000,-5000],xlat,'fontsize',11.5,'FontWeight','bold');
% yticks(yt);
% yticklabels(yl);
text(xxx(100),0,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(xxx(100),0,'B','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(280),0,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(xxx(280),0,'16 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(end)+0.6,15000,'Total Attenuated Backscatter 532 (km^{-1}sr^{-1})','fontsize',11.5,'FontWeight','bold','Rotation',-90);

ax6 = axes('Position',[0.63,0.9,0.25 0.05],'Box','on');
axes(ax6)
contour(lon_coast,lat_coast,Coast,[1,1],'-k')
hold on
plot(xxx,ll,'or','MarkerSize',2,'MarkerFaceColor','r');
xlim([90 200]);ylim([-28,-4]);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);

ax3 = subplot('position',[0.05,0.1,0.4,0.4]);
xxx = cell2mat(lon_BS(3));
xxx = flipud(xxx);
xt = 200:floor((length(xxx)-200-200)/4):length(xxx)-200;
xi = xxx(xt);
ll = cell2mat(lat_BS(3));
ll = flipud(ll);
lt = 200:floor((length(ll)-200-200)/4):length(ll)-200;
xlon = {'138.48°E','139.17°E','139.88°E','140.60°E','141.35°E'};
xlat = {'12.72°S','15.81°S','18.89°S','21.97°S','25.05°S'};
temp = flipud(cell2mat(BackScatter(3)));
t = wiener2(temp,[5 5]);
pcolor(xxx,Height,t'),shading flat
colormap(color);
caxis([1*10^-4 0.1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
xticks(xi);
xticklabels(xlon);
text(xi-0.2,[-5000,-5000,-5000,-5000,-5000],xlat,'fontsize',11.5,'FontWeight','bold');
yticks(yt);
yticklabels(yl);
text(xxx(100),0,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(xxx(100),0,'C','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(340),0,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(xxx(340),0,'17 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
ax7 = axes('Position',[0.2,0.45,0.25 0.05],'Box','on');
axes(ax7)
contour(lon_coast,lat_coast,Coast,[1,1],'-k')
hold on
plot(xxx,ll,'or','MarkerSize',2,'MarkerFaceColor','r');
xlim([90 200]);ylim([-28,-4]);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);

ax4 = subplot('position',[0.48,0.1,0.4,0.4]);
xxx = cell2mat(lon_BS(4));
xxx = flipud(xxx);
xt = 300:floor((length(xxx)-300-300)/4):length(xxx)-300;
xi = xxx(xt);
ll = cell2mat(lat_BS(4));
ll = flipud(ll);
lt = 300:floor((length(ll)-300-300)/4):length(ll)-300;
xlon = {'93.93°E','94.82°E','95.73°E','96.67°E','97.64°E'};
xlat = {'7.69°S','11.64°S','15.59°S','19.54°S','23.48°S'};
temp = flipud(cell2mat(BackScatter(4)));
t = wiener2(temp,[3 3]);
pcolor(xxx,Height,t'),shading flat
colormap(color);
caxis([1*10^-4 0.1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
xticks(xi);
xticklabels(xlon);
text(xi-0.255,[-5000,-5000,-5000,-5000,-5000],xlat,'fontsize',11.5,'FontWeight','bold');

% yticks(yt);
% yticklabels(yl);
text(xxx(100),0,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(xxx(100),0,'D','Color','k','fontsize',11.5,'FontWeight','bold');
text(xxx(410),0,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(xxx(410),0,'19 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
set(gcf,'Position',[100 100 900 900])
ax8 = axes('Position',[0.63,0.45,0.25 0.05],'Box','on');
axes(ax8)
contour(lon_coast,lat_coast,Coast,[1,1],'-k')
hold on
plot(xxx,ll,'or','MarkerSize',2,'MarkerFaceColor','r');
xlim([90 200]);ylim([-28,-4]);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);

colormap(ax8,[c0;c1;c2;c3]/256);
% caxis([1,34]);
cb3 = colorbar(ax8,'eastoutside');
set(cb3,'position',[0.9,0.125,0.01,0.8]);
ct3 = (0:(2-0)/34:2)+2/(2*34);
ct3label = ["1×10^{-4}","2","3","4","5","6","7","8","9",...
            "1×10^{-3}","1.5","2","2.5","3","3.5","4","4.5","5","5.5","6","6.5","7","7.5","8",...
            "1×10^{-2}","2","3","4","5","6","7","8","9",...
            "1×10^{-1}"];
set(cb3,'Ticks',ct3(1:end-1),'TickLabels',ct3label);
