clear all;clc;close all;
%% ---------------------------------------------------------
load('Figure2.mat');

f1=figure(1);
set(f1,'Position',[200 150 1100 700])

s1=subplot(2,2,1)
[nonduct.rtd.gtime,non_rep]=unique(nonduct.rtd.gtime);
nonduct.rtd.x=nonduct.rtd.x(non_rep);
nonduct.rtd.y=nonduct.rtd.y(non_rep);
nonduct.rtd.z=nonduct.rtd.z(non_rep);
nonduct.rtd.rrr=nonduct.rtd.rrr(non_rep);
nonduct.rtd.rkxbig=nonduct.rtd.rkxbig(non_rep);
nonduct.rtd.rkybig=nonduct.rtd.rkybig(non_rep);
nonduct.rtd.rkzbig=nonduct.rtd.rkzbig(non_rep);

xm=sqrt(nonduct.n3d.x.^2+nonduct.n3d.y.^2);
xy=sqrt(nonduct.rtd.x.^2+nonduct.rtd.y.^2);
pcolor(xm,nonduct.n3d.z,log10(nonduct.n3d.dent1))
shading flat
hold on
% contour(xm,hout.n3d.z,hout.n3d.flhr,[0 0]+freq0*1d-3,'m')
plot_Earth('2M','l',[1.5:0.5:3])
xlabel('R_x_y (R_E)')
ylabel('Z (R_E)')
axis equal
c1=colorbar('eastoutside');
colormap(jet)
c1.Label.String='log(density(m^-^3))';
caxis([8.5,10])
xlim([0,3])
ylim([-1.5,1.5])
hold on
plot(xy,nonduct.rtd.z,'color','m','linewidth',2); % ray path

phi_xy=atan2(nonduct.rtd.x,nonduct.rtd.y);
phi_k=atan2(nonduct.rtd.rkxbig,nonduct.rtd.rkybig);
rkxybig=(nonduct.rtd.rkxbig.^2+nonduct.rtd.rkybig.^2).^0.5;
rkxybig_proj=rkxybig.*cos(phi_k-phi_xy);
k_unit_xy=rkxybig_proj./((rkxybig_proj.^2+nonduct.rtd.rkzbig.^2).^0.5);
k_unit_z=nonduct.rtd.rkzbig./((rkxybig_proj.^2+nonduct.rtd.rkzbig.^2).^0.5);
timing=0:nonduct.rtd.gtime(end)/10:nonduct.rtd.gtime(end);
timing_index=interp1(nonduct.rtd.gtime,(1:length(nonduct.rtd.rrr)),timing);
xy_timing=interp1((1:length(nonduct.rtd.rrr)),xy,timing_index);
z_timing=interp1((1:length(nonduct.rtd.rrr)),nonduct.rtd.z,timing_index);
k_unit_xy_timing=interp1((1:length(nonduct.rtd.rrr)),k_unit_xy,timing_index);
k_unit_z_timing=interp1((1:length(nonduct.rtd.rrr)),k_unit_z,timing_index);
for i=1:length(timing)
%for i=1:1
    hold on
%     line([xy(i) xy(i)+k_unit_xy(i)*0.15],[nonduct.rtd.z(i) nonduct.rtd.z(i)+k_unit_z(i)*0.15], ...
%         'Color','k','linewidth',1.5)
    line([xy_timing(i) xy_timing(i)+k_unit_xy_timing(i)*0.15],[z_timing(i) z_timing(i)+k_unit_z_timing(i)*0.15], ...
        'Color','k','linewidth',1.5)
end

title('Nonducted')
text(-0.7,1.5,'(a)','FontSize',16)
set(gca,'TickDir','out')
set(gca,'FontSize',12)

s2=subplot(2,2,2)
[duct.rtd.gtime,non_rep]=unique(duct.rtd.gtime);
duct.rtd.x=duct.rtd.x(non_rep);
duct.rtd.y=duct.rtd.y(non_rep);
duct.rtd.z=duct.rtd.z(non_rep);
duct.rtd.rrr=duct.rtd.rrr(non_rep);
duct.rtd.rkxbig=duct.rtd.rkxbig(non_rep);
duct.rtd.rkybig=duct.rtd.rkybig(non_rep);
duct.rtd.rkzbig=duct.rtd.rkzbig(non_rep);

xm=sqrt(duct.n3d.x.^2+duct.n3d.y.^2);
xy=sqrt(duct.rtd.x.^2+duct.rtd.y.^2);
pcolor(xm,duct.n3d.z,log10(duct.n3d.dent1))
shading flat
hold on
% contour(xm,hout.n3d.z,hout.n3d.flhr,[0 0]+freq0*1d-3,'m')
plot_Earth('2M','l',[1.5:0.5:3])
xlabel('R_x_y (R_E)')
ylabel('Z (R_E)')
axis equal
c2=colorbar('eastoutside');
colormap(jet)
c2.Label.String='log(density(m^-^3))';
caxis([8.5,10])
xlim([0,3])
ylim([-1.5,1.5])
hold on
plot(xy,duct.rtd.z,'color','m','linewidth',2); % ray path

