%% Read first variable
clear all

filename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_updatemech_updateinflowU/data_changeUin_2.380E-02';
prodfilename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_updatemech_updateinflowU/prodRates_changeUin_2.380E-02';
% filename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_noTD_restartNicholas/data_noTD_2.720E-02';
% prodfilename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_noTD_restartNicholas/prodRates_noTD_2.720E-02';
% filename_prod = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion/cori_data_updatemech/prodRates_data_1.200E-02';
% filename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_noTD_restartNicholas/data_noTD_2.000E-02';
% prodfilename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_noTD_restartNicholas/prodRates_noTD_2.000E-02';
% filename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_updatemech_updateinflowU_noturb/data_noTD_3.120E-01';
% prodfilename = '/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_updatemech_updateinflowU_noturb/prodRates_3.120E-01';
T = NGAdatareader_large(filename,16);
% H = NGAdatareader_large(filename,8);
prod = NGAdatareader_large(prodfilename,13);
% peakprod = 168.1;
% prod = NGAdatareader_large(filename_prod,13);
[x,y,z,xm,ym,zm,nx,ny,nz,dx,dy,dz] = NGA_grid_reader('/media/jschlup/Storage/Research/H2/3D_Detailed/Papers/thermal_diffusion_updatemech_updateinflowU_noturb/config');
L = y(end);

oneD_withTD = importdata('/media/jschlup/Storage/Research/H2/1D_Detailed/WSSCI_Meeting/2015/followup/trans_nobindiff_output/MA_with_TD_UpdatedMech/H2_p01_0phi0_4000tu0298.kg');
[peaksource_withTD, r] = max(oneD_withTD.data(:,32));
peaksource_withTD = peaksource_withTD./oneD_withTD.data(r,3); % Divide by density to get units to match NGA

%% Plots
close all
colormap jet
[xg,yg,zg] = meshgrid(x(1:end-1),y(1:end-1),z(1:end-1));

for i=1:nx
    xm(i) = (x(i+1)+x(i))/2;
end
for j=1:ny
    ym(j) = (y(j+1)+y(j))/2;
end
for k=1:nz
    zm(k) = (z(k+1)+z(k))/2;
end

% figure()
% xslice = 0;
% yslice = [0,y(end-1)/L];
% zslice = [0,z(end-1)/L];
% slice(xg/L,yg/L,zg/L,permute(T,[2,1,3]),xslice,yslice,zslice);
% xlabel('x/L');
% ylabel('y/L');
% zlabel('z/L');
% set(gca,'FontSize',42);
% axis equal
% shading interp
% colormap jet
% % h = colorbar; title(h, '$T [K]$','Interpreter','Latex')
% % caxis([0 6])
% box on

% figure()
% contourf(xg(:,100:900,1)'/L,yg(:,100:900,1)'/L,H(100:900,:,120)*10^4,50,'EdgeColor','none') ; colormap jet; axis image;
% set(gca,'XTick',[]); set(gca,'YTick',[]);
% ax=gca;
% pos=get(gca,'pos');
% set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*0.95]);
% pos=get(gca,'pos');
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.15*pos(4) pos(3) 0.03]); title(h, '$Y_{H}\times 10^{4}$','Interpreter','Latex')
% set(gca,'FontSize',16);
% axis equal
% shading interp

% figure()
% contourf(xg(:,100:900,1)'/L,yg(:,100:900,1)'/L,T(100:900,:,120),50,'EdgeColor','none') ; colormap jet; axis image;
% set(gca,'XTick',[]); set(gca,'YTick',[]);
% ax=gca;
% pos=get(gca,'pos');
% set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*0.95]);
% pos=get(gca,'pos');
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.15*pos(4) pos(3) 0.03]); title(h, '$T [K]$','Interpreter','Latex')
% set(gca,'FontSize',16);
% axis equal
% shading interp

% figure()
% contourf(xg(:,100:900,1)'/L,yg(:,100:900,1)'/L,prod(100:900,:,120)/peaksource_withTD,50,'EdgeColor','none') ; colormap jet; axis image;
% set(gca,'XTick',[]); set(gca,'YTick',[]);
% ax=gca;
% pos=get(gca,'pos');
% set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*0.95]);
% pos=get(gca,'pos');
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.15*pos(4) pos(3) 0.03]); title(h, '$\dot{\omega}_{H_2O}^{TD}/\dot{\omega}_{H_2O}^{1D}$','Interpreter','Latex')
% set(gca,'FontSize',16);
% axis equal
% shading interp
%
close all

startval = 1520/4; endval = 3*1520/4;
% startval = 1; endval = 1520;

figure('units','normalized','outerposition',[0 0 1 1])
contourf(xg(:,startval:endval,1)'/L,yg(:,startval:endval,1)'/L,T(startval:endval,:,120),50,'EdgeColor','none') ; colormap jet; axis image;
ax=gca;
set(ax,'XTick',[]); set(ax,'YTick',[]);
pos=get(gca,'pos');
set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*1.00]);
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.3*pos(4) pos(3) 0.03]); xlabel(h, '$T$ [K]','Interpreter','Latex')
% set(gca,'FontSize',36);
shading interp
export_fig('figures/3D_contourT_full.png');

