%% SAPS paper figure: 4x3: split from fig3_3, show the 2D status of ISp.
% 1: 2D ionospheric Enflux and convection with DMSP trajectory, 2: potential+fac, 3: TEC.
% Actual Figure 4 revision in paper, add SigP contour on top of TEC. 
% Use a separate row for WID.
close all; clear;
% addpath(genpath('~binzheng/matlab_source'));
% addpath('./subaxis/');
model_type = 'MIX';

%% Directory and runs
% rootdir = '/glade/scratch/ldong/'; subdir='LTR/';
srun0 = 'LTR_qq_W17';
% plt_dir0 = ['/glade/work/ldong/Plots/',subdir,srun0,'/'];
% plt_dir1 = [plt_dir0,'/SAPS_DMSP/'];
plt_dir = [pwd,'/'];%[plt_dir0,'/SAPS_Paper1/'];
if(~exist(plt_dir,'dir')); mkdir(plt_dir); end

efxns = cell(1,5); facns = efxns; potns = efxns; TEC_MIXs = efxns; SigP_MIXs = efxns;
Vdx_smps = efxns; Vdy_smps = efxns;
LTR_WIDs = efxns; LTR_efxs = efxns; LTR_orbits = efxns;
sruns = ['LTR_qq_W17_1000_1100'; 'LTR_qq_W17_1130_1230'; 'LTR_qq_W17_1500_1530'; 'LTR_qq_W17_2000_2030'; 'LTR_qq_W17_2130_2230'];
t_start0s= ['2013-03-17T10-12-00'; '2013-03-17T11-53-00'; '2013-03-17T15-17-00'; '2013-03-17T20-20-00'; '2013-03-17T22-02-00'];
t_stop0s = ['2013-03-17T10-22-00'; '2013-03-17T12-03-00'; '2013-03-17T15-27-00'; '2013-03-17T20-30-00'; '2013-03-17T22-12-00'];
t_starts = ['2013-03-17T10-17-00'; '2013-03-17T11-58-00'; '2013-03-17T15-21-00'; '2013-03-17T20-23-00'; '2013-03-17T22-05-00'];
% for ii=1:3
%     srun = sruns(ii,:); t_start0 = t_start0s(ii,:); t_stop0 = t_stop0s(ii,:);
%     hdf_fldr = [rootdir,subdir,srun0,'/'];
%     load([plt_dir1,srun,'_LTR_sampled_cmp_',t_start0,'_',t_stop0(end-7:end-3),'.mat']);%,...
%     %     'xi','yi','x_smp','y_smp','efxn','facn','potn','TEC_MIX','Vdx_smp','Vdy_smp',...
%     %     'LTR_WID','LTR_efx','LTR_orbit');
%     efxns{ii} = efxn; facns{ii} = facn; potns{ii} = potn; TEC_MIXs{ii} = TEC_MIX;
%     Vdx_smps{ii} = Vdx_smp; Vdy_smps{ii} = Vdy_smp; 
%     LTR_WIDs{ii} = LTR_WID; LTR_efxs{ii} = LTR_efx; LTR_orbits{ii} = LTR_orbit;
% end
for ii=1:3
    srun = sruns(ii,:); t_start0 = t_start0s(ii,:); t_stop0 = t_stop0s(ii,:);
%     hdf_fldr = [rootdir,subdir,srun0,'/'];
    load([srun,'_LTR_sampled_cmp_',t_start0,'_',t_stop0(end-7:end-3),'_MIX_TIE.mat']);% generated by SAPS_paper1_fig3_2.m
