In this post I want to present an easy method to obtain a discrete parametrization of a surface S defined implicitly (f(x,y,z)=0).

This problem was discussed here several times, the most recent post is

http://www.mapleprimes.com/posts/207661-Isolation-Of-Sides-Of-The-Surface-On-The-Graph

S is supposed to be the boundary of a convex body having (x0,y0,z0) an interior point and contained in a ball of radius R centered at (x0,y0,z0).

Actually, the procedure also works if the body is only star-shaped with respect to the interior point, and it is also possible to plot only a part of the surface

inside a solid angle centered at (x0,y0,z0).

*Usage:*

Par3d(f, x=x0, y=y0, z=z0, R, m, n, theta1 .. theta2, phi1 .. phi2)

f is an expression depending on the variables x, y, z

x0, y0, z0 are the coordinates of the interior point

R is the radius of the ball which contains the surface,

m, n are the numbers of the grid lines which will be generated

The last two parameters are optional and are used when only a part of S will be parametrized.

The procedure **Par3d** returns a **MESH** structure M, which can be plotted with **PLOT3D(M)**.

Par3d :=proc(f,x::`=`,y::`=`,z::`=`,R,m,n,th:=0..2*Pi,ph:=0..Pi)
local A,i,j, rij,fij,Cth,Sth,Cph,Sph, theta,phi, r;
A:=Array(1..m+1,1..n+1,1..3,datatype=float[8]);
for i from 0 to m do for j from 0 to n do
theta:=op(1,th)+i/m*(op(2,th)-op(1,th));
phi:=op(1,ph)+j/n*(op(2,ph)-op(1,ph));
Cth:=evalf(cos(theta)); Sth:=evalf(sin(theta));
Cph:=evalf(cos(phi)); Sph:=evalf(sin(phi));
fij:= eval(f, [lhs(x)=rhs(x)+r*Sph*Cth, lhs(y)=rhs(y)+r*Sph*Sth, lhs(z)=rhs(z)+r*Cph]);
rij:=fsolve(fij,r=0..R);
if [rij]=[] or not(type(rij,numeric)) then print(['i'=i,'j'=j], fij); rij:=undefined fi;
A[i+1,j+1,1]:=evalf(rhs(x)+rij*Sph*Cth);
A[i+1,j+1,2]:=evalf(rhs(y)+rij*Sph*Sth);
A[i+1,j+1,3]:=evalf(rhs(z)+rij*Cph);
od;od:
MESH(A);
end:

The procedure is not optimized, e.g.

- Cth, etc could be Vectors computed outside the loops

- Some small changes to use evalhf.

**###### EXAMPLES ######**

**f1 := x^2+3*y^2+4*z^2 - x*y - 2*y*z - 10:**

plots:-implicitplot3d(f1, x=-5..5, y=-5..5, z=-2..2);

**M:=Par3d(f1, x=0,y=0,z=0,5,40,40):**

PLOT3D(M);

**f2 := x^4+y^4+z^4-1:**

M:=Par3d(f2, x=0,y=0,z=0,5,40,40):

PLOT3D(M);

**M:=Par3d(f2, x=0,y=0,z=0, 5,40,40, 0..Pi, 0 .. Pi/3): #Plot half of the top only**

plots:-display(PLOT3D(M), scaling=constrained);

**M:=Par3d(f2, x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):**

N:=Par3d(f2+0.01, x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):

plots:-display(PLOT3D(M), color=red):

plots:-display(PLOT3D(N), color=green):

plots:-display(%,%%, orientation=[-40,65,10]);

**f3 := (x^2+y^2-1)^2+(z+sin(x*y+z))^4-120:**

plots:-implicitplot3d(f3, x=-4..4,y=-4..4,z=-5..5, numpoints=10000);

**Par3d(f3, x=0,y=0,z=0,5, 30,30):**

PLOT3D(%);

**Note.**

The procedure could be used to plot locally around a point (x0,y0,z0)

One may use the spherical coordinates (theta0,phi0) and then call the procedure taking theta0-a .. theta0+a, phi0-b, .. phi0+b for the trailing parameters

The spherical coordonates can be computed using:

ThetaPhi :=proc(x,y,z, X,Y,Z)
local r:=sqrt((X-x)^2+(Y-y)^2+(Z-z)^2);
['theta'=arctan(Y-y,X-x), 'phi'=arccos((Z-z)/r)]
end:

**ThetaPhi(10,20,30, 11,21,28);evalf(%);**