1 function [CX,CH]=center_mesh(X)
3 % compute centers of mesh cells by averaging the coordinates of the cell corners
5 % X {x,y,z} 3D arrays of coordinates of mesh corners
7 % CX {x,y,z} 3D arrays of coordinates of mesh cell centers
8 % CW {x,y,z} 3D arrays of coordinates of mesh horizontal face centers
13 CX{k} = (X{k}(1:end-1,1:end-1,1:end-1)+X{k}(2:end,1:end-1,1:end-1)+...
14 X{k}(1:end-1,2:end,1:end-1)+X{k}(1:end-1,1:end-1,2:end)+...
15 X{k}(2:end,2:end,1:end-1)+X{k}(2:end,1:end-1,2:end)+...
16 X{k}(1:end-1,2:end,2:end)+X{k}(2:end,2:end,2:end))/8;
19 CB = X{3}(1:end-1,1:end-1,1)+X{3}(2:end,1:end-1,1)+...
20 X{3}(1:end-1,2:end,1)+X{3}(1:end-1,1:end-1,1);
21 % center height above terrain
22 CH = zeros(size(CX{3}));
24 CH(:,:,k)=CX{3}(:,:,k)-CB;