%         'xi','yi','x_smp','y_smp','efxn','facn','potn','TEC_MIX','Vdx_smp','Vdy_smp',...
%         'LTR_WID_MIX','LTR_efx_MIX','LTR_ped_MIX','LTR_hal_MIX',...
%         'LTR_WID_TIE','LTR_efx_TIE','LTR_ped_TIE','LTR_hal_TIE','LTR_fac','LTR_tec','LTR_orbit');
    efxns{ii} = efxn; facns{ii} = facn; potns{ii} = potn; TEC_MIXs{ii} = TEC_MIX; SigP_MIXs{ii} = SigP_MIX;
    Vdx_smps{ii} = Vdx_smp; Vdy_smps{ii} = Vdy_smp; 
    LTR_WIDs{ii} = LTR_WID_MIX; LTR_efxs{ii} = LTR_efx_MIX; 
    LTR_orbits{ii} = LTR_orbit;
end
% Load DMSP data.
dsec = 10;
load(['Vionh_F18_20130317.mat']);%,'ttime','Vionh','slat','slon','salt','ele_totl_engy');
DMSP_WID = Vionh(1:end); % m/s
DMSP_lat = slat(1:dsec:end-1); % Geo lat deg
DMSP_lon = slon(1:dsec:end-1); % Geo lon deg
DMSP_time=ttime(1:end); % epoch day
DMSP_efx = ele_totl_engy(1:1:end);
LTR_time=ttime(1:dsec:end);

%% Plot 2x2 2D: 1. West drift on enflux; 2. Phi contour on FAC; 3. Sigma_P; 4. NmF2.
qvfac=1/20; 
spc = 0.01; pdn = 0.01; mgn = 0.025; lwid=2; font=16; 
cl1=[0, 0.4470, 0.7410]; cl2='k';
Xlim=[-1,1]*0.866; Ylim=Xlim;
qminmax = [0.01, 0.99];
nrow=4; ncol=3;

fig = figure; 
set(fig,'DefaultAxesColorOrder',zeros(2,3));
ddftL = [20,40,20,30,-30]/86400;
ddftD = [-230,-180,-140,-10,-60]/86400;
hcells = cell(nrow,ncol);
% xcols = [0.0350, 0.6750, 0.6350];
% xcols1= [0.0464, 0.3577, 0.6615];
% yrows = [0.8621, 0.5864, 0.3107, 0.0350]-[0.0325, 0.0325, 0.03, 0.025];
dxcol = 0.3; % panel distance along x.
xplab = 1.2;
dxcol1 = 0.25; dyrow = 0.2407; dyrow1 = 0.1029+0.03;
panlabs = {'(a)','(b)','(c)';'(d)','(e)', '(f)';'(g)','(h)','(i)';'(j)', ...
    '(k)','(l)'};
for ii=1:ncol
    efxn = efxns{ii}; potn = potns{ii}; facn = facns{ii}; 
    LTR_orbit = LTR_orbits{ii}; LTR_WID = LTR_WIDs{ii}; LTR_efx = LTR_efxs{ii}; 
    TEC_MIX = TEC_MIXs{ii}; SigP_MIX = SigP_MIXs{ii};
    t_start0 = t_start0s(ii,:); t_stop0 = t_stop0s(ii,:);
    stime = datenum(t_start0,'yyyy-mm-ddTHH-MM-ss');
    etime = datenum(t_stop0 ,'yyyy-mm-ddTHH-MM-ss');
    tlim=[stime,etime]; ttick = (stime*1440:2:etime*1440)/1440;
    iplt1 = find(LTR_time<=etime&LTR_time>=stime);
    iplt2 = find(DMSP_time<=etime&DMSP_time>=stime);
    dft0L = datenum(t_starts(ii,:),'yyyy-mm-ddTHH-MM-ss')+ddftL(ii); i2dshowL= find(LTR_time==dft0L); % the moment shown in 2D.
    [pkefxL,ipkefxL] = max(LTR_efx(iplt1)); 
    [pkefxD,ipkefxD] = max(medfilt1(DMSP_efx(iplt2),3)); 
    jebL = 0; while(LTR_efx(iplt1(1)+jebL)<=0.1*pkefxL); jebL = jebL+1; end; iebL = iplt1(1) + jebL;
    jebD = 0; while(DMSP_efx(iplt2(1)+jebD)<=0.1*pkefxD); jebD = jebD+1; end; iebD = iplt2(1) + jebD;
    
    h1=subaxis(nrow,ncol,ii,'Spacing',spc,'Padding',pdn,'Margin',mgn); hold on; box on; grid on;
    hold on; box on; axis equal tight off; % Ion Drift
