function Ve=ETPMS_volume(resolution,a,b,c)
nGrid=10;
h=1.0/(resolution*nGrid);
X = h/2:h:1-h/2;
[x,y,z]=meshgrid(X);
Phi = cos(2*pi*x-pi)+a.*cos(2*pi*y-pi)+b.*cos(2*pi*z-pi)-c;
k=1e6;
H = 1./(1+exp(-k.*Phi));

Ve = zeros(resolution,resolution,resolution);
for i = 1:resolution
    for j = 1:resolution
        for k = 1:resolution
%             dVe(i,j,k) = sum(sum(sum(f((i-1)*nGrid+1:i*nGrid,(j-1)*nGrid+1:j*nGrid,(k-1)*nGrid+1:k*nGrid))));
            Ve(j,i,k) = sum(sum(sum(H((i-1)*nGrid+1:i*nGrid,(j-1)*nGrid+1:j*nGrid,(k-1)*nGrid+1:k*nGrid))));
        end
    end
end
% Ve = Ve./(resolution*nGrid)^3;
Ve = Ve./(nGrid)^3;
end