%% 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;
% matname = [dash_dir,'fig4_saps_3d.mat'];
% save(matname,'surf_efxIB_x','surf_efxIB_y','surf_efxIB_z','con12_mix','con12_rcm','iefxc2','ifacc2','iefxc_3d','ifacc_3d',...
%     'xi_glb','yi_glb','zi_glb','Rinset','WIDn_glb','surf_facIB_x','surf_facIB_y','surf_facIB_z','iefxc1','ifacc1','rho_eq_plt','xe','ye','ze','fltEq2Ni','fltraceEq2N');

% xi, yi, WIDn, % 361x45
% xe, ye, rho_eq, % 96x193
% surf_efxIB_x,surf_efxIB_y,surf_efxIB_z, % 100x296
% surf_facIB_x,surf_facIB_y,surf_facIB_z, % 100x173
% iaccN2S, % 1817x1
% Ncc1, % 780
% con12_rcm, con12_mix % 1826x3

    % Prepare for variables to be written to nc file.
    ncname = [dash_dir,'fig5_saps_3d.nc']; % generated with SAPS_ReMIX_RCM_plt_3Din.m, called by SAPS_ReMIX_RCM_plt_SAPS2.m, called by SAPS_ReMIX_RCM.m
    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,'WIDn');   
    WIDn = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'xe');   
    xe = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'ye');   
    ye = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'ze');   
    ze = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'rho_eq');   
    rho_eq = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_efxIB_x');   
    surf_efxIB_x = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_efxIB_y');   
    surf_efxIB_y = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_efxIB_z');   
    surf_efxIB_z = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_facIB_x');   
    surf_facIB_x = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_facIB_y');   
    surf_facIB_y = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'surf_facIB_z');   
    surf_facIB_z = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'iaccN2S');   
    iaccN2S = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'Ncc1');   
    Ncc1 = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'con12_rcm');   
    con12_rcm = netcdf.getVar(ncid,varid);
    varid = netcdf.inqVarID(ncid,'con12_mix');   
    con12_mix = netcdf.getVar(ncid,varid);
    %% Plot
    ifacLLB = (644:823)+Ncc1; %128:570;
    ifacHLB = (823:999)+Ncc1; %128:570;
    iefxLLB = (35:331);%(735:1133)+Ncc1;
    iefxHLB = (331:692);%(735:1133)+Ncc1;

    iefxc_3d = find(iaccN2S<=iefxHLB(end)&iaccN2S>=iefxLLB(1)); 
    iefxc_3dL= find(iaccN2S<=iefxLLB(end)&iaccN2S>=iefxLLB(1)); 
    iefxc_3dH= find(iaccN2S<=iefxHLB(end)&iaccN2S>=iefxHLB(1)); 
    ifacc_3d = find(iaccN2S<=ifacHLB(end)&iaccN2S>=ifacLLB(1));
    ifacc_3dL= find(iaccN2S<=ifacLLB(end)&iaccN2S>=ifacLLB(1));
    ifacc_3dH= find(iaccN2S<=ifacHLB(end)&iaccN2S>=ifacHLB(1));
    
    %% Expand variables for remix plot
    Rin = 1.02; % radius of ISP surface.
    [nLon,nLat] = size(xi);
    zi = sqrt(1-xi.^2-yi.^2);
    mix_colat = acos(zi)*180/pi; mix_lon = atan2(yi,xi)*180/pi;
    mix_colat_glb = zeros(nLon,nLat*2); mix_lon_glb = mix_colat_glb;
    mix_colat_glb(:,1:nLat) = mix_colat; 
    mix_colat_glb(:,nLat+1:end) = repmat(45.6+(1:nLat),nLon,1);
    mix_lon_glb(:,1:nLat) = mix_lon; 
    mix_lon_glb(:,nLat+1:end) = repmat(mix_lon(:,1),1,nLat);
    xi_glb = sind(mix_colat_glb).*cosd(mix_lon_glb)*Rin;
    yi_glb = sind(mix_colat_glb).*sind(mix_lon_glb)*Rin;
    zi_glb = cosd(mix_colat_glb)*Rin;
    WIDn_glb = zeros(size(xi_glb));
    WIDn_glb(:,1:nLat) = WIDn;
    WIDn_glb(:,nLat-3:nLat) = 0;
    rho_eq_plt = log10(rho_eq);
    %
    pngname = [plt_dir,'fig5_saps_3d'];
    disp(['Plotting ',pngname]);
    nrow=1; ncol=1; fpos=[500,39,936,800];
    zfac=1;
    widlim = [-1,1]*2; %Plim=[0,50];
    Xeqlim=[-15,15]; 
    Yeqlim=[-30,30]/2;  
    Zlim=[0,15*zfac];
    az = -200; 
    el = 40;
    cl_efx = 'm'; % efx
    cl_fac = [0.4660, 0.6740, 0.1880]; % fac 
    cl_wid = 'b'; % wid