%     pcolor(xi,yi,efxn); shading flat; caxis(quantile(efxn(:),qminmax)); caxis([0,1.]);
%     quiver([x_smp(:);0.7],[y_smp(:);0.7],[Vdx_smp(:);0]*qvfac,[Vdy_smp(:);-2]*qvfac,'k','autoscale','off'); 
    quiver(x_smp(:),y_smp(:),Vdx_smp(:)*qvfac,Vdy_smp(:)*qvfac,'k','autoscale','off'); 
    hold on;
    plotion1('k',[50 60 70 80],'sat','MLT',1,font-2); view(-90,90);
    iorbit = LTR_orbit(:,1)~=0;
    plot(LTR_orbit(iorbit,1),LTR_orbit(iorbit,2),'k','LineWidth',2);
    [~,iorbit0] = min(abs(LTR_time-dft0L));
    plot(LTR_orbit(iorbit0,1),LTR_orbit(iorbit0,2),'rx','LineWidth',3,'MarkerSize',10);
    xlim(Xlim); ylim(Ylim);
    set(gca,'fontsize',font);
%     if ii==ncol
%         c2=colorbar; title(c2,{'EnFlux','[mW/m^2]'}); %set(c2,'position',[0.9662    yrows(2)    0.0108    dyrow-0.05]);
%         set(c2,'position',[0.9436    0.7205    0.0108    0.1907]);
%     end
%     ht2=title(datestr(dft0L,'HH:MM:ss')); %set(ht2,'position',[1.0,1.2,0])
    if ii==2
        ht1=text(0.78,0.85,'2 km/s','fontsize',font-4); 
        quiver(0.7,0.7,0*qvfac,-2*qvfac,'k','autoscale','off'); 
    end
    text(0,xplab,panlabs{1,ii},'FontSize',font);
    ht2=title([datestr(dft0L,'HH:MM:ss'),' UT']); %set(ht2,'position',[1.0,1.2,0])
    hcells{1,ii} = h1; set(h1,'position',[0.0350+(ii-1)*dxcol,0.7550,0.2900,0.2100]);
    
    h2=subaxis(nrow,ncol,ii+ncol,'Spacing',spc,'Padding',pdn,'Margin',mgn); hold on; box on; grid on;
    hold on; box on; axis equal tight off; % EnFlux
    pcolor(xi,yi,efxn); shading flat; caxis(quantile(efxn(:),qminmax)); caxis([0,1.]);
    hold on;
    plotion1('k',[50 60 70 80],'sat','MLT',1,font-2); view(-90,90);
%     iorbit = LTR_orbit(:,1)~=0;
%     plot(LTR_orbit(iorbit,1),LTR_orbit(iorbit,2),'k','LineWidth',2);
%     [~,iorbit0] = min(abs(LTR_time-dft0L));
    plot(LTR_orbit(iorbit0,1),LTR_orbit(iorbit0,2),'rx','LineWidth',3,'MarkerSize',10);
    xlim(Xlim); ylim(Ylim);
    set(gca,'fontsize',font);
    if ii==ncol
        c2=colorbar; title(c2,{'EnFlux','[mW/m^2]'}); %set(c2,'position',[0.9662    yrows(2)    0.0108    dyrow-0.05]);
        set(c2,'position',[0.92    0.5262    0.0108    0.18]);
    end
    text(0,xplab,panlabs{2,ii},'FontSize',font);
