%% Plot figure 4 of geospace concussion.
% Version 2, move time series to the top.
% Three snapshots of WI_ExB in a row on top.
% two columns of time history for two LT, one at PGR, the other at HOK.
% In each column from top to bottom are REMIX keo of Zonal E field, WI, and
% time series of SuperDARN and TIEGCM sampling at that beam.

close all; clear;
strdatefmt = 'yyyy-mm-ddTHH-MM-ss';
MJD0 = datenum('2000-01-01T00-00-00',strdatefmt)-51544;
constants;
addpath('~ldong/mcode/color_brewer/'); %colormap(brewermap([],'YlGnBu'))
dash_dir = '/glade/work/ldong/DASH/Shi21_concussion/';
plt_dir = dash_dir;

    % Prepare for variables to be written to nc file.
    ncname = [dash_dir,'fig3_WI.nc']; % generated with SAPS_paper2_fig5.m and SAPS_paper2_fig5_plt.m
    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
    
    %% Load all data.
    t_start = '2011-10-24T18-28-00';
    t_stop  = '2011-10-24T18-38-00';
    smin = datenum(t_start,strdatefmt)*1440;
    emin = datenum(t_stop ,strdatefmt)*1440;
    imins = smin+[1.5,3,5];
    strdates = {};
    for ii=1:3
        imin = imins(ii);
        strdate = datestr(imin/1440,strdatefmt);
        strdates{ii} = strdate;
    end
    
        UTdays_sd = {UTdays_sd_001,UTdays_sd_002};
        WI_SD_sd = {WI_SD_sd_001,WI_SD_sd_002};
        stats = {'PGR B12 Median','HOK B04 Median'};
        MLAT_tie_plt = {MLAT_tie_plt_001,MLAT_tie_plt_002};
        Vper_plt = {Vper_plt_001,Vper_plt_002};
        MLAT_mix_plt = {MLAT_mix_plt_001,MLAT_mix_plt_002};
        Ezon_plt = {Ezon_plt_001,Ezon_plt_002};
        WI_plts = {WI_plts_001,WI_plts_002,WI_plts_003};
        
        SLT = LON/15+12;
        p_slts = mod(p_lons/15+(imin-floor(imin/1440)*1440)/60,24);
        p_lonp = p_lons*0;
        for ii=1:2
            [~,islt] = min(abs(SLT(:,1)-p_slts(ii)));
            p_lonp(ii) = LON(islt,1);
        end
    %% Plot
    % for 2d WI plot
    slt2lon = 15;
    latlim = [-1,1]*90;
    sltlim = [0,24];
    lonlim = sltlim*slt2lon-180;
    lattick = -60:30:60;
    slttick = 3:3:21;
    lontick = slttick*slt2lon-180;
    sltticklab = {};
    for ii=1:length(slttick)
        sltticklab{ii} = num2str(slttick(ii),'%02d');
    end
    ttls = strdates;
    latlab = 'MLat [deg]';
    cttl_wi = 'V_{Vertical} [m/s]';
    wilim = [-1,1]*100;
    cl_grd = 'k';
    % for keo plot.
    cttl_eaz = 'E_{Zonal} [mV/m]';
    eazlim = [-1,1]*50;
    cmap_WI = '*RdBu';
    cmap_Ea = '*RdBu';
    tlim = [smin,emin]/1440;
    ttick = (smin:2:emin)/1440;
    mlatlim = [45,90];
    mlattick = 50:10:80;
    % for SuperDARN sampling
    cl_wi_ts = 'k';
    wilim_sd = [-1,1]*50;
    cl_sds = {'m',cl_gre}; % [0.9290 0.6940 0.1250]
%     cl_sds = {[0.4940 0.1840 0.5560], [0.4660 0.6740 0.1880]};
    mk_sds = {'^','o'};
    msize = 10;
    dt0 = 0.5/1440;
    
    xpan1 = -170;
    ypan1 = 75;
    xlab1 = tlim(1)+0.02*(tlim(2)-tlim(1));
    ylab1 = 38;
    xpan = xlab1;
    ypan = 85;
%     panlabs = {'(a)','(b)','(c)','(d)','(f)','(h)','(e)','(g)','(i)'};
    panlabs = {'(a)','(c)','(e)','(b)','(d)','(f)','(g)','(h)','(i)'};
    
    lwid = 2; font = 14;
    hps = {}; cls = {};
    icls = 1;
    
    pngname = ncname(1:end-3);
    disp(['Plotting ',pngname]);
    nrow = 5; ncol = 3; fpos = [10,10,1200,820];
    fig = figure('position',fpos);
    
    % bottom left 3x2
    isub = 0;
    nrow = 5; ncol = 2;
    for icol=1:2
        % plot SuperDARN sampling.
        isub = isub+1;
        hps{isub} = subplot(nrow,ncol,icol); %disp([nrow,ncol,icol]);
        hold on; box on; grid on;
        set(gca,'Fontsize',font,'LineWidth',lwid);
        plot(UTdays_ts,WI_SD_ts(:,icol),'color',cl_wi_ts,'LineWidth',lwid); 
        plot(UTdays_sd{icol}-dt0,-WI_SD_sd{icol},'color',cl_sds{icol},'LineWidth',lwid*1.2);
        xlim(tlim);
