nx = 256; ny = 256; nz = 256; ratio = 2;
mx = nx/ratio; my = ny/ratio; mz = nz/ratio;
x = [0.5:ratio:nx-0.5]; y = [0.5:ratio:ny-0.5]; z = [0.5:ratio:nz-0.5];
xx = [0.5:1:nx-0.5]; yy = [0.5:1:ny-0.5];
rhol = 5.116; rhog = 2.027; midrho = 0.5*(rhol + rhog);
[X,Y,Z] = meshgrid(y,x,z); [XX,YY] = meshgrid(yy,xx);
filename1a = 'density3D00.dat';
filename1b = 'velocity3D00.dat';
filename2a = 'densityCut00.dat';
filename2b = 'velocityCut00.dat';
density3D = reshape(density3D,mx,my,mz);
velocity3D = reshape(velocity3D,mx,my,mz);
density2D = reshape(density2D,nx,ny);
velocity2D = reshape(velocity2D,nx,ny);
isosurface(X,Y,Z,density3D,midrho,velocity3D);
colormap(othercolor('BuDRd_12'));
c.Location = 'southoutside';
originalColorLimits = caxis(gca);
contour(YY,XX,density2D,midrho,'edgecolor','black','Linewidth',1.0,'Parent',HG);
contourf(YY,XX,velocity2D,30,'edgecolor','none','Parent',HG);
HG.Matrix = makehgtform('xrotate',pi/2);
caxis(gca, originalColorLimits);
axis([0 256 0 256 0 256]);
xlabel('$y$','fontname','Times','Fontsize',20,'interpreter','Latex');
ylabel('$x$','fontname','Times','Fontsize',20,'interpreter','Latex');
zlabel('$z$','fontname','Times','Fontsize',20,'interpreter','Latex');
set(gca,'XTick',aaa,'fontname','Times','fontsize',20)
set(gca,'YTick',bbb,'fontname','Times','fontsize',20)
set(gcf, 'Position', [100, 100, 800, 800])