dx=[0.5*ones(1,120) ];
dy=[0.5*ones(1,120) ];
dz=[2*ones(1,1)];
nx=size(dx,2);
ny=size(dy,2);
nz=size(dz,2);
x_size = 60;
y_size = 60;
coordinates_pre = zeros(nx*ny*nz,2);
select_t = 120;
for k = 1:nz
    for j = 1:ny
        for i = 1:nx
            if i == 1
                coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,1) = sum(dx(1:1),2)/2;
            else
                coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,1) = (sum(dx(1:i-1),2)+sum(dx(1:i),2))/2;
            end
            if j == 1
                coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,2) = sum(dy(1:1),2)/2;
            else
                coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,2) = (sum(dy(1:j-1),2)+sum(dy(1:j),2))/2;
            end
            %             if k == 1
            %                 coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,3) = sum(dz(1:1),2)/2;
            %             else
            %             coordinates_pre(i+(j-1)*nx+(k-1)*nx*ny,3) = (sum(dz(1:k-1),2)+sum(dz(1:k),2))/2;
            %             end
        end
    end
end
% select_t = 508;
zData = OutputRs{select_t, 1}.p(1:nx*ny*nz,1);
zData1 = 1-OutputRs{select_t, 1}.sw(1:nx*ny*nz,1);
xData = coordinates_pre(:,1);
yData = coordinates_pre(:,2);
% zData = DATA_plot(:,21);
% zData1 = DATA_plot(:,22);

% minx = min(xData);
% maxx = max(xData);
% miny = min(yData);
% maxy = max(yData);
minx = 0;
maxx = x_size;
miny = 0;
maxy = y_size;

tx = (minx:0.05:maxx)'; 
ty = (miny:0.05:maxy)'; 

[XI,YI] = meshgrid(tx,ty);
% 非画图区域标记
% 计算域多边形刻画
% Comput_domain_x = [0 40 200 160 200 40 0];
% Comput_domain_y = [40 80 80 40 0 0 40];
Comput_domain_x = [0 x_size x_size 0 0];
Comput_domain_y = [0 0 y_size y_size 0];
% load('coordinates_boundary_node_order.mat');
% coordinates_boundary_node_order(:,1) = coordinates_boundary_node_order(:,1) + 200;
% coordinates_boundary_node_order(:,2) = coordinates_boundary_node_order(:,2) + 130;
% Comput_domain_x = coordinates_boundary_node_order(:,1);
% Comput_domain_y = coordinates_boundary_node_order(:,2);
%
NX = size(XI,1);
NY = size(YI,1);
% for i = 1:NX
%     for j = 1:NY
%         in=inpolygon(XI,YI,Comput_domain_x,Comput_domain_x);
%     end
% end
in=inpolygon(XI,YI,Comput_domain_x',Comput_domain_y');
XI(~in) = nan; YI(~in) = nan;

% for i = 1:NX
%     for j = 1:NY
%         in=inpolygon(XI(i,j),YI(i,j),Comput_domain_x',Comput_domain_y');
%         if ~in    XI(i,j) = nan; YI(i,j) = nan; end
%     end 
% end
ZI = griddata(xData,yData,zData,XI,YI);
ZII = griddata(xData,yData,zData1,XI,YI);
figure('color','w');surf(XI,YI,ZI), hold
% plot3(xData,yData,zData,'o'), hold off
grid off
set(gca,'FontName','Times New Roman','FontSize',10,'FontWeight','bold');
xlabel( 'x, m' );
ylabel( 'y, m' );
zlabel( 'z, m' );
% grid on
view( 0, 90 );
axis equal;
xlim([0,60]);
ylim([0,60]);
colorbar
colormap jet
% caxis([10,15]);
hold off
shading interp
% set('edgecolor','none')
figure('color','w');surf(XI,YI,ZII), hold
% plot3(xData,yData,zData1,'o'), hold off
set(gca,'FontName','Times New Roman','FontSize',10,'FontWeight','bold');
xlabel( 'x, m' );
ylabel( 'y, m' );
zlabel( 'z, m' );
view( 0, 90 );
axis equal;
xlim([0,60]);
ylim([0,60]);
grid off
colorbar
colormap jet
caxis([0.2,0.8]);
hold off
shading interp