load('Figure3-5.mat');

% colored figure of wave normal: percentage per mlat & ray tracing
% comparison
f2=figure(2);
set(gcf,'Position',[100,150,1800,800])

valued=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5);
wna_bin=0:10:180;
mlat_bin=-20:1:10;
mlat_wna=nan(31,19);
for i=1:30
    for j=1:18
        if length(find(mlat(valued)>mlat_bin(i) & mlat(valued)<mlat_bin(i+1)))>100
            mlat_wna(i,j)=length(find(mlat(valued)>mlat_bin(i) & mlat(valued)<mlat_bin(i+1) & ...
                Wave_3B3E(valued)>wna_bin(j) & Wave_3B3E(valued)<wna_bin(j+1))) / ...
                length(find(mlat(valued)>mlat_bin(i) & mlat(valued)<mlat_bin(i+1)));
        end
    end
end
s21=subplot(1,2,1);
pcolor(mlat_bin,wna_bin,mlat_wna');
hold on
text(-25,180,'(a)','FontSize',25)
shading flat
colormap(jet)
xlim([-20 9])
xlabel('MLAT (°)')
ylabel('Polar angle of k vector (°)')
%title('Alpha signals distribution')
set(gca,'XMinorTick','on','YMinorTick','on')
set(gca,'TickDir','out')
ax1=gca;
ax1.YRuler.MinorTickValues=0:10:180;
c21=colorbar('eastoutside');
c21.Label.String='Normalized signal number distribution';
set(gca,'FontSize',20)


% plot the ray-tracing results with the quantiles
mlat_number=[];
quan_wna_mlat=[];
mlat_min=min(mlat(valued));
mlat_max=max(mlat(valued));
for i=floor(mlat_min):1:floor(mlat_max)
    mlat_index=find(mlat(valued)>i & mlat(valued)<i+1);
    if length(mlat_index)>100
        mlat_number=[mlat_number;i];
        quan_wna_mlat=[quan_wna_mlat;quantile(Wave_3B3E(valued(mlat_index)),3)];
    end
end

load('Figure4_rt.mat');

s22=subplot(1,2,2);
for i=1:3
    for j=1:36
        hold on
        index=find(output{i}.hout.rtd{j}.lat>-20 & output{i}.hout.rtd{j}.lat<10);
        index1=find(output{i}.hout1.rtd{j}.lat>-20 & output{i}.hout1.rtd{j}.lat<10);
        index2=find(output{i}.hout2.rtd{j}.lat>-20 & output{i}.hout2.rtd{j}.lat<10);
        
        if i==1
            p1=plot(output{i}.hout.rtd{j}.lat(index),output{i}.hout.rtd{j}.psid(index),'r');
        elseif i==2
            p2=plot(output{i}.hout.rtd{j}.lat(index),output{i}.hout.rtd{j}.psid(index),'g');
        else
            p3=plot(output{i}.hout.rtd{j}.lat(index),output{i}.hout.rtd{j}.psid(index),'b');
        end       
    end 
end

p4=plot(mlat_number+0.5,quan_wna_mlat(:,1),'LineStyle','--','Linewidth',1.5,'Color','k');
p5=plot(mlat_number+0.5,quan_wna_mlat(:,2),'LineStyle','-','Linewidth',1.5,'Color','k');
p6=plot(mlat_number+0.5,quan_wna_mlat(:,3),'LineStyle',':','Linewidth',1.5,'Color','k');
text(-26,150,'(b)','FontSize',25)
text(-18,125,'75%','Color','k','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
text(-18,112,'50%','Color','k','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
text(-18,94,'25%','Color','k','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
text(12,146,'NOV','Color','g','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
text(12,139,'KHA','Color','r','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
text(12,141,'KRA','Color','b','FontSize',18,'HorizontalAlignment','center','FontWeight','bold')
xlim([-20 14])
ylim([90 150])
xlabel('MLAT (°)')
ylabel('Polar angle of k vector (°)')

set(gca,'FontSize',20)
s21.Position=[0.08 0.11 0.38 0.815];
s22.Position=[0.60 0.11 0.38 0.815];
ah1=axes('position',get(gca,'position'),'visible','off');