%     ht2=title(datestr(dft0L,'HH:MM:ss')); %set(ht2,'position',[1.0,1.2,0])
    hcells{2,ii} = h2; set(h2,'position',[0.0350+(ii-1)*dxcol,0.5150,0.2900,0.2100]);
    
    h3=subaxis(nrow,ncol,ii+ncol*2,'Spacing',spc,'Padding',pdn,'Margin',mgn); hold on; box on; grid on;
    hold on; box on; axis equal tight off; % Ion Drift on EnFlux
    pcolor(xi,yi,facn); shading flat; caxis(quantile(facn(:),qminmax)); caxis([-1,1]*0.6);
%     dplvl=10; plvl = (ceil(min(potn(:))/dplvl)*dplvl) : dplvl : (floor(max(potn(:))/dplvl)*dplvl);
%     contour(xi,yi,potn,plvl,'w','LineWidth',lwid/2,'LineStyle','--'); hold on;
    ipos = potn>0;
    potn_p = potn; potn_n = potn; potn_p(~ipos) = NaN; potn_n(ipos) = NaN;
    contour(xi,yi,potn_p,0:20:180,'k','LineWidth',lwid/2,'LineStyle','-'); hold on;
    contour(xi,yi,potn_n,-180:20:0,'k','LineWidth',lwid/2,'LineStyle','--'); hold on;
    plotion1('k',[50 60 70 80],'sat','MLT',1,font-4); view(-90,90);
    plot(LTR_orbit(iorbit0,1),LTR_orbit(iorbit0,2),'rx','LineWidth',3,'MarkerSize',10);
    xlim(Xlim); ylim(Ylim);
    set(gca,'fontsize',font);
    if ii==ncol
        c3=colorbar; hc3=title(c3,{'J_{//} [\muA/m^2]'}); %set(c3,'position',[0.9662    yrows(3)    0.0108    dyrow-0.05]);
        set(c3,'position',[0.92    0.2813    0.0108    0.18]);
    end
    ht2_1=text(-0.7, 1.5,{[num2str(min(potn(:)),'%.2f'),' kV'],...
        [num2str(max(potn(:)),'%.2f'),' kV']},'fontsize',font-4);
    ht2_2=text( 0.7, 1.5,{[num2str(min(facn(:)),'%.2f'),' \muA/m^2'],...
        [num2str(max(facn(:)),'%.2f'),' \muA/m^2']},'fontsize',font-4);
    text(0,xplab,panlabs{3,ii},'FontSize',font);
    hcells{3,ii} = h3; set(h3,'position',[0.0350+(ii-1)*dxcol,0.2750,0.2900,0.2100]);
    
    h4=subaxis(nrow,ncol,ii+ncol*3,'Spacing',spc,'Padding',pdn,'Margin',mgn); hold on; box on; grid on;
    hold on; box on; axis equal tight off; % Ion Drift on EnFlux
    pcolor(xi,yi,TEC_MIX*1e-12); shading flat; caxis(quantile(TEC_MIX(:)*1e-12,qminmax)); caxis([0,15]);
    contour(xi,yi,SigP_MIX,0:3:18,'k','LineWidth',lwid/2,'LineStyle','-'); hold on;
    plotion1('k',[50 60 70 80],'sat','MLT',1,font-4); view(-90,90);
    plot(LTR_orbit(iorbit0,1),LTR_orbit(iorbit0,2),'rx','LineWidth',3,'MarkerSize',10);
    xlim(Xlim); ylim(Ylim);
    set(gca,'fontsize',font);
    if ii==ncol
        c4=colorbar; title(c4,{'TEC','[TECU]'}); %set(c4,'position',[0.9662    yrows(4)    0.0108    dyrow-0.05]);
        set(c4,'position',[0.92    0.0238    0.0108    0.18]);
    end
    text(0,xplab,panlabs{4,ii},'FontSize',font);
    hcells{4,ii} = h4; set(h4,'position',[0.0350+(ii-1)*dxcol,0.0350,0.2900,0.2100]);
end
set(fig,'position',[1028          49         829         804]);
% print(fig,[plt_dir,'fig03_LTR_DMSP_MIXTIE_r2'],'-depsc');
% close(fig);
%% END of CODE