% changed from SAPS_paper2_fig5_plt2.m, add IMF/Dst above with four vertical
%% test written nc files.
close all; clear;

addpath('~ldong/mcode/color_brewer/'); %colormap(brewermap([],'YlGnBu'))
dash_dir = '/glade/work/ldong/DASH/Lin20_DEP_SAPS/';
plt_dir = dash_dir;

    % Prepare for variables to be written to nc file.
    ncname = [dash_dir,'fig3_tail.nc'];
    ncid = netcdf.open(ncname,'NC_NOWRITE'); 
    info = ncinfo(ncname);
    varnames = {}; % store all available variable names.
    for ii=1:size(info.Variables,2)
        varnames{ii} = info.Variables(ii).Name;
        disp([num2str(ii,'%02d'),' ',info.Variables(ii).Name]); 
    end
    
    for ii=1:length(varnames)
        varid = netcdf.inqVarID(ncid,varnames{ii});   
        eval([varnames{ii},' = netcdf.getVar(ncid,varid);']);
    end
    
    WID_clines_eqs_06UT = {WID_clines_eqs_06UT_01,WID_clines_eqs_06UT_02}; 
    WID_clines_eqs_07UT = {WID_clines_eqs_07UT_01,WID_clines_eqs_07UT_02}; 
    WID_clines_eqs_09UT = {WID_clines_eqs_09UT_01,WID_clines_eqs_09UT_02,WID_clines_eqs_09UT_03,WID_clines_eqs_09UT_04}; 
    WID_clines_eqs_11UT = {WID_clines_eqs_11UT_01,WID_clines_eqs_11UT_02};
    WID_clines_eqs = {WID_clines_eqs_06UT,WID_clines_eqs_07UT,WID_clines_eqs_09UT,WID_clines_eqs_11UT}; 
    FAC_clines_eqs_06UT = {FAC_clines_eqs_06UT_01,FAC_clines_eqs_06UT_02};
    FAC_clines_eqs_07UT = {FAC_clines_eqs_07UT_01,FAC_clines_eqs_07UT_02};
    FAC_clines_eqs_09UT = {FAC_clines_eqs_09UT_01,FAC_clines_eqs_09UT_02,FAC_clines_eqs_09UT_03,FAC_clines_eqs_09UT_04,FAC_clines_eqs_09UT_05,FAC_clines_eqs_09UT_06};
    FAC_clines_eqs_11UT = {FAC_clines_eqs_11UT_01,FAC_clines_eqs_11UT_02};
    FAC_clines_eqs = {FAC_clines_eqs_06UT,FAC_clines_eqs_07UT,FAC_clines_eqs_09UT,FAC_clines_eqs_11UT};
    EFX_clines_eqs_06UT = {EFX_clines_eqs_06UT_01};
    EFX_clines_eqs_07UT = {EFX_clines_eqs_07UT_01,EFX_clines_eqs_07UT_02};
    EFX_clines_eqs_09UT = {EFX_clines_eqs_09UT_01};
    EFX_clines_eqs_11UT = {EFX_clines_eqs_11UT_01,EFX_clines_eqs_11UT_02};
    EFX_clines_eqs = {EFX_clines_eqs_06UT,EFX_clines_eqs_07UT,EFX_clines_eqs_09UT,EFX_clines_eqs_11UT};

    %%
    t_start = '2013-03-17T00-00-00';
    strdates = {'2013-03-17T06-00-00','2013-03-17T07-00-00','2013-03-17T09-00-00','2013-03-17T11-00-00'};
    
