load('Figure3-5.mat');

% duct signals with L and Kp, south of MLAT=-5deg
valued=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5 & Wave_3B3E>90 & Wave_P>90 & mlat<-5);
in20=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5 & Wave_3B3E>160 & Wave_P>90 & mlat<-5);
in30=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5 & Wave_3B3E>150 & Wave_P>90 & mlat<-5);
in40=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5 & Wave_3B3E>140 & Wave_P>90 & mlat<-5);
in35=find((log10(B_str)>-3.75 | log10(E_str)>-2.75) & Plan_3B3E>0.75 & ...
    log10(B_str)<5 & log10(E_str)<5 & Wave_3B3E>145 & Wave_P>90 & mlat<-5);


f3=figure(3);
set(gcf,'Position',[250,300,1400,1100])

lshell_bin=1.1:0.2:3.7;
in20_l_num=nan(1,length(lshell_bin));
in30_l_num=nan(1,length(lshell_bin));
in40_l_num=nan(1,length(lshell_bin));
in35_l_num=nan(1,length(lshell_bin));
total_l_num=nan(1,length(lshell_bin));
for i=1:length(lshell_bin)
    in20_l_num(i)=length(find(lshell(in20)>lshell_bin(i) & lshell(in20)<lshell_bin(i)+0.2));
    in30_l_num(i)=length(find(lshell(in30)>lshell_bin(i) & lshell(in30)<lshell_bin(i)+0.2));
    in40_l_num(i)=length(find(lshell(in40)>lshell_bin(i) & lshell(in40)<lshell_bin(i)+0.2));
    in35_l_num(i)=length(find(lshell(in35)>lshell_bin(i) & lshell(in35)<lshell_bin(i)+0.2));
    total_l_num(i)=length(find(lshell(valued)>lshell_bin(i) & lshell(valued)<lshell_bin(i)+0.2));
    if total_l_num(i)<100
        in20_l_num(i)=nan;
        in30_l_num(i)=nan;
        in40_l_num(i)=nan;
        in35_l_num(i)=nan;
    end
end

s31=subplot(2,2,1);
yyaxis left
hold on
p311=plot(lshell_bin+0.1,in35_l_num./total_l_num,'b','Linewidth',1.5);
ymax=ylim;
text(1.0,ymax(2),'(a)','FontSize',20)
ylabel('Ducted pulses occurrence rate')
hold off
yyaxis right
p312=plot(lshell_bin+0.1,total_l_num,'r','LineStyle','-','Linewidth',1.5);
xlim([1.4 3.4])
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(2).Color = 'r';
ylabel('Total signal number (MLAT < -5°)')
xlabel('L')
xticks(1.2:0.2:3.6)
set(gca,'FontSize',15)

Kp_max=floor((Kp_max+4)/10);
Kp_max_bin=0:8;
in20_Kp_max_num=nan(1,length(Kp_max_bin));
in30_Kp_max_num=nan(1,length(Kp_max_bin));
in40_Kp_max_num=nan(1,length(Kp_max_bin));
in35_Kp_max_num=nan(1,length(Kp_max_bin));
total_Kp_max_num=nan(1,length(Kp_max_bin));
for i=1:length(Kp_max_bin)
    in20_Kp_max_num(i)=length(find(Kp_max(in20)==Kp_max_bin(i)));
    in30_Kp_max_num(i)=length(find(Kp_max(in30)==Kp_max_bin(i)));
    in40_Kp_max_num(i)=length(find(Kp_max(in40)==Kp_max_bin(i)));
    in35_Kp_max_num(i)=length(find(Kp_max(in35)==Kp_max_bin(i)));
    total_Kp_max_num(i)=length(find(Kp_max(valued)==Kp_max_bin(i)));
    if total_Kp_max_num(i)<100
        in20_Kp_max_num(i)=nan;
        in30_Kp_max_num(i)=nan;
        in40_Kp_max_num(i)=nan;
        in35_Kp_max_num(i)=nan;
    end
end

s32=subplot(2,2,2);
yyaxis left
hold on
p321=plot(Kp_max_bin,in35_Kp_max_num./total_Kp_max_num,'b','Linewidth',1.5);
ymax=ylim;
text(-1.3,ymax(2),'(b)','FontSize',20)
ylabel('Ducted pulses occurrence rate')
hold off
yyaxis right
p322=plot(Kp_max_bin,total_Kp_max_num,'r','LineStyle','-','Linewidth',1.5);
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(2).Color = 'r';
ylabel('Total signal number (MLAT < -5°)')
xlim([0 7])
xlabel('Maximum Kp in previous 48h')
set(gca,'FontSize',15)

% get the ducted and nonducted power using signals from the south of
% MLAT=5deg
wna_bond=35;
ducted=find(mlat(valued)<-5 & Wave_3B3E(valued)>180-wna_bond & Wave_P(valued)>90 & Wave_3B3E(valued)>90);
nonducted=find(mlat(valued)<-5 & Wave_3B3E(valued)<180-wna_bond & Wave_P(valued)>90 & Wave_3B3E(valued)>90);
mean_all_B_ducted=mean(B_str_mean(valued(ducted)));
mean_all_B_nonducted=mean(B_str_mean(valued(nonducted)));
ratio_all_B_ducted=(mean_all_B_ducted^2)*length(ducted)/((mean_all_B_ducted^2)*length(ducted) ...
    +(mean_all_B_nonducted^2)*length(nonducted));