phi_xy=atan2(duct.rtd.x,duct.rtd.y);
phi_k=atan2(duct.rtd.rkxbig,duct.rtd.rkybig);
rkxybig=(duct.rtd.rkxbig.^2+duct.rtd.rkybig.^2).^0.5;
rkxybig_proj=rkxybig.*cos(phi_k-phi_xy);
k_unit_xy=rkxybig_proj./((rkxybig_proj.^2+duct.rtd.rkzbig.^2).^0.5);
k_unit_z=duct.rtd.rkzbig./((rkxybig_proj.^2+duct.rtd.rkzbig.^2).^0.5);
timing=0:duct.rtd.gtime(end)/10:duct.rtd.gtime(end);
timing_index=interp1(duct.rtd.gtime,(1:length(duct.rtd.rrr)),timing);
xy_timing=interp1((1:length(duct.rtd.rrr)),xy,timing_index);
z_timing=interp1((1:length(duct.rtd.rrr)),duct.rtd.z,timing_index);
k_unit_xy_timing=interp1((1:length(duct.rtd.rrr)),k_unit_xy,timing_index);
k_unit_z_timing=interp1((1:length(duct.rtd.rrr)),k_unit_z,timing_index);
for i=1:length(timing)
%for i=1:1
    hold on
%     line([xy(i) xy(i)+k_unit_xy(i)*0.15],[duct.rtd.z(i) duct.rtd.z(i)+k_unit_z(i)*0.15], ...
%         'Color','k','linewidth',1.5)
    line([xy_timing(i) xy_timing(i)+k_unit_xy_timing(i)*0.15],[z_timing(i) z_timing(i)+k_unit_z_timing(i)*0.15], ...
        'Color','k','linewidth',1.5)
    
end

title('Ducted')
text(-0.7,1.5,'(b)','FontSize',16)
set(gca,'TickDir','out')
set(gca,'FontSize',12)


s3=subplot(2,2,3)
patch([-20 -20 10 10], [90 180 180 90], [0.8 0.8 0.8],'EdgeColor','none')
hold on
p31=plot(nonduct.rtd.lat,nonduct.rtd.psid,'b','Linewidth',1.2);
thetap=atan2(abs(nonduct.rtd.rvgp),nonduct.rtd.rvgz)*180/pi;
p32=plot(nonduct.rtd.lat,mod(thetap+360,360),'r','Linewidth',1.2);
p33=plot(nonduct.rtd.lat,nonduct.rtd.Gendrin,'k--','Linewidth',1.2);
ylim([90 180])
xticks(-50:10:50)
%xticklabels({'','-40','','-20','','0','','20','','40',''})
yticks(90:10:180)
yticklabels({'90','','','120','','','150','','','180'})
legend([p31 p32 p33],{'\theta_k','\theta_p','\theta_G'},'Location','southeast')
xlabel('MLAT (°)')
ylabel('Polar angle (°)')
text(-72,180,'(c)','FontSize',16)
set(gca,'FontSize',12)
set(gca,'TickDir','out')
set(gca, 'Layer', 'top')

s4=subplot(2,2,4)
patch([-20 -20 10 10], [90 180 180 90], [0.8 0.8 0.8],'EdgeColor','none')
hold on
p41=plot(duct.rtd.lat,duct.rtd.psid,'b','Linewidth',1.2);
thetap=atan2(abs(duct.rtd.rvgp),duct.rtd.rvgz)*180/pi;
p42=plot(duct.rtd.lat,mod(thetap+360,360),'r','Linewidth',1.2);
p43=plot(duct.rtd.lat,duct.rtd.Gendrin,'k--','Linewidth',1.2);
ylim([90 180])
xticks(-50:10:50)
yticks(90:10:180)
yticklabels({'90','','','120','','','150','','','180'})
legend([p41 p42 p43],{'\theta_k','\theta_p','\theta_G'},'Location','southeast')
xlabel('MLAT (°)')
ylabel('Polar angle (°)')
text(-72,180,'(d)','FontSize',16)
set(gca,'FontSize',12)
set(gca,'TickDir','out')
set(gca, 'Layer', 'top')

s1.Position=[0.1 0.45 0.3 0.47];
s2.Position=[0.58 0.45 0.3 0.47];
s3.Position=[0.1 0.1 0.3 0.25];
s4.Position=[0.58 0.1 0.3 0.25];