%% Figure 4 and Figure 6 of Lin21 JGR SAPS paper.
close all; clear;

addpath('~ldong/mcode/color_brewer/'); %colormap(brewermap([],'YlGnBu'))
dash_dir = '/glade/work/ldong/DASH/Lin20_DEP_SAPS/';
plt_dir = dash_dir;
% matname = [dash_dir,'fig3_dep_on_off.mat'];
% save(matname,'xi','yi','efxn1','facn1','widn1','UI_corot','efxn2','facn2','widn2','LATsmp','iMLT0');
ifig = 4;
if ifig==4
    ncfn = 'fig4_dep_on_off.nc'; % generated with SAPS_paper2_fig3.m with ifig=4;
    srunttl1 = 'Diffuse on'; 
    srunttl2 = 'Diffuse off';
    xr1 = -0.25;
    xr2 = -0.25;
else
    ncfn = 'fig6_fed_10.nc'; % generated with SAPS_paper2_fig3.m with ifig=5;
    srunttl1 = 'Fedder95';
    srunttl2 = {'RCM+Zhang15','1/10 loss rate'};
    xr1 = -0.25;
    xr2 = -0.4;
end
yr1 = 0.95;

    % Prepare for variables to be written to nc file.
    ncname = [dash_dir,ncfn]; %'fig4_dep_on_off.nc'
    ncid = netcdf.open(ncname,'NC_NOWRITE'); 
    
    varid = netcdf.inqVarID(ncid,'xi');   
    xi = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'yi');   
    yi = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'efxn1');   
    efxn1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'facn1');   
    facn1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'widn1');   
    widn1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'efxn2');   
    efxn2 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'facn2');   
    facn2 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'widn2');   
    widn2 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'LATsmp');   
    LATsmp = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'UI_corot');   
    UI_corot = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'TFACn1');   
    TFACn1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'TENGn1');   
    TENGn1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'TFACn2');   
    TFACn2 = netcdf.getVar(ncid,varid);
    varid= netcdf.inqVarID(ncid,'TENGn2');   
    TENGn2 = netcdf.getVar(ncid,varid);
    varid= netcdf.inqVarID(ncid,'iday1');   
    iday1 = netcdf.getVar(ncid,varid);
    strdate1 = datestr(iday1,'yyyy-mm-ddTHH-MM-ss');
    %% Plot
    MLTi = atan2(yi,xi)*180/pi/15+12;
    LATi = acos(sqrt(xi.^2+yi.^2))*180/pi;
    MLTsmp = 20;
    [~,iMLT0] = min(abs(MLTi(:,1)-MLTsmp));
    LATsmp = LATi(iMLT0,:);
    
    lwid=2; 
    font=16;
    efxlim = [0,20];
    faclim = [-1,1];
    widlim = [-1,1]*2;
    efxlim4 = [0,24];
    efxtick = 0:3:24;
    faclim4 = [-1,1]*2;
    widlim4 = [-1,1]*2;
    latlim = [45,90];
    efxlvl = [1,1]*1;
    faclvl = [1,1]*0.2;
    lats = 50:10:90;
    dpot=10;
    cl_efx = 'm'; % efx
    cl_fac = [0.4660, 0.6740, 0.1880]; % fac 
    cl_wid = 'b'; % wid
%     cl_lat1= [0.3010, 0.7450, 0.9330];
    cl_lat1= [0.8500, 0.3250, 0.0980];
    xp1 = 0.7;
    yp1 = 0.7;
    xp2 = 46;
    yp2 = 22;
    xtr = 0.70; 
    ytr = -0.35; 
    xtt2= 0.70; 
    ytt2= -0.35; 
    xbl = -0.65; % min lab bottom left
    ybl =  0.70; 
    xbr = -0.65; 
    ybr = -0.40; 