%         ylim(wilim);
        set(gca,'xtick',ttick,'xticklabel','');
        if icol==1
            ylabel(cttl_wi);
        end
        hlg = legend({'TIEGCM',stats{icol}(1:7)});
        if icol==1
            set(hlg,'location','southeast');
        end
%         text(xlab1,ylab1,['MLT=',num2str(MLT_SD(icol),'%.1f')],'Fontsize',font);
        ylime = get(gca,'ylim');
        text(xpan,ylime(1)+0.9*(ylime(2)-ylime(1)),panlabs{isub},'Fontsize',font);
        
        % plot remix keo of WI.
        isub = isub+1;
        hps{isub} = subplot(nrow,ncol,icol+ncol); %disp([nrow,ncol,icol+ncol]);
        hold on; box on; grid on;
        set(gca,'Fontsize',font,'LineWidth',lwid);
        colormap(hps{isub},brewermap([],cmap_WI)); 
        pcolor(UTdays_tie,MLAT_tie_plt{icol},Vper_plt{icol}); 
        shading flat; 
        if icol==2
            cls{icls} = colorbar;
            ylabel(cls{icls},cttl_wi);
            icls = icls+1;
        end
        caxis(wilim);
        xlim(tlim);
        ylim(mlatlim);
        set(gca,'xtick',ttick,'xticklabel','','ytick',mlattick);
        plot_grid(ttick,mlattick,cl_grd);
        if icol==1
            ylabel(latlab);
        end
%         text(xlab1,ylab1,['MLT=',num2str(MLT_SD(icol),'%.1f')],'Fontsize',font);
        text(xpan,ypan,panlabs{isub},'Fontsize',font);
        text(xlab1,ylab1,['MLT=',num2str(MLT_SD(icol),'%.1f')],'Fontsize',font);
        
        % plot tie keo of Eazi
        isub = isub+1;
        hps{isub} = subplot(nrow,ncol,icol+ncol*2); %disp([nrow,ncol,icol+ncol*2]);
        hold on; box on; grid on;
        set(gca,'Fontsize',font,'LineWidth',lwid);
        colormap(hps{isub},brewermap([],cmap_Ea)); 
        pcolor(UTdays_mix,MLAT_mix_plt{icol},Ezon_plt{icol}); 
        shading flat; 
        if icol==2
            cls{icls} = colorbar;
            ylabel(cls{icls},cttl_eaz);
            icls = icls+1;
        end
        caxis(eazlim);
        xlim(tlim);
        ylim(mlatlim);
        set(gca,'xtick',ttick,'ytick',mlattick);
        datetick('x','HH:MM','keepticks');
        plot_grid(ttick,mlattick,cl_grd);
%         xlabel([t_start(1:10),' UT']);
        if icol==1
            ylabel(latlab);
        end
        text(xpan,ypan,panlabs{isub},'Fontsize',font);
    end
    
    % first row of WI(GLat, SLT)
    nrow = 5; ncol = 3;
    for icol=1:3
        isub = icol+6;
        hps{isub} = subplot(nrow,ncol,icol+9+[0,ncol]);
        hold on; box on; grid on; axis tight; % equal
        set(gca,'Fontsize',font,'Linewidth',lwid);
        colormap(brewermap([],cmap_WI));
        pcolor(LON,LAT,WI_plts{icol});
        shading flat; 
%         cls{isub} = colorbar; 
        caxis(wilim);
        xlim(lonlim);
        ylim(latlim);
%         xlabel('SLT [Hour]');
        if icol==1
            ylabel('Geo Lat [deg]');
            set(gca,'ytick',lattick);
        else
            set(gca,'ytick',lattick,'yticklabel','');
        end
        set(gca,'xtick',lontick,'ytick',lattick,'xticklabel',sltticklab);
        plot_grid(lontick,lattick,cl_grd);
        title(ttls{icol});
        if icol==3
            cls{icls} = colorbar;
            ylabel(cls{icls},cttl_wi);
            icls = icls+1;
        end
        % plot SD location
        for isd=1:length(p_lons)
            plot(p_lonp(isd),p_lats(isd),'color',cl_sds{isd},'Marker',mk_sds{isd},'Markersize',msize,'LineWidth',lwid*1.5);
        end
        if icol==1
            text(lonlim(1)-0.1*(lonlim(2)-lonlim(1)),-108,'SLT','Fontsize',font);
        end
        text(xpan1,ypan1,panlabs{isub},'Fontsize',font);
    end
    
    %
% %     for isub=1:2
% %         hpos = get(hps{isub},'position');
% %         set(hps{isub},'position',hpos+[0,0.1,0,0]);
% %     end
% %     set(hps{3},'position',[0.6916,0.7280,0.2134,0.2970]);
    
    for isub=1:6
        [irow,icol] = ind2sub([3,2],isub);
        hpos = get(hps{isub},'position');
        set(hps{isub},'position',hpos+[0,0.05*(3-irow)-0.1,0.3347-hpos(3),0.05]);
    end
    
    xpos = [0.13, 0.397, 0.665];
    ypos = 0.05;
    dpx = 0.24;
    dpy = 0.20;
    for isub=7:9
        set(hps{isub},'position',[xpos(isub-6)    ypos    dpx    dpy]);
    end
    
    cpos1 = get(cls{1},'position');
    cpos2 = get(cls{2},'position');
    set(cls{2},'position',cpos2+[cpos1(1)-cpos2(1),0,0,0]);
    
    set(gcf,'renderer','painters');
    print(fig,pngname,'-dpng');
    close(fig);
    