f=figure('units','normalized','outerposition',[0 0 1 1]);
% subplot(2,1,1)
contourf(xg(:,startval:endval,1)'/L,yg(:,startval:endval,1)'/L,T(startval:endval,:,120),50,'EdgeColor','none') ; colormap jet; axis image;
pos=get(gca,'pos');
set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*1.00]);
% xlim([0 8]); ylim([0 1+100*eps]);
caxis([300 1600])
h = colorbar; h.Label.String = '$T [K]$'; h.Label.Interpreter = 'Latex'; h.TickLabelInterpreter = 'Latex'; h.Ticks = [500 1000 1500]; h.LineWidth = 2; %title(h, '$T$ [K]','Interpreter','Latex')
ax=gca;
set(ax,'XTick',[]); set(ax,'YTick',[]);
set(gca,'FontSize',40);
c = get(f, 'children');
set(c(2), 'visible', 'off');
f.CurrentAxes.Children.Visible = 'off';
h.Position = [pos(4) pos(2) pos(3)*.02 pos(4)*0.4];
export_fig('figures/3D_contourT_colorbar.png');

% subplot(2,1,2)
figure('units','normalized','outerposition',[0 0 1 1])
contourf(xg(:,startval:endval,1)'/L,yg(:,startval:endval,1)'/L,prod(startval:endval,:,120)./peaksource_withTD,50,'EdgeColor','none') ; colormap jet; axis image;
ax=gca; caxis([0 15])
set(ax,'XTick',[]); set(ax,'YTick',[]);
pos=get(gca,'pos');
set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*1.00]);
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.3*pos(4) pos(3) 0.03]); xlabel(h, '$T$ [K]','Interpreter','Latex')
% set(gca,'FontSize',36);
shading interp
export_fig('figures/3D_contoursrc_full.png');

% subplot(2,1,2)
f=figure('units','normalized','outerposition',[0 0 1 1]);
contourf(xg(:,startval:endval,1)'/L,yg(:,startval:endval,1)'/L,prod(startval:endval,:,120)./peaksource_withTD,50,'EdgeColor','none') ; colormap jet; axis image;
pos=get(gca,'pos');
set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*1.00]);
caxis([0 16]);
h = colorbar; h.Label.String = '$\dot{\omega}_{\rm H_2O}/\dot{\omega}_{\rm H_2O}^{\rm 1D}$'; h.Label.Interpreter = 'Latex'; h.TickLabelInterpreter = 'Latex'; h.Ticks = [5 10 15]; h.LineWidth = 2; %title(h, '$\dot{\omega}_{H_2O}/\dot{\omega}_{H_2O}^{1D}$','Interpreter','Latex')
ax=gca; caxis([0 15])
set(ax,'XTick',[]); set(ax,'YTick',[]);
set(gca,'FontSize',40);
c = get(f, 'children');
set(c(2), 'visible', 'off');
f.CurrentAxes.Children.Visible = 'off';
h.Position = [pos(4) pos(2) pos(3)*.02 pos(4)*0.4];
export_fig('figures/3D_contoursrc_colorbar.png');

% f = figure('units','normalized','outerposition',[0 0 1 1]);
% contourf(xm/lf,ym/lf,cart_noTD_prodrates{7}{7}'./peaksource_noTD,50,'EdgeColor','none') ; colormap jet; axis image;
% caxis([0 3.5]);
% % pos=get(gca,'pos');
% h = colorbar; h.Label.String = '$\dot{\omega}_{H_2O}/\dot{\omega}_{H_2O}^{1D}$'; h.Label.Interpreter = 'Latex'; h.TickLabelInterpreter = 'Latex'; h.Ticks = [1 2 3]; h.LineWidth = 2; %title(h, '$\dot{\omega}_{H_2O}/\dot{\omega}_{H_2O}^{1D}$','Interpreter','Latex')
% set(gca,'pos',[pos(1) pos(2) pos(3)*0.8 pos(4)*0.8]);
% xlabel('$x/l_F$'); ylabel('$y/l_F$');
% xlim([0 xmax]); ylim([0 ymax]);
% ax = gca;
% set(ax,'XTick',[0 20 40 60 80 100]); set(ax,'YTick',[0 10 20]);
% set(ax,'FontSize',36,'LineWidth',2); set(ax,'XMinorTick','on'); ax.XAxis.MinorTickValuesMode = 'manual'; ax.XAxis.MinorTickValues = 0:5:120;
% c = get(f, 'children');
% set(c(2), 'visible', 'off');
% f.CurrentAxes.Children.Visible = 'off';
% export_fig('figures/2DFP_contoursrc_colorbar.png');

%
% figure()
% contourf(xg(:,:,1)'/L,yg(:,:,1)'/L,prod(:,:,120)/peaksource_withTD,50,'EdgeColor','none') ; colormap jet; axis image;
% set(gca,'XTick',[]); set(gca,'YTick',[]);
% ax=gca;
% pos=get(gca,'pos');
% set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)*0.95]);
% pos=get(gca,'pos');
% h =colorbar('location','southoutside','position',[pos(1) pos(2)+0.3*pos(4) pos(3) 0.03]); title(h, '$\dot{\omega}_{H_2O}^{TD}/\dot{\omega}_{H_2O}^{1D}$','Interpreter','Latex')
% set(gca,'FontSize',12);
% axis equal
% shading interp

figure()
p = patch(isosurface(xm/L,ym/L,zm/L,permute(T,[2,1,3]),1180));
isonormals(xm/L,ym/L,zm/L,permute(T,[2,1,3]),p)
isocolors(xm/L,ym/L,zm/L,permute(prod/peaksource_withTD,[2,1,3]),p)
p.FaceColor = 'interp';
p.EdgeColor = 'none';
view([-60 15]); axis tight
camlight left
lighting gouraud
colormap jet; box on; axis image; caxis([0 15]);
ax = gca; ax.Visible = 'off'; ax.FontSize = 24;
export_fig('figures/3D_isoT_laminar.png')