%     xlv1 = -0.5490;
%     ylv1 = -0.5249;
%     xlv2 = -0.3779;
%     ylv2 = 1.07; 
    xlv1 =-0.7011;
    ylv1 =-0.2884;
    xlv2 = 0.4448;
    ylv2 = 1.09;
    
    nrow=2; ncol=4; fpos = [1,39,1740,830];
    pngname = [plt_dir,ncfn(1:end-3)];
    disp(['Plotting ',pngname]);
    panlabs = {'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)'};
    panlabs = reshape(panlabs,nrow,ncol);
    fig = figure('position',fpos,'visible','on');
    
    hp1 = subplot(nrow,ncol,1);
    colormap(hp1,brewermap([],'*YlGnBu'));
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,efxn1); shading flat; c1=colorbar('southoutside'); caxis(efxlim);
    contour(xi,yi,efxn1,efxlvl,'color',cl_efx,'LineWidth',lwid*1.2,'LineStyle','-');
    plotion(cl_lat1,lats);
    ylabel(c1,'EnFlux [mW/m^2]');
    text(xr1,yr1,srunttl1,'fontsize',font,'fontweight','bold','rotation',90);
    hp1t2=text(xtr,    ytr,         0,[num2str(TENGn1,'%.1f'),' GW'],'fontsize',font-1); 
    hp1t3=text(xlv1,ylv1,[num2str(efxlvl(1),'%.1f'),' mW/m^2'],'color',cl_efx,'fontsize',font-1,'fontweight','bold');
%     set(c1,'position');
    htt1=text(xp1,yp1,panlabs{1},'Fontsize',font,'fontweight','bold');
    
    hp2 = subplot(nrow,ncol,2);
%     colormap(hp2,redblue);
    colormap(hp2,brewermap([],'*PuOr'));
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,facn1); shading flat; c2=colorbar('southoutside'); caxis(faclim);
    contour(xi,yi,facn1,faclvl,'color',cl_fac,'LineWidth',lwid*1.2,'LineStyle','-');
%     plot_pot(xi,yi,potn1,dpot,cl1,lwid/2);
    plotion('k',lats);
    ylabel(c2,'FAC [\muA/m^2]');
    hp2t2=text(xtt2,    ytt2,         0,[num2str(TFACn1,'%.1f'),' MA'],'fontsize',font-1);
    hp2t3=text(xbl, ybl,         0,{'min: ',[num2str(min(facn1(:)),'%.1f'),' \muA/m^2']},'fontsize',font-1);
    hp2t4=text(xbr, ybr,         0,{'max: ',[num2str(max(facn1(:)),'%.1f'),' \muA/m^2']},'fontsize',font-1);
    hp2t5=text(xlv2,ylv2,[num2str(faclvl(1),'%.1f'),' \muA/m^2'],'color',cl_fac,'fontsize',font-1,'fontweight','bold');
    ht2=title(strdate1);
    set(ht2,'position',[0.85, -0.8, 0]);
    htt2=text(xp1,yp1,panlabs{2},'Fontsize',font,'fontweight','bold');
    
    hp3 = subplot(nrow,ncol,3);
    colormap(hp3,brewermap([],'*RdBu'));%RdYlBu  PuOr
%     colormap(hp3,redblue);
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,widn1-UI_corot); shading flat; c3=colorbar('southoutside'); caxis(widlim);
    contour(xi,yi,efxn1,efxlvl,'color',cl_efx,'LineWidth',lwid*1.2,'LineStyle','-');
    contour(xi,yi,facn1,faclvl,'color',cl_fac,'LineWidth',lwid*1.2,'LineStyle','-');
%     plotion('k',lats);
    plotion('k',lats);
    ylabel(c3,'WID [km/s]');
    htt3=text(xp1,yp1,panlabs{3},'Fontsize',font,'fontweight','bold');
    
    hp4 = subplot(nrow,ncol,4); 
    yyaxis left;
    hold on; box on; grid on;
    set(gca,'YColor',cl_efx);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    hl41=plot(LATsmp,efxn1(iMLT0,:),'color',cl_efx,'LineWidth',lwid,'LineStyle','-');
    xlim(latlim);
    ylim(efxlim4);
    set(gca,'ytick',efxtick);