% dashed lines to indicate the times.
    var_limMs = {[-2,2], [-1,1]*2,[-1,1],[0,10]};
    plvl = [0,5:5:50];
    widlvl = [1,1]*0.5;
    faclvl = [1,1]*0.2;
    efxlvl = [1,1]*1.0;
    var_lvls = {plvl,widlvl,faclvl,efxlvl};
    ivec = 1;
    arrl0 = 10; % in x/y
    vecl0s = [500,10]; % 500 km/s, 10 mV/m
    vscales = arrl0./vecl0s;
    Xlim = [-45,15];
    Ylim = [-1,1]*25;
    dxyt = 10;
    dxyv = 2;
    xtick = Xlim(1):dxyt:Xlim(2);
    ytick = Ylim(1):dxyt:Ylim(2);
    xv0 = xtick(1); % location to plot vector scale length.
    yv0 = ytick(1);
    vttls = {'km/s','mV/m'};
    cttls = {'log_{10} P [nPa]','WID [km/s]','FAC [\mu A/m^2]','EnFlux [mW/m^2]'};
    xlab = 'X [R_E]';
    ylab = 'Y [R_E]';
    Xlimi = [-1,1]*0.8;
    Ylimi = [-1,1]*0.8;
    tlim = [4,14];
    ttick = 4:1:14;
    tcols = [6,7,9,11];
    xtext = 13.5;
    imflim = [-1,1]*20;
    dstlim = [-150,50];
    imftick = -20:10:20;
    dsttick = -150:50:50;
    
    plt_grd = 0;
    plt_con = [[1,0,0,0]; ... % p
        [1,1,0,0]; ... % wid
        [1,0,1,0]; ... % fac
        [1,1,1,1]]; % efx

    cl_grd = 'k';
    cl_fac = [0.4660, 0.6740, 0.1880];
    cl_wid = 1-cl_fac;
    cl_efx = 'm';
    cl_vec = 'k';
    cl_peq = 'k';
    cl_cons= {cl_peq,cl_wid,cl_fac,cl_efx};
    cl_lat1= [0.8500, 0.3250, 0.0980];
    cmaps = {'*RdYlBu','*RdBu','*PuOr','*YlGnBu'};
    lats = 50:10:90;
    hps = {};
    cls = {};
    xttl = 18;
    yttl = 0;
    ltlabss = {[],[1,1,0,1],[0,1,0,1],[0,1,1,1]};
    font = 14; lwid = 2;
    pngname = [plt_dir,'fig3_tail'];
    disp(['Plotting ',pngname]);
    panlabs = {'(a)','(b)','(c)','(d)',...
        '(e)','(f)','(g)','(h)',...
        '(i)','(j)','(k)','(l)',...
        '(m)','(n)','(o)','(p)'};
    xp1 = 10;
    yp1 = 25;
    xp2 = 0.6;
    yp2 = 0.75;
    xc1 = -0.6;
    yc1 = -0.7;
    conlabs = {[num2str(widlvl(1),'%.1f'),' km/s'],[num2str(faclvl(1),'%.1f'),' \muA/m^2'],[num2str(efxlvl(1),'%.1f'),' mW/m^2']};
    concls = {cl_wid,cl_fac,cl_efx};
    
    nrow = 5; ncol = 4; fpos = [1,39,1160,830];
    fig = figure('position',fpos,'visible','on');
    
    hp0=subplot(nrow,ncol,1:ncol);
    yyaxis left;
    hold on; box on; grid on;
    set(gca,'Fontsize',font-2,'LineWidth',lwid);
    plot(Hours,IMFBy,'g-','LineWidth',lwid); 
    plot(Hours,IMFBz,'b-','LineWidth',lwid);
    plot(xlim,[0,0],'k--');
    text(xtext,15,'B_Y','color','g','Fontsize',font);
    text(xtext, 5,'B_Z','color','b','Fontsize',font);
    xlim(tlim);
    ylim(imflim);
    xlabel([t_start(1:10),' UT [Hours]']);
    ylabel('IMF');
    set(gca,'xtick',ttick,'ytick',imftick);

    yyaxis right;
    hold on; box on; grid on;
    set(gca,'Fontsize',font-2,'LineWidth',lwid);
    plot(Hours,SYMH,'r','LineWidth',lwid);
    plot(xlim,[0,0],'r--');
    for ii=1:length(tcols)
        plot([1,1]*tcols(ii),dstlim,'k-','LineWidth',lwid+1);
    end
    xlim(tlim);
    ylim(dstlim);
    xlabel('2013-03-17 UT','Fontsize',font-2);
    ylabel('SYMH [nT]');
    set(gca,'xtick',ttick,'ytick',dsttick);
    
    for icol = 1:ncol
        strdate = strdates{icol};
        if icol==1
            ttl = strdate;
        else
            ttl = strdate(12:16);
        end
        ttl = [strrep(strdate(12:13),'-',':'),' UT'];
        peq = peqsa(:,:,icol); 
        V_ExB_az = V_ExB_azsa(:,:,icol);
        facn = facnsa(:,:,icol);
        efxn = efxnsa(:,:,icol);
        var_pltMs = {log10(peq),V_ExB_az,facn,efxn};
        vxeq = vxeqsa(:,:,icol);
        vyeq = vyeqsa(:,:,icol);
        bzeq = bzeqsa(:,:,icol);
        if 1
            WID_clines_eq = WID_clines_eqs{icol};
            FAC_clines_eq = FAC_clines_eqs{icol};
            EFX_clines_eq = EFX_clines_eqs{icol};
        end
        
        if 1
            [xvec,yvec] = meshgrid(Xlim(1):2.5:0,-15:2.5:15);
            vxvec = griddata(xe,ye,vxeq,xvec,yvec);
            vyvec = griddata(xe,ye,vyeq,xvec,yvec);