%     cl_sun = [1.0,0.5,0];
    lwid = 2;
    xa0 = Xeqlim(1);
    ya0 = Yeqlim(1);
    za0 = 0;
    axl = 12;
    
    fig = figure('position',fpos); 
    % Plot MHD density distribution in eq. plane.
    ax0 = axes;
    hold on; axis equal tight off; box off; view(3);
    colormap(ax0,brewermap([],'*PuOr'));
    hs0=surf(xe,ye,ze*zfac,rho_eq_plt); shading interp; alpha(hs0,1.0);
    xlim(Xeqlim); 
    ylim(Yeqlim); 
    zlim(Zlim);
    view(az,el);
    
    % Plot FAC/EnFlux contour lines in eq. plane.
    ax1 = axes;
    hold on; axis equal tight off; box off; view(3);
    iefxc1 = [iefxc_3d; iefxc_3d(3)];
    ifacc1 = [ifacc_3d; ifacc_3d(8)];
    con12_rcm_efx = con12_rcm(iefxc1,:);
    con12_rcm_fac = con12_rcm(ifacc1,:);
    plot3(con12_rcm_efx(:,1),con12_rcm_efx(:,2),con12_rcm_efx(:,3)*0,'color',cl_efx,'LineStyle','-','LineWidth',lwid);
    plot3(con12_rcm_fac(:,1),con12_rcm_fac(:,2),con12_rcm_fac(:,3)*0,'color',cl_fac,'LineStyle','-','LineWidth',lwid);
    ht1=text(-6,6,0,{'Electron','precipitation','boundary'},'color',cl_efx,'fontweight','bold');
    set(ht1,'position',[8,15,5]);
    ht2=text(-7,7,0,{'Duskside','downward','FAC boundary'},'color',cl_fac,'fontweight','bold');