ratio_all_B_nonducted=(mean_all_B_nonducted^2)*length(nonducted)/((mean_all_B_ducted^2)*length(ducted) ...
    +(mean_all_B_nonducted^2)*length(nonducted));

% get the case and power ratio for lshells
lshell_bin=1.1:0.2:3.7;
ratio_L_ducted_power=nan(length(lshell_bin),1);
ratio_L_ducted_case=nan(length(lshell_bin),1);
L_total_case=nan(length(lshell_bin),1);
L_total_power=nan(length(lshell_bin),1);
for i=1:length(lshell_bin)
    ducted_list=find(lshell(valued(ducted))>lshell_bin(i) & lshell(valued(ducted))<lshell_bin(i)+0.2);
    nonducted_list=find(lshell(valued(nonducted))>lshell_bin(i) & lshell(valued(nonducted))<lshell_bin(i)+0.2);
    L_total_case(i)=length(nonducted_list)+length(ducted_list);
    ratio_L_ducted_case(i)=length(ducted_list)/(length(nonducted_list)+length(ducted_list));
    ratio_L_ducted_power(i)=mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)/ ...
        (mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)+mean(B_str_mean(valued(nonducted(nonducted_list))).^2)*length(nonducted_list));
    L_total_power(i)=mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)*0.4+mean(B_str_mean(valued(nonducted(nonducted_list))).^2)*length(nonducted_list)*0.4;
    if L_total_case(i)<100
        ratio_L_ducted_case(i)=nan;
        ratio_L_ducted_power(i)=nan;
    end
end

s33=subplot(2,2,3);
yyaxis left
hold on
p331=plot(lshell_bin+0.1,ratio_L_ducted_power,'b','LineStyle','-','Linewidth',1.5);
text(1.0,0.8,'(c)','FontSize',20)
ylabel('Ducted energy density proportion')
yyaxis right
p332=plot(lshell_bin+0.1,L_total_power,'r','LineStyle','-','Linewidth',1.5);
hold off
ylabel('Total observed energy density (nT^2s) \newline (MLAT < -5°)')
xlabel('L')
xlim([1.4 3.4])
xticks(1.2:0.2:3.6)
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(2).Color = 'r';
box on
set(gca,'FontSize',15)

Kp_max_bin=0:8;
ratio_Kp_max_ducted_power=nan(length(Kp_max_bin),1);
ratio_Kp_max_ducted_case=nan(length(Kp_max_bin),1);
Kp_max_total_case=nan(length(Kp_max_bin),1);
Kp_max_total_power=nan(length(Kp_max_bin),1);
for i=1:length(Kp_max_bin)
    ducted_list=find(Kp_max(valued(ducted))==Kp_max_bin(i));
    nonducted_list=find(Kp_max(valued(nonducted))==Kp_max_bin(i));
    Kp_max_total_case(i)=length(nonducted_list)+length(ducted_list);
    ratio_Kp_max_ducted_case(i)=length(ducted_list)/(length(nonducted_list)+length(ducted_list));
    ratio_Kp_max_ducted_power(i)=mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)/ ...
        (mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)+mean(B_str_mean(valued(nonducted(nonducted_list))).^2)*length(nonducted_list));
    Kp_max_total_power(i)=mean(B_str_mean(valued(ducted(ducted_list))).^2)*length(ducted_list)*0.4+mean(B_str_mean(valued(nonducted(nonducted_list))).^2)*length(nonducted_list)*0.4;
    if Kp_max_total_case(i)<100
        ratio_Kp_max_ducted_case(i)=nan;
        ratio_Kp_max_ducted_power(i)=nan;
    end
end

s34=subplot(2,2,4);
yyaxis left
hold on
p341=plot(Kp_max_bin,ratio_Kp_max_ducted_power,'b','LineStyle','-','Linewidth',1.5);
text(-1.3,0.6,'(d)','FontSize',20)
ylabel('Ducted energy density proportion')
yyaxis right
plot(Kp_max_bin,Kp_max_total_power,'r','LineStyle','-','Linewidth',1.5);
hold off
ylabel('Total observed energy density (nT^2s) \newline (MLAT < -5°)')
xlim([0 7])
ax = gca;
ax.YAxis(1).Color = 'b';
ax.YAxis(2).Color = 'r';
box on
xlabel('Maximum Kp in previous 48h')
set(gca,'FontSize',15)

s31.Position=[0.08 0.55 0.35 0.42];
s32.Position=[0.57 0.55 0.35 0.42];
s33.Position=[0.08 0.06 0.35 0.42];
s34.Position=[0.57 0.06 0.35 0.42];
leg1.Position=[0.2900 0.8716 0.0779 0.0880];
leg2.Position=[0.8100 0.8716 0.0779 0.0880];