%             exvec = griddata(xe,ye,exeq,xvec,yvec);
%             eyvec = griddata(xe,ye,eyeq,xvec,yvec);
        end
        
        Vecxs = {vxvec};
        Vecys = {vyvec};
        
        for irow=1
            isub = (irow-1)*ncol+icol;
            hps{isub} = subplot(nrow,ncol,isub+ncol);
            hold on; box on; grid on; axis equal tight; view(-90,90);
            set(gca,'Fontsize',font-2,'LineWidth',lwid);
            colormap(hps{isub},brewermap([],cmaps{irow}));
            pcolor(xe,ye,var_pltMs{irow});
            shading interp;
            if icol==ncol
                cls{irow}=colorbar;
            end
            caxis(var_limMs{irow});
            if irow==1
                contour(xe,ye,peq,plvl,cl_peq);
            end
            if ivec && irow==1 % no vec for rho.
                quiver(xvec,yvec,Vecxs{ivec}*vscales(ivec),Vecys{ivec}*vscales(ivec),'color',cl_vec,'LineWidth',lwid*0.6,'autoscale','off');
            end
            if ivec && irow==1 && icol==1 % no vec for rho.
                quiver(xv0,yv0,arrl0,0,'color',cl_vec,'Linewidth',lwid*3,'autoscale','off');
                text(xv0,yv0-7,[num2str(vecl0s(ivec),'%.0f'),' ',vttls{ivec}],'color',cl_vec,'Fontsize',font,'fontweight','bold','rotation',90);
            end
            if plt_grd
                plot_grid(xtick,ytick,cl_grd);
            end
            for iline=1:length(FAC_clines_eq)
                temp = FAC_clines_eq{iline};
                plot(temp(:,1),temp(:,2),'color',cl_fac,'LineWidth',lwid*0.6)
            end
            for iline=1:length(EFX_clines_eq)
                temp = EFX_clines_eq{iline};
                plot(temp(:,1),temp(:,2),'color',cl_efx,'LineWidth',lwid*0.6)
            end
            for iline=1:length(WID_clines_eq)
                temp = WID_clines_eq{iline};
                plot(temp(:,1),temp(:,2),'color',cl_wid,'LineWidth',lwid*0.6)
            end
%             contour(xe,ye,bzeq,[1,1]*0.001,'LineWidth',lwid);
            xlim(Xlim);
            ylim(Ylim);
            set(gca,'xtick',xtick,'ytick',ytick);
            if icol==1
                xlabel(xlab);
                hyl=ylabel(ylab); % [-54.7828   -0.0000   -1.0000]
                set(hyl,'position',[-47.3 43   -1.0000]);
            end
            if icol==ncol
                ylabel(cls{irow},cttls{irow});
            end
            if irow==1
                htt = title(ttl);
                set(htt,'position',[xttl,yttl,0]);
            end
            text(xp1,yp1,panlabs{isub},'Fontweight','Bold','Fontsize',font);
        end
        
        % 2-4 row
        for irow=2:nrow-1
            isub = (irow-1)*ncol+icol;
            hps{isub} = subplot(nrow,ncol,isub+ncol);
            hold on; axis equal tight off; view(-90,90);
            set(gca,'Fontsize',font,'LineWidth',lwid);
            colormap(hps{isub},brewermap([],cmaps{irow}));
            pcolor(xi,yi,var_pltMs{irow}); 
            for iisub=2:nrow-1
                if plt_con(iisub,irow)
                contour(xi,yi,var_pltMs{iisub},var_lvls{iisub},'color',cl_cons{iisub},'LineWidth',lwid);
                end
            end
            shading interp;
            if icol==ncol
                cls{irow}=colorbar;
            end
            caxis(var_limMs{irow});
            plotion(cl_lat1,lats,'sat','MLT',1,12,ltlabss{irow});
            xlim(Xlimi);
            ylim(Ylimi);
            if icol==ncol
                ylabel(cls{irow},cttls{irow});
            end
            text(xp2,yp2,panlabs{isub},'Fontweight','Bold','Fontsize',font);
            if icol==1
                text(xc1,yc1,conlabs{irow-1},'color',concls{irow-1},'Fontweight','Bold','Fontsize',font-2);
            end
        end
    end
        
    %%
    pos = get(hp0,'position');
    set(hp0,'position',pos+[0,0.06,0,0]);
    
    dpx = -0.05;
%     dpy0= [-0.05,-0.12,-0.11,-0.10];
    dpy0= [-0.04,-0.10,-0.10,-0.10];
    dpy = oneD(repmat(dpy0,nrow-1,1));
    dwx = 0.07;
    dwy = 0.07;
    for isub=1:nrow*ncol-ncol
        pos=get(hps{isub},'position');
        set(hps{isub},'position',pos+[dpx,dpy(isub),dwx,dwy]);
    end
        
    for irow=1:nrow-1
        isub = irow*ncol;
        pos=get(cls{irow},'position');
        set(cls{irow},'position',pos+[0.01,0.025,-0.005,-0.05]);
        pos=get(hps{isub},'position');
        set(hps{isub},'position',pos+[0.03,0,0,0]);
    end
    
    set(gcf,'renderer','painters');
    print(fig,pngname,'-dpng');
    print(fig,pngname,'-depsc');
%% END of CODE