%% Plot Figure 2 of MAGE-DMSP comparison. nc files generated in SAPS_paper2_fig2.
close all; clear;

strdatefmt = 'yyyy-mm-ddTHH-MM-ss';
dash_dir = '/glade/work/ldong/DASH/Lin20_DEP_SAPS/';
plt_dir = dash_dir;
ttimeAd0 = datenum('20130317','yyyymmdd');
Je2PhiE = 1.6e-12*pi;
SCs={'F16', 'F17', 'F18'};
% matname = [dash_dir,'fig2_mage_dmsp.mat'];
% save(matname,'ttimeA','DMSP_efxs','LTR_efxs','DMSP_Vihs','LTR_Vihs','DMSP_facs','LTR_facs','LTR_peds');

    % Prepare for variables to be written to nc file.
    ncname = [dash_dir,'fig2_mage_dmsp.nc'];
    
    ncid = netcdf.open(ncname,'NC_NOWRITE'); 
    varid = netcdf.inqVarID(ncid,'ttimeA');   
    ttimeA = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'DMSP_efxs');   
    DMSP_efxs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'DMSP_facs');   
    DMSP_facs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'DMSP_Vihs');   
    DMSP_Vihs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'LTR_efxs');   
    LTR_efxs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'LTR_facs');   
    LTR_facs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'LTR_Vihs');   
    LTR_Vihs = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'LTR_peds');   
    LTR_peds = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'MLAT');   
    MLAT = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'MLT');   
    MLT = netcdf.getVar(ncid,varid);
    
        %% Plot frame
        % Variables
    % 1. determine time range for one column.
    twin=10; 
    SCts = {...
            3,'2013-03-17T08-28-00','2013-03-17T08-32-00','2013-03-17T08-34-00'; ...
            3,'2013-03-17T11-18-00','2013-03-17T11-24-00','2013-03-17T11-26-00'; ...
            3,'2013-03-17T13-00-00','2013-03-17T13-06-30','2013-03-17T13-08-30'; ...
            };
    jdowns = {'2013-03-17T08-30-30','2013-03-17T08-33-10'; ...
            '2013-03-17T11-23-20','2013-03-17T11-26-50'; ...
            '2013-03-17T13-06-10','2013-03-17T13-09-00';};
    tjdowns = reshape(datenum(jdowns(:),strdatefmt),size(jdowns));

    % Parameters
    nrow=4; ncol=3;
    fpos = [1,35,1294,830];
    pngname = [plt_dir,'fig2_mage_dmsp'];
    disp(['Plotting ',pngname]);
    lwid=2; 
    lwidr = lwid/1.0;
    msize=8;
    font=14; 
    efxlim1 = [0,25]; 
    widlim1 = [-1,1]*2.5;
    cl_efx = 'm';
    cl_fac = [0.4660, 0.6740, 0.1880];
    cl_saps= 'b';
    cl_aubd= [0.3010, 0.7450, 0.9330];%[0.4940, 0.1840, 0.5560];
    cl_mage= [0.8500, 0.3250, 0.0980];
    ylabEfx = {'EnFlux','[mW/m^2]'};
    ylabVhd = {'V_{HD}','[km/s]'}; 
    ylabFac = {'FAC','[\mu A/m^2]'}; 
    ylabPed = {'\Sigma_P','[S]'};
    yrlim = 0.85;
    lrlim = 0.03;
    trlim = 0.1;
    panlabs = {'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)','(i)','(j)','(k)','(l)'};
    panlabs = reshape(panlabs,nrow,ncol);
    % Plot
    fig = figure('position',fpos,'visible','on');
    for icol=1:ncol
        iSC = SCts{icol,1};
        tmin = datenum(SCts{icol,2},strdatefmt);
        tmlocal1 = SCts{icol,3}; % to determine SAPS location as local peak.
        tmlocal2 = SCts{icol,4};
        iredi = localmaxr(ttimeA,LTR_efxs(:,iSC),datenum(tmlocal1,strdatefmt),datenum(tmlocal2,strdatefmt),0.1);
        iredj = localmaxr(ttimeA,DMSP_efxs(:,iSC),datenum(tmlocal1,strdatefmt),datenum(tmlocal2,strdatefmt),0.1);
        tmax = tmin+twin/1440;
        tlim = [tmin,tmax]; 
        implt = find(ttimeA<=tmax & ttimeA>=tmin); 
        if(abs(MLAT(implt(end),iSC))>abs(MLAT(implt(1),iSC)))
            xl = ttimeA(implt(1));
        else
            xl = ttimeA(implt(end));
        end
        xr = ttimeA(iredi);
        itick = find(mod(round((ttimeA-ttimeAd0)*86400),120)==0 & ttimeA<=tlim(2) & ttimeA>=tlim(1));
        ttick = ttimeA(itick);
        tlabs = datestr(ttick,'HH:MM');
        tlab1 = {};
        for ii=1:size(tlabs,1)
            tlab1{ii} = tlabs(ii,:);
        end
        tlab2 = compose('%.1f',MLAT(itick,iSC)');
%         tlab2{1} = ['MLAT ',tlab2{1}];
        tlab3 = compose('%.1f',MLT(itick,iSC)');
        labelArray = [tlab1; tlab2; tlab3]; 
        tlabs = strtrim(sprintf('%s\\newline%s\\newline%s\n', labelArray{:}));
        tlab0 = {'UT','MLAT','MLT'};
        
        xtext=tlim(2)-1.5/1440; 
    
        hp1=subplot(nrow,ncol,icol); hold on; box on; grid on;
        set(gca,'LineWidth',lwid,'Fontsize',font);
        hl(1)=plot(ttimeA(implt),DMSP_efxs(implt,iSC)*Je2PhiE,'k','LineWidth',lwid);
        hl(2)=plot(ttimeA(implt) ,LTR_efxs(implt,iSC),'color',cl_mage,'LineWidth',lwid); %colors{iSC},
        plot([1,1]*ttimeA(iredi),ylim,'color',cl_aubd,'Linestyle','--','LineWidth',lwidr); disp([MLAT(iredi,iSC),MLAT(iredj,iSC)])
        temp = real([DMSP_efxs(implt,iSC)*Je2PhiE;LTR_efxs(implt,iSC)]);
        ylimp = [min(temp),max(temp)*1.05];
        ylim(ylimp);
        patch([xl,xl,xr,xr],ylimp([1,2,2,1]),cl_efx,'FaceAlpha',0.2);
        set(gca,'xtick',ttick,'xticklabel',''); %datetick('x','HH:MM','keepticks');
        xlim(tlim); 
%         ylim(efxlim1); 
        if icol==1 
            ylabel(ylabEfx,'Fontsize',font);
        end
        ht0 = title(datestr(ttimeA(iredi),'yyyy-mm-dd/HH:MM:ss'));
        htt1=text(xtext,ylimp(1)+yrlim*(ylimp(2)-ylimp(1)),panlabs{1,icol},'Fontsize',font,'fontweight','bold');
        hl(3)=plot([xl,xr],[1,1]*ylimp(1)+lrlim*(ylimp(2)-ylimp(1)),'color',cl_efx,'LineWidth',4);
        htt2a = text(min(xl,xr),ylimp(1)-trlim*(ylimp(2)-ylimp(1)),'SubAuroral','color',cl_efx,'Fontsize',font,'fontweight','bold');
        if icol==2
            hlg_ab=legend(hl(1:2),{SCs{iSC},'MAGE'},'Location','Northwest','Fontsize',font,'Box','on');
        end
        
        hp2=subplot(nrow,ncol,ncol+icol); hold on; box on; grid on;
        set(gca,'LineWidth',lwid,'Fontsize',font);
        plot(ttimeA(implt),DMSP_Vihs(implt,iSC)*1e-3,'k','LineWidth',lwid);
        plot(ttimeA(implt),LTR_Vihs(implt,iSC)*1e-3,'color',cl_mage,'LineWidth',lwid); %,colors{iSC}
        plot([1,1]*ttimeA(iredi),ylim,'color',cl_aubd,'Linestyle','--','LineWidth',lwidr);
        temp = real([DMSP_Vihs(implt,iSC);LTR_Vihs(implt,iSC)])*1e-3;
        ylimp = [min(temp),max(temp)]*1.2;
        ylim(ylimp); 
        patch([xl,xl,xr,xr],ylimp([1,2,2,1]),cl_efx,'FaceAlpha',0.2);
        patch(tjdowns(icol,[1,1,2,2]),ylimp([1,2,2,1]),cl_fac,'FaceAlpha',0.2);
        xl1 = max(min(xl,xr),min(tjdowns(icol,:)));
        xr1 = min(max(xl,xr),max(tjdowns(icol,:)));
%         plot([xl1,xr1],[1,1]*ylimp(1)+lrlim*(ylimp(2)-ylimp(1)),'color',cl_saps,'LineWidth',4);
        set(gca,'xtick',ttick,'xticklabel',''); %datetick('x','HH:MM','keepticks');
        xlim(tlim); 
%         ylim(widlim1); 
        if icol==1
            ylabel(ylabVhd,'Fontsize',font);
        end
        htt2=text(xtext,ylimp(1)+yrlim*(ylimp(2)-ylimp(1)),panlabs{2,icol},'Fontsize',font,'fontweight','bold');
        htt21 = text(xl1,ylimp(1)-trlim*(ylimp(2)-ylimp(1)),'SAPS','color',cl_saps,'Fontsize',font,'fontweight','bold');
        
        hp3=subplot(nrow,ncol,ncol*2+icol); hold on; box on; grid on;
        set(gca,'LineWidth',lwid,'Fontsize',font);
        plot(ttimeA(implt),DMSP_facs(implt,iSC)*sign(mean(MLAT(implt,iSC))),'k','LineWidth',lwid);
        plot(ttimeA(implt),LTR_facs(implt,iSC)*sign(mean(MLAT(implt,iSC))),'color',cl_mage,'LineWidth',lwid); %,colors{iSC} % *sign(mean(MLAT(implt,iSC)))
        plot([1,1]*ttimeA(iredi),ylim,'color',cl_aubd,'Linestyle','--','LineWidth',lwidr);
        temp = real([DMSP_facs(implt,iSC);LTR_facs(implt,iSC)]*sign(mean(MLAT(implt,iSC))));
        ylimp = [min(temp),max(temp)]*1.05;
        ylim(ylimp); 
%         patch([xl,xl,xr,xr],[Ylim(1),Ylim(2),Ylim(2),Ylim(1)],'g','FaceAlpha',0.2);
        patch(tjdowns(icol,[1,1,2,2]),ylimp([1,2,2,1]),cl_fac,'FaceAlpha',0.2);
        plot(tjdowns(icol,:),[1,1]*ylimp(1)+lrlim*(ylimp(2)-ylimp(1)),'color',cl_fac,'LineWidth',4);
        set(gca,'xtick',ttick,'xticklabel',''); %datetick('x','HH:MM','keepticks');
        xlim(tlim); 
%         ylim(widlim1); 
        if icol==1
            ylabel(ylabFac,'Fontsize',font);
        end
        htt3=text(xtext,ylimp(1)+yrlim*(ylimp(2)-ylimp(1)),panlabs{3,icol},'Fontsize',font,'fontweight','bold');
        htt2b = text(min(tjdowns(icol,:)),ylimp(1)-trlim*(ylimp(2)-ylimp(1)),'FAC Down','color',cl_fac,'Fontsize',font,'fontweight','bold');
        
        hp4=subplot(nrow,ncol,ncol*3+icol); hold on; box on; grid on;
        set(gca,'LineWidth',lwid,'Fontsize',font);
        plot(ttimeA(implt),LTR_peds(implt,iSC),'color',cl_mage,'LineWidth',lwid); %,colors{iSC}
        plot([1,1]*ttimeA(iredi),ylim,'color',cl_aubd,'Linestyle','--','LineWidth',lwidr);
%         patch([xl,xl,xr,xr],[Ylim(1),Ylim(2),Ylim(2),Ylim(1)],'g','FaceAlpha',0.2);
        set(gca,'xtick',ttick,'xticklabel',tlabs,'Fontsize',font-1); 
%         datetick('x','HH:MM','keepticks');
        xlim(tlim); 
%         ylim(pedlim1); 
        temp = real(LTR_peds(implt,iSC));
        ylimp = [min(temp)*0.95,max(temp)*1.05];
        ylim(ylimp); 
        if icol==1
            ylabel(ylabPed,'Fontsize',font);
        end
        htt4=text(xtext,ylimp(1)+yrlim*(ylimp(2)-ylimp(1)),panlabs{4,icol},'Fontsize',font,'fontweight','bold');
%         xlabel('UT');
        if icol==1
            htt4a = text(tlim(1)-3/1440,ylimp(1)-0.285*(ylimp(2)-ylimp(1)),tlab0,'Fontsize',font);
%             htt4a = text(tlim(1)-3/1440,ylimp(1)-0.2*(ylimp(2)-ylimp(1)),tlab0,'Fontsize',font);
        end
        
        print(fig,pngname,'-dpng');
        print(fig,pngname,'-depsc');
%         close(fig);
    end