%     set(ht2,'position',[-7.2287,    7.3234,   -0.2688]);
    set(ht2,'position',[10.5,    2,   -0.2688]);
    xlim(Xeqlim); 
    ylim(Yeqlim); 
    zlim(Zlim);
    view(az,el);
    
    % Plot WID in 3D ISp.
    ax2 = axes;
    hold on; axis equal tight off; box off; view(3);
    colormap(ax2,brewermap([],'*RdBu'));
    [xei,yei]=meshgrid(-1:0.5:1,-1:0.5:1);
    zei = xei*0;
    hs0a=surf(xei,yei,zei,zei+1); shading interp; alpha(hs0,1.0);
    hsf2=surf(ax2,xi_glb,yi_glb,zi_glb,WIDn_glb); shading flat; caxis(widlim); hold on;
    iefxc2 = iefxc_3d(2:end);
    ifacc2 = ifacc_3d(8:end);
    plot3(con12_mix(iefxc2,1)*Rin,con12_mix(iefxc2,2)*Rin,(con12_mix(iefxc2,3))*Rin,'color',cl_efx,'LineStyle','-','LineWidth',lwid/3);
    plot3(con12_mix(ifacc2,1)*Rin,con12_mix(ifacc2,2)*Rin,(con12_mix(ifacc2,3))*Rin,'color',cl_fac,'LineStyle','-','LineWidth',lwid/3);
    if 1
    quiver3(1,0,0,axl,0,0,'k','LineWidth',lwid);
    quiver3(0,1,0,0,axl+2,0,'k','LineWidth',lwid);
    quiver3(0,0,1,0,0,axl+2.5,'k','LineWidth',lwid);
    text(axl+1,0,0,'X');
    text(0,axl+2,0,'Y');
    text(-1,0,axl+2,'Z');
    xyzbox0 = [-1,1,0; -1,-1,0; 1,-1,0; 1,1,0; -1,1,0];
    plot3(xyzbox0(:,1),xyzbox0(:,2),xyzbox0(:,3),'k','LineWidth',lwid/2);
    else
    quiver3(xa0+zeros(1,3),ya0+zeros(1,3),za0+zeros(1,3),[2,0,0],[0,2,0],[0,0,2],'k','LineWidth',lwid);
    text(xa0+2,ya0,0,'X');
    text(xa0,ya0+2,0,'Y');
    text(xa0,ya0,za0+2,'Z');
    end
    xlim(Xeqlim); 
    ylim(Yeqlim); 
    zlim(Zlim);
    view(az,el);
    
    % Plot sampled field lines from eq. to ISP. Try to make it two surfaces.
    ax3 = axes;
    hold on; axis equal tight off; box off; view(3);
    colormap(ax3,brewermap([],'*YlGnBu')); caxis([0,2.5]);
    hs2=surf(surf_facIB_x,surf_facIB_y,surf_facIB_z*zfac,surf_facIB_z+1); shading interp; alpha(hs2,0.7);
    plot3(surf_facIB_x(end,:),surf_facIB_y(end,:),surf_facIB_z(end,:)*zfac,'color',cl_fac);
    xlim(Xeqlim); 
    ylim(Yeqlim); 
    zlim(Zlim);
    view(az,el); 
    
    % Plot inset of WID in ISP.
    ax2a = axes('Position',[0.55,0.5,0.3,0.3]);
    Rinset = 1;
    hold on; axis equal tight off; box on; view(3);
    colormap(ax2a,brewermap([],'*RdBu'));
    [xei,yei]=meshgrid(-1:0.5:1,-1:0.5:1);
    zei = xei*0;
    dz = 0.5;
    hs0a=surf(xei,yei,zei,zei+1); shading interp; alpha(hs0,1.0);
    hsf2a=surf(ax2a,xi_glb*Rinset,yi_glb*Rinset,zi_glb*Rinset-dz,WIDn_glb); shading flat; caxis(widlim); hold on;
    iefxc2 = iefxc_3d(2:end);
    ifacc2 = ifacc_3d(8:end);
    plot3(con12_mix(iefxc2,1)*Rin*Rinset,con12_mix(iefxc2,2)*Rin*Rinset,(con12_mix(iefxc2,3))*Rin*Rinset-dz,'color',cl_efx,'LineStyle','-','LineWidth',lwid/1.6);
    plot3(con12_mix(ifacc2,1)*Rin*Rinset,con12_mix(ifacc2,2)*Rin*Rinset,(con12_mix(ifacc2,3))*Rin*Rinset-dz,'color',cl_fac,'LineStyle','-','LineWidth',lwid/1.6);
    xyzInlim = [-1,1]*Rin^2;
    xlim(xyzInlim); 
    ylim(xyzInlim); 
    zlim([0,2]);
    view(az,el);
    xyzbox0 = [-1,1,0; -1,-1,0; 1,-1,0; 1,1,0; -1,1,0];
    plot3(xyzbox0(:,1),xyzbox0(:,2),xyzbox0(:,3),'k','LineWidth',lwid/2);
%     lightangle(gca,az,el);
%     h.FaceLighting = 'gouraud';
%     h.AmbientStrength = 0.3;
%     h.DiffuseStrength = 0.8;
%     h.SpecularStrength = 0.9;
%     h.SpecularExponent = 25;
%     h.BackFaceLighting = 'lit';
    
    % Plot sampled field lines from eq. to ISP. Try to make it two surfaces.
    ax4 = axes;
    hold on; axis equal tight off; box on; view(3);
    colormap(ax4,redblue);
    hs1=surf(surf_efxIB_x,surf_efxIB_y,surf_efxIB_z*zfac); shading interp; alpha(hs1,0.5);
    plot3(surf_efxIB_x(end,:),surf_efxIB_y(end,:),surf_efxIB_z(end,:)*zfac,cl_efx);
    xlim(Xeqlim); 
    ylim(Yeqlim); 
    zlim(Zlim);
    view(az,el);
    
    % Add labels for SAPS.
    plot3([ 1,-7],[ 1,1.8],[0,6],'k','LineWidth',lwid/2);
    plot3([ 1,-7.6],[-1,-7.7],[0,5],'k','LineWidth',lwid/2);
    xsa = -10;
    ysa = 6;
    zsa = 2;
    ht4=text(xsa,ysa,zsa,'SAPS','color',cl_wid,'fontweight','bold','fontsize',10);
%     quiver3(xsa,ysa,zsa,-5*0.94,-13*0.94,3.8,'color',cl_wid,'LineWidth',1,'Markersize',8);
%     % for not subtract 0.5 when plotting inset.
    quiver3(xsa,ysa,zsa,-3.2,-9.5,3.8,'color',cl_wid,'LineWidth',1,'Markersize',8);
    quiver3(xsa,ysa,zsa,11,0,0,'color',cl_wid,'LineWidth',1,'Markersize',8);
    
    print(fig,pngname,'-dpng');
    print(fig,pngname,'-depsc');
    
    %% END of CODE