%     xlabel('Lat [deg]');
    ylabel('EnFlux [mW/m^2]');
    htt4=text(xp2,yp2,panlabs{4},'Fontsize',font,'fontweight','bold');
    %
    yyaxis right;
    axy1 = gca;
    hold on; box on; grid on;
    set(gca,'YColor',cl_fac);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    hl42=plot(LATsmp,facn1(iMLT0,:),'color',cl_fac,'LineWidth',lwid,'LineStyle','-');
    xlim(latlim);
    ylim(faclim4);
    ht41=ylabel('FAC [\mu A/m^2]');
    set(ht41,'position',[94,0,-1]);
    ht4 = title(['MLT=',num2str(MLTsmp,'%.1f')]);
    set(ht4,'position',[67.9487   -7.5736   -0.5000]);
    % add a third y axis to the right for WID.
    hold on;
    posy1 = get(gca,'position');
    dyr = 0.04;
    axy3 = axes('Position',posy1+[0,0,0.04,0],'box','off','Color','none','XColor','none','YColor',cl_wid,'xtick',[],...
        'xlim',[latlim(1),latlim(2)+11],'ylim',widlim4,'yaxislocation','right');
    set(axy3,'LineWidth',lwid,'Fontsize',font);
    hl43 = line(LATsmp,widn1(iMLT0,:)-UI_corot(iMLT0,:),'Color',cl_wid,'LineWidth',lwid,'LineStyle','-','Parent',axy3);
    ylabel(axy3,'WID [km/s]');
    
    hp5 = subplot(nrow,ncol,5);
    colormap(hp5,brewermap([],'*YlGnBu'));
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,efxn2); shading flat; %c5=colorbar; 
    caxis(efxlim);
    contour(xi,yi,efxn2,efxlvl,'color',cl_efx,'LineWidth',lwid*1.2,'LineStyle','-');
    plotion(cl_lat1,lats);
%     ylabel(c5,'EnFlux [mW/m^2]');
    text(xr2,yr1,srunttl2,'fontsize',font,'fontweight','bold','rotation',90);
    hp5t2=text(xtr,    ytr,         0,[num2str(TENGn2,'%.1f'),' GW'],'fontsize',font-1); 
    hp5t3=text(xlv1,ylv1,[num2str(efxlvl(1),'%.1f'),' mW/m^2'],'color',cl_efx,'fontsize',font-1,'fontweight','bold');
    htt5=text(xp1,yp1,panlabs{5},'Fontsize',font,'fontweight','bold');
    
    hp6 = subplot(nrow,ncol,6);
%     colormap(hp6,redblue);
    colormap(hp6,brewermap([],'*PuOr'));
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,facn2); shading flat; %c6=colorbar; 
    caxis(faclim);
    contour(xi,yi,facn2,faclvl,'color',cl_fac,'LineWidth',lwid*1.2,'LineStyle','-');
    plotion('k',lats);
%     ylabel(c6,'FAC [\mu A/m^2]');
    hp6t2=text(xtt2,    ytt2,         0,[num2str(TFACn2,'%.1f'),' MA'],'fontsize',font-1); 
    hp6t3=text(xbl, ybl,         0,{'min: ',[num2str(min(facn2(:)),'%.1f'),' \muA/m^2']},'fontsize',font-1);
    hp6t4=text(xbr, ybr,         0,{'max: ',[num2str(max(facn2(:)),'%.1f'),' \muA/m^2']},'fontsize',font-1);
    hp6t5=text(xlv2,ylv2,[num2str(faclvl(1),'%.1f'),' \muA/m^2'],'color',cl_fac,'fontsize',font-1,'fontweight','bold');
    htt6=text(xp1,yp1,panlabs{6},'Fontsize',font,'fontweight','bold');
    
    hp7 = subplot(nrow,ncol,7);
