Making a 3D plot and rotating it

Hi everyone, i have a code that plots 3D according to some parametric equations and rotates that 3D plot first around y-axis then around z-axis. But my problem is i get this wiggly 3D side surfaces instead of smooth one. I think it is because for loop plots 100 times. How can i get smooth side surfaces? Also is there any efficient method to write this code?

a = 0.0201; b = 0.0924;
alpha = 0.1063; beta = 0;

u = linspace(0,pi);
v = linspace(0,2*pi);
r = linspace(25,75);
[U V R] = meshgrid(u,v,r);

Z = -a*R.*cos(U).*cos(V);
Y = b*R.*cos(U).*sin(V);
X = sqrt((R.^2)-(Z.^2)-(Y.^2));

figure
hold on

for i = 1:100
    s = surf(X(:,:,i),Y(:,:,i),Z(:,:,i));
    direction = [0 1 0];
    rotate(s,direction,-rad2deg(alpha));
    direction = [0 0 1];
    rotate(s,direction,rad2deg(beta));
end

view(3)
xlabel('x');
ylabel('y');
zlabel('z');
grid on
shading interp
colorbar
hold off

[Sam Chak]

@Ahmet Turan, it appears that you have plotted 100 pieces circular disks, and then projected them in a 3D plot.

This probably explains non-smooth surface that you are seeing.

n = 100;
a = 0.0201; b = 0.0924; alpha = 0.1063; beta  = 0;
u = linspace( 0,   pi, n);
v = linspace( 0, 2*pi, n);
r = linspace(25,   75, n);
[U V R] = meshgrid(u, v, r);

Z = -a*R.*cos(U).*cos(V);
Y = b*R.*cos(U).*sin(V);
X = sqrt((R.^2) - (Z.^2) - (Y.^2));

figure
hold on
for i = 1:numel(u)
    s = surf(X(:,:,i), Y(:,:,i), Z(:,:,i));
    direction = [0 1 0];
    rotate(s, direction, -rad2deg(alpha));
    direction = [0 0 1];
    rotate(s, direction,  rad2deg(beta));
end
view(180, 0)
xlabel('x'); ylabel('y'); zlabel('z'); grid on, shading interp, colorbar, hold off

[Ahmet Turan]

@Sam Chak Yes, that is because i couldn't do it that way. How can i implement this in 3D?

[Sam Chak]

My geometry is relatively weak. 😅

Are you trying to generate a closed surface? Something like the Tower of Hanoi.

Answer 1

[Karim]

@Sam Chak as requested a seperate answer to show the "patch" approach

As mentioned in the comments, there remains a part to combine the two patches and close them up (i.e. we need to define a bit of logic to determine the connectivity).

Below i did it quickly by creating 2 patches. However, i needs a bit of further thinking about the logic to build the grid and connect everything properly. However i hope i was able to show the idea/concept.

alpha = 0.1063;
beta = 0;

a = 0.0201;
b = 0.0924;

numU = 30;
numV = 30;
numR = 15;

Rmin = 25;
Rmax = 75;

[GridShell,FacesShell] = GetPatchDataShell(numU,numV,numR,Rmin,Rmax,a,b);
[GridTop,FacesTop] = GetPatchDataTop(numU,numV,numR,Rmin,Rmax,a,b);

