Help plotting equations in 3D

I am looking to plot the following in MATLAB to create a figure similar to the picture attached.

The equations are:

x=a*cos(phi);
y=a*sin(phi);
z=sqrt((c.^2./(r-(x.^2 + y.^2))).*(r-((1+epsilon).*(x.^2 + y.^2))));
where r=(a-A.*((y.^2.*(3.*x.^2-y.^2))./(x.^2+y.^2)).^B).^2;

And constants are:

a=1;
A=0.1339;
B=0.397;
c=1;
phi=0:0.001:2*pi;
epsilon=10.^-5;

NOTE: The original z equation is attached (I just solved for z and made a variable r to save space):

I have tried using plot3(X,Y,Z) using

[X,Y]=meshgrid(x,y);
z=peaks(6284);

because I found something similar in another tread.

The result is:

I appreciate any help you can offer!

[Star Strider]

The plot3 function requires three column vectors for its arguments. It doesn’t use meshgrid.

It will help if you post your code, or attach it if it is longer than about 20 lines. (Use the ‘paperclip’ icon and complete both steps.)

[Walter Roberson]

plot3() does accept meshgrid arguments.

Answer 1

[Mike Garrity]

Actually, I feel like we're missing a piece of the question here.

In your code, meshgrid and peaks(6284) are creating 2D arrays. That's what you want when you're drawing a parametric surface where each point is a function of 2 parameter values. But it doesn't look like you've got a parametric surface.

From the picture, it kind of looks like you're trying to draw a parametric curve where each point is a function of a single parameter value (perhaps phi?). In that case, you're going to be generating 1D arrays for your X, Y, and Z coordinates. As Star Strider and Walter said, once you have these 1D arrays, you're probably going to pass them to the plot3 function.

But that doesn't look like what you've got either. In your code, you're generating Z as a function of X and Y, rather than as a function of the parameter value phi. But that doesn't match your picture, because in your picture there are multiple Z values for the same X & Y pair.

On the other hand, your function actually looks like an implicit surface. This is where you're drawing a surface through all of the points where the left hand side of the equation are exactly equal to 1. But, if I code that up:

cla
n = 100;
[y,x,z] = ndgrid(linspace(-.5,.5,n), ...
                 linspace(-.5,.5,n), ...
                 linspace(-1,1,n));             
r = (a-A.*((y.^2.*(3.*x.^2-y.^2))./(x.^2+y.^2)).^B).^2;
f = ((x.^2+y.^2)./r).*(1+epsilon-z.^2/c.^2) + z.^2./c.^2;
isosurface(x,y,z,f,1)
view(3)
camlight
xlabel('X')
ylabel('Y')
zlabel('Z')

Then I get this:

which doesn't look at all like your picture.

Perhaps you could give us a bit more background?

BTW, you might find these examples from my blog helpful in figuring out how to draw these various types of curves and surfaces ( link1, link2, link3, link4). Perhaps one of those will give you an idea about what you're trying to draw.

[Joseph Nettesheim]

Dear Mike and others,

Thank you for taking the time to look at this. Mike, I am trying to create an implicit surface (I didn't know it was called that until your answer, thank you!). The links are very helpful in introducing me to isosurface.

I have attached the paper I am referencing for the equations. Equations 5 and 6 are of interest to me. I would eventually like to extract the data to make an STL file.

I tried to recreate Fig 2 using the code you wrote however then I get the error:

Error using matlab.graphics.primitive.Patch/set While setting the 'VertexNormals' property of Patch: This is not a valid PatchNormals value. Complex inputs are not supported

Error in isonormals (line 89) set(p, 'vertexnormals', n)

Error in isosurface (line 127) isonormals(x,y,z,data, p);

Error in BBC_iso (line 15) isosurface(x,y,z,f,1)

I have also attached the original code I wrote in trying to plot it with plot3, though I think I will take an approach similar to Mike's.

I will spend some time trying to figure this out with isosurface, and going through your blog examples to gain some more confidence in plotting implicit surfaces. Thanks again!

[Mike Garrity]

Thanks, the paper helps. You are definitely trying to draw an implicit surface here.

The error message is because the values of f aren't all real. I'm afraid that some places check for this in MATLAB Graphics and some don't. You could just add an "abs" or "real". That can be appropriate when the imaginary parts are very small, and have just resulted from roundoff errors. But I'm not sure that's the case here. I feel like this is telling us that we don't understand something about this equation.

[Joseph Nettesheim]

Hmm I am still having a tough time with this. Even using Gustavo Morales' Ezimplot3 function ( ezimplot3) I still cannot generate the intended geometries.

[Joseph Nettesheim]

I have sort of figured out a solution. It doesn't look as pretty as the figure from the paper, and for some reason at z=.2 (where z^2/c^2=1) the surface is defined for all x,y. Also the shape itself doesn't quite look right.

if true
    cla

    a=1;
    A=-3;
    B=1;
    c=.2;
    phi=0:0.01:2*pi;
    epsilon=10^-5;

    n = 100;

    x = linspace(-2*pi, 2*pi, n);
    y = linspace(-2*pi, 2*pi, n);
    z = linspace(-1, 1, n);

    [X, Y, Z] = meshgrid(x, y, z);

    term1num = (X.^2 + Y.^2);
    term1den = (a-A.*((Y.^2.*(3.*X.^2-Y.^2))./(X.^2+Y.^2).^3).^B).^2;

    term1 = term1num ./ term1den;
    term2 = (1+ epsilon -(Z.^2./c.^2));
    term3 = (Z.^2./c.^2);
    vals = (term1 .* term2) + term3;

    p = patch(isosurface(X, Y, Z, vals, 1));
    isosurface(X, Y, Z, vals,1);
    camlight
    view(3);

    xlabel('X');
    ylabel('Y');
    zlabel('Z');

end