%     colormap(hp7,redblue);
    colormap(hp7,brewermap([],'*RdBu'));%RdYlBu PuOr
    hold on; axis equal tight off; view(-90,90);
    set(gca,'LineWidth',lwid,'Fontsize',font);
    pcolor(xi,yi,widn2-UI_corot); shading flat; %c7=colorbar; 
    caxis(widlim);
    contour(xi,yi,efxn2,efxlvl,'color',cl_efx,'LineWidth',lwid*1.2,'LineStyle','-');
    contour(xi,yi,facn2,faclvl,'color',cl_fac,'LineWidth',lwid*1.2,'LineStyle','-');
    plotion('k',lats);
%     ylabel(c7,'WID [km/s]');
    htt7=text(xp1,yp1,panlabs{7},'Fontsize',font,'fontweight','bold');
    
    hp8 = subplot(nrow,ncol,8); 
    yyaxis left;
    hold on; box on; grid on;
    set(gca,'LineWidth',lwid,'Fontsize',font);
    hl81=plot(LATsmp,efxn2(iMLT0,:),'color',cl_efx,'LineWidth',lwid,'LineStyle','-');
    xlim(latlim);
    ylim(efxlim4);
    set(gca,'ytick',efxtick);
    ht82=xlabel('MLAT [deg]');
    set(ht82,'position',[67.5,-2.5,-1]);
    ylabel('EnFlux [mW/m^2]');
    set(gca,'YColor',cl_efx);
    htt8=text(xp2,yp2,panlabs{8},'Fontsize',font,'fontweight','bold');
    %
    yyaxis right;
    hold on; box on; grid on;
    set(gca,'LineWidth',lwid,'Fontsize',font);
    hl82=plot(LATsmp,facn2(iMLT0,:),'color',cl_fac,'LineWidth',lwid,'LineStyle','-');
%     hl83=plot(LATsmp,widn2(iMLT0,:)-UI_corot(iMLT0,:),'color',cl_wid,'LineWidth',lwid,'LineStyle','-');
    xlim(latlim);
    ylim(faclim4);
    ht81=ylabel('FAC [\mu A/m^2]');
    set(ht81,'position',[94,0,-1]);
    set(gca,'YColor',cl_fac);
    % add a third y axis to the right for WID.
    hold on;
    posy1 = get(gca,'position');
    dyr = 0.04;
    axy3 = axes('Position',posy1+[0,-0.03,0.04,0],'box','off','Color','none','XColor','none','YColor',cl_wid,'xtick',[],...
        'xlim',[latlim(1),latlim(2)+11],'ylim',widlim4,'yaxislocation','right');
    set(axy3,'LineWidth',lwid,'Fontsize',font);
    hl83 = line(LATsmp,widn2(iMLT0,:)-UI_corot(iMLT0,:),'Color',cl_wid,'LineWidth',lwid,'LineStyle','-','Parent',axy3);
    ylabel(axy3,'WID [km/s]');
    
    set(hp1,'position',[0.1100    0.5838    0.1566    0.3412]);
    set(hp2,'position',[0.3161    0.5838    0.1566    0.3412]);
    set(hp3,'position',[0.5222    0.5838    0.1566    0.3412]);
    set(c1 ,'position',[0.1100    0.5174    0.1566    0.0264]);
    set(c2 ,'position',[0.3161    0.5174    0.1566    0.0264]);
    set(c3 ,'position',[0.5222    0.5174    0.1566    0.0264]);
    
    set(hp5,'position',[0.1100    0.0800    0.1566    0.3412]);
    set(hp6,'position',[0.3161    0.0800    0.1566    0.3412]);
    set(hp7,'position',[0.5222    0.0800    0.1566    0.3412]);
    set(hp8,'position',[0.7484    0.0800    0.1566    0.3412]);
%     set(c5 ,'position',[0.2663    0.1100    0.0111    0.3412]);
%     set(c6 ,'position',[0.4683    0.1100    0.0111    0.3412]);
%     set(c7 ,'position',[0.6744    0.1100    0.0111    0.3412]);
%     legend(hl82,'FAC [\mu A/m^2]');
    set(gcf,'renderer','painters');
    print(fig,pngname,'-dpng');
    print(fig,pngname,'-depsc');
%     close(fig);