function VOXEL = ETPMSVoxel(num_size,a,b,c)
h = 1.0/num_size;
xx = h/2:h:1-h/2;
[x,y,z]=meshgrid(xx,xx,xx);
Phi = cos(2*pi*x-pi)+a.*cos(2*pi*y-pi)+b.*cos(2*pi*z-pi)-c;
% VOXEL = P;
% VOXEL(P>0) = 1;
% VOXEL(P<0) = 1e-9;

k=1e6;
% H = 1./(1+exp(-k.*Phi)).*(Phi>=0)+exp(k.*Phi)./(1+exp(k.*Phi)).*(Phi<0);
H = 1./(1+exp(-k.*Phi));
VOXEL = max(H,1e-9);

% alpha = 0.005;
% eps = 1e-9;
% VOXEL=3*(1-eps)/4*(Phi./alpha-Phi.^3./(3.*alpha.^3))+(1+eps)/2;
% VOXEL(Phi>alpha) = 1;
% VOXEL(Phi<-alpha) = eps;
end