function f = TPMS()
f1 = @(x1,y1,z1) cos(x1) + cos(y1) + cos(z1) + 0.2 ; % P-element
f2 = @(x2,y2,z2) cos(x2).*cos(y2).*cos(z2) - sin(x2).*sin(y2).*sin(z2) + 0.2; % D-element
f3 = @(x3,y3,z3) sin(x3).*cos(y3) + sin(y3).*cos(z3) + sin(z3).*cos(x3);  % G-element
f4 = @(x4,y4,z4) 2.*(cos(x4).*cos(y4) + cos(y4).*cos(z4) + cos(z4).*cos(x4)) ...
                 - (cos(2.*x4) + cos(2.*y4) + cos(2.*z4));  % I-WP-element
f5 = @(x5,y5,z5) -2.*(cos(2.*x5).*cos(2.*y5) + cos(2.*y5).*cos(2.*z5) + cos(2.*z5).*cos(2.*x5))...
                 + 8.*cos(x5).*cos(y5).*cos(z5);  % F-RD-element
f6 = @(x6,y6,z6) cos(2.*pi.*x6 - pi) + cos(2.*pi.*y6 - pi) + cos(2.*pi.*z6 - pi);  % P-element(diff domain)
f7 = @(x7,y7,z7) 0.8 -(f1(x7,y7,z7) + 0.2);  % P-element-opposite
f8 = @(x8,y8,z8) sin(x8 + y8)./sqrt(x8.^2 + y8.^2 + z8.^2 + 1);  % Phi(x,y,z,c) = A(x,y,z) * RBF(x,y,z,c)
f9 = @(x9,y9,z9) f1(x9,y9,z9) - f8(x9,y9,z9);  % P-element(x,y,z,c) - Phi(x,y,z,c)
f10 = @(x10,y10,z10) f7(x10,y10,z10) + f8(x10,y10,z10);  % P-element-oppsite(x,y,z,c) + Phi(x,y,z,c)
Ralpha = 0;
f11 = @(x11,y11,z11) 1/(1+Ralpha).*(f9(x11,y11,z11) + f10(x11,y11,z11)...
        - power(power(f9(x11,y11,z11),2) + power(f10(x11,y11,z11),2)...
                - 2.* Ralpha.* f9(x11,y11,z11).* f10(x11,y11,z11),1/2));  % P-shell(with RBF)

[x,y,z] = meshgrid(-3.2:.1:3.2,-3.2:.1:3.2,-3.2:.1:3.2);
%[x,y,z] = meshgrid(-1.58:.1:4.72,-1.58:.1:4.72,-1.58:.1:4.72);
%[x,y,z] = meshgrid(0:.1:6.3,0:.1:6.3,0:.1:6.3);
%[x,y,z] = meshgrid(-5:.1:5,-5:.1:5,-5:.1:5);
%[x,y,z] = meshgrid(-10:.1:10,-10:.1:10,-10:.1:10);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% v = f8(x,y,z);
% h1 = patch(isosurface(x,y,z,v,0));
% %h2 = patch(isocaps(x,y,z,v,0));
% set(h1,'FaceColor','[0.5 0.5 0.5]','EdgeColor','none');
% %set(h2,'FaceColor','green','EdgeColor','none');
% xlabel('x');ylabel('y');zlabel('z'); 
% alpha(1); grid on; view([1,1,1]); axis equal; camlight; lighting gouraud

% hold on
% v = f7(x,y,z);
% h1 = patch(isosurface(x,y,z,v,0));
% %h2 = patch(isocaps(x,y,z,-v,0));
% set(h1,'FaceColor','red','EdgeColor','none');
% %set(h2,'FaceColor','blue','EdgeColor','none');
% xlabel('x');ylabel('y');zlabel('z'); 
% alpha(1); grid on; view([1,1,1]); axis equal; camlight; lighting gouraud
 
% P-Shell
v = f11(x,y,z);
h1 = patch(isosurface(x,y,z,v,0));
%h2 = patch(isocaps(x,y,z,v,0));
set(h1,'FaceColor','red','EdgeColor','none');
%set(h2,'FaceColor','blue','EdgeColor','none');
xlabel('x');ylabel('y');zlabel('z'); 
alpha(1); grid on; view([1,1,1]); axis equal; camlight; lighting gouraud

f = f11;