% rotate the grid
GridShell = (roty(-rad2deg(alpha)) * GridShell')';
GridShell = (rotz( rad2deg( beta)) * GridShell')';

GridTop = (roty(-rad2deg(alpha)) * GridTop')';
GridTop = (rotz( rad2deg( beta)) * GridTop')';

% plot the patch
figure
hold on
patch('Faces',FacesShell,'Vertices',GridShell,'FaceColor','r','FaceAlpha',0.25)
patch('Faces',FacesTop,'Vertices',GridTop,'FaceColor','r','FaceAlpha',0.25)
hold off
grid on
title('Combined plot')
view([-65 -10])

% seperate plot to only display the "top face"
figure
patch('Faces',FacesTop,'Vertices',GridTop,'FaceColor','r','FaceAlpha',0.25)
grid on
view([-55 10])
title('Top face')

% seperate plot to only display the "outer shell"
figure
patch('Faces',FacesShell,'Vertices',GridShell,'FaceColor','r','FaceAlpha',0.25)
grid on
view([-50 30])
title('Outer shell')

[Sam Chak]

+1 for you. 👍

[Ahmet Turan]

Thanks, i grasped the idea behind it!

Answer 2

[Karim]

To create a figure like the tower of hanoi, you can try something like the code below.

Have fun :)

edit: added the rotation

alpha = 0.1063;
beta = 0;

x0 = 0; % start at x = 0, y = 0, z = 0
y0 = 0;
z0 = 0;

R1 = 100; % starting ring radius
R2 = 20; % radius for the cross section

numRings = 5;

figure
hold on
for i = 1:numRings
    [X,Y,Z] = Create_3D_ring(x0,y0,z0,R1,R2);
    s = surf(X,Y,Z,'EdgeAlpha',0);

    rotate(s,[0 1 0],-rad2deg(alpha));
    rotate(s,[0 0 1], rad2deg(beta));
    
    % reduce radius of the ring
    R1 = R1-R2;
    % increase the z step
    z0 = z0 + 1.5*R2;
end
hold off
axis equal
grid on
view([25 20])

[Sam Chak]

+1 A Nice one! 👍

[Ahmet Turan]

It looks pretty cool!

But unfortunately i am not trying to make hanoi towers. I realized that i described my problem realy plain. I need to go into detail. So first of all i am trying to plot a 3d closed surface with these parametric equations:

Also R changes between [25e3 75e3]. For my case a = 0.0201 and b = 0.0924. When i implement these equations it needs to plot a shape like this:

Then finally i need to rotate this shape first around y axis ccw with angle of -alpha (alpha = 0.2241rad = 40.338 degree), then around z axis ccw with angle of beta (beta = 0) which can be implemented with these code:

direction = [0 1 0];
rotate(s,direction,-rad2deg(alpha))
direction = [0 0 1];
rotate(s,direction,rad2deg(beta))

Then after this rotation my shape needs to look like these images:

But as you guys told before instead of 3d plotting i got these image by placing elliptical disks together, this is why my shape is not smooth. So how can i plot this 3d closed surface?

[Karim]

one aproach is by using the "patch" command. For this you need to construct the set of vetices (i.e. coordinates for the grid points) and connectivity matrix (i.e. the 4 points that make up each of the little rectangles)

with these you can build the figure like the one you reffer to, i do not have the time to build the full figure, but here you can see the concept for the outer shell:

alpha = 0.1063;
beta = 0;

a = 0.0201;
b = 0.0924;

% number of segments over different dirrections
% this determines the number of rectangles
numU = 1; % 1 in the u parameter to construct the outer shell
numV = 20;
numR = 10;

% parameter values to evaluate the grid
U = linspace(0,pi,numU);
V = linspace(0,2*pi,numV+1); V(end) = [];
R = linspace(25,75,numR);

% allocate the grid matrix
MyGrid = zeros(numU*numV*numR,3);

% evaluate the grid at the different points
% this can be vectorized, i made a loop to not over complicate things
i = 0;
for vi = 1:numV
    for ui = 1:numU
        for ri = 1:numR
            
            Z = -a*R(ri).*cos(U(ui)).*cos(V(vi));
            Y =  b*R(ri).*cos(U(ui)).*sin(V(vi));
            X = sqrt((R(ri).^2)-(Z.^2)-(Y.^2));
            
            i = i+1;
            MyGrid(i,:) = [X Y Z];
            
        end
    end
end

% use logic to build the conectivity matrix
Faces = [1:(numR-1);2:numR;numR+(2:numR);numR+(1:(numR-1))]';
Faces = repmat(Faces,numV,1) + repmat(repelem(numR:numR:(numR*numV),numR-1)',1,4);
OffsetIdx = Faces > (numU*numV*numR); % quick work arround to close the figure
Faces(OffsetIdx) = Faces(OffsetIdx) - (numU*numV*numR);

% rotate the grid
MyGrid = (roty(-rad2deg(alpha)) * MyGrid')';
MyGrid = (rotz( rad2deg( beta)) * MyGrid')';

% plot the patch
figure
patch('Faces',Faces,'Vertices',MyGrid,'FaceColor','r','FaceAlpha',0.25)
grid on
axis equal
view([-35 35])

% EDIT: extra figure to show the result witouth "axis equal", since this
% seems to be he case on your other figures
figure
patch('Faces',Faces,'Vertices',MyGrid,'FaceColor','r','FaceAlpha',0.25)
grid on
view([-50 30])

[Sam Chak]

Hi @Karim. Thanks for showing how to use the patch() function.

Unsure if the math describes the wrapped part in the circled area.

By the way, I have found one MVP answered an almost similar question here:

https://www.mathworks.com/matlabcentral/answers/1744660-3d-plotting-of-given-parametric-equations

[Karim]

Hi @Sam Chak, yes you can obtain the "top" face in a similair manner. I will try to add the code later when i have a bit more time. To obtain the "top" face, the trick will be to loop again over the equations to evaluate the grid points on that "face" and determine the logic behind the connectivity. They can then be combined to form a single "patch".

The reason why i did it using patch (and the loop which could be vectorized) is that meshgrid also returns the inner faces which i guessed the OP didn't want. I assumed he only wanted the outer shell.

[Karim]

@Sam Chak below you can find the idea where i also added the top face...

However there remains a part to combine the two patches and close them up (i.e. we need to define a bit of logic to determine the connectivity).

Below i did it quick by creating 2 patches. However, i needs a bit of further thinking about the logic to build the grid and connect everything properly. However i hope i was able to show the idea/concept.

alpha = 0.1063;
beta = 0;

a = 0.0201;
b = 0.0924;

numU = 50;
numV = 50;
numR = 20;

Rmin = 25;
Rmax = 75;

[GridShell,FacesShell] = GetPatchDataShell(numU,numV,numR,Rmin,Rmax,a,b);
[GridTop,FacesTop] = GetPatchDataTop(numU,numV,numR,Rmin,Rmax,a,b);

% rotate the grid
GridShell = (roty(-rad2deg(alpha)) * GridShell')';
GridShell = (rotz( rad2deg( beta)) * GridShell')';

GridTop = (roty(-rad2deg(alpha)) * GridTop')';
GridTop = (rotz( rad2deg( beta)) * GridTop')';

% plot the patch
figure
hold on
patch('Faces',FacesShell,'Vertices',GridShell,'FaceColor','r','FaceAlpha',0.25)
patch('Faces',FacesTop,'Vertices',GridTop,'FaceColor','r','FaceAlpha',0.25)
hold off
grid on
view([-55 10])

% seperate plot to only display the "top face"
figure
patch('Faces',FacesTop,'Vertices',GridTop,'FaceColor','r','FaceAlpha',0.25)
grid on
view([-55 10])

[Sam Chak]

Hi @Karim, really appreciate your contributions on this topic.

You should post in New Answer so that I (and others) can vote it 👍👍👍.