how to create a volume from the revolution of a variable area trapezoid

I'm looking for the work flow of how to create a volume by revolving a variable area trapezoid around the z axis. I have a formula for the area of the trapezoid and it is a function of the radial distance from the origin. Both the height and this radial distance vary with the angle of revolution.

The figure immediately below shows the elliptical path that the rotation takes. Because the trapezoid area varies with the angle of revolution, the path traced in the x-y plane is an ellipse and not a circle. The figure below shows the path in the x-y plane:

This next figure shows the x-z plane which shows the variable-area trapezoid at phi = 0 degrees and phi = 180 degrees. The rotation occurs around the vertical axis that goes through the point "S". Notice the large reduction in area of the trapezoid.

The radial line length, p, to the ellipse edge is given by,

where θ is a constant value as are w and R. At θ = 45 degrees and R = 0.0006, w = 0.0007908.

NOTE: Figure 2 is not to scale with Fig. 1.

What I am trying to do is to verify an equation that I have developed for calculating this volume. I want to make sure that my equation, in fact, produces the volume of the revolved variable-area trapezoid. The area of the trapezoid is given by

and h is given by:

where B is a constant and is a constant where p is evaluated at ϕ = 90 degrees. For the values of θ, w and R given above, B = 0.867.

The volume then is given by,

The next figure shows an incremental volume diagram which is then integrated as in the equation immediately above,

And so the final result should look something like the blue portion of each drawing below. How would I go about implementing this in Matlab?

[John D'Errico]

Sorry. I was going to write some code for surfaces of revolution a few years ago. Did not finish it, probably because I was not sure what I wanted the result to be, and what degree of flexibility it needed. Sadly, there are a lot of tools I've partly built over the years, so I have a moderately large bone yard of code, round tuits. :(

You can get what you want, (I think, if Im guessing right about what you want) by starting with a cylinder, then a simple linear transformation of it would give you the blue thingy.

Anyway, what result would you want to see? I.e., a triangulated surface? Just a pretty picture? Both? Something else?

[Robert Demyanovich]

I have updated the post which should answer a number of questions. However, I don't know how to use latex in these posts. Perhaps you could provide me with some help on that. i'll hold off providing an equation for p because it is quite complicated. I guess I could insert images if latex doesn't work in this forum.

[Walter Roberson]

It looks like @Matt J must have edited the latex into existence.

To insert latex, start in text mode anywhere on a line, and either at the beginning of a line or after a space, type a $ and then immediately the latex code, and finish by typing a $ again. The first thing after the initial $ must not be a space, and the last thing before the closing $ must not be a space. Upon typing the closing $ the latex will immediately be rendered.

If there is some problem with the latex entry or if you want to edit it, position your cursor immediately after the rendered latex, and press your backspace key . That will change the rendered version to plain text version again. To re-render, position immediately after the closing $ and press space.

Latex mode does not work in code blocks, or in "code example" blocks.

The supported latex in latex mode here exceeds that which is supported by MATLAB 'latex' interpreter, or which is supported by "insert equation" in Live Script. Unfortunately there does not appear to be documentation about what exactly is supported.

[Robert Demyanovich]

Thank you for that brief tutorial on latex.

[Robert Demyanovich]

I have changed the first figure to more clearly show the radial line, p, which starts at the point "S".

[Matt J]

But is not a function of ϕ? It is constant?

[Robert Demyanovich]

Yes. That is correct. is not a function of ϕ. It is constant and is evaluated at a value of ϕ equal to 90 degrees.

[Robert Demyanovich]

Values of constants are now included in the original post.

[Mathieu NOE]

maybe I'm reading the post a bit too fast , but if the goal is just to generate a plot, the blue section is basically just and elliptic tube that intersect two planes ... do we really need all this complicated maths ?

[William Rose]

@Robert Demyanovich,

[Edit: fix typos in the text. Code and formulas are not changed.]

[Edit: add equation for, and calculation of, hs.]

Thank you for posting an interesting problem.

You say in your original post that you want to "create a volume by revolving a variable area trapezoid around the y axis". However, you say that Figure 1 (figure with an ellipse) in the original post shows the x-y plane. In that case, the +z axis points upward out of that plane of Figure 1, and goes through point S. Is that correct? (It would help if each 2D figure included labels for axes.) You say Fig. 2 in the original post shows the x-z plane, and you say that the axis of rotation passes vertically through point S in Figure 2. Therefore it seems the trapezoid revolves about the z axis. But that contradicts what you said earlier. Please clarify.

You write in the original post

It seems to me, based on your other explanations and comments, that the area given by the equaiton above is the area when , because it is only when that . It seems to me that a more general formula for the area is

Am I understanding the problem correctly?

Your equation for p/R appears on my screen as a bunch of text that looks like it should be a Latex equation. So I will try to make it appear as a formatted equation:

And w, R, and θ are simply constants:, R=0.0006, w=0.0007908. Then and the equation above can be written

(Check my math there: did I do that right?) Then we can plot of versus ϕ, and we can plot in the x,y plane:

phi=0:pi/20:2*pi; w=7.908e-4;  R=6e-4;
% Compute p(phi):
p=((w/2)*cos(phi)+(R^2*(sin(phi).^2+0.5*cos(phi).^2)-(w^2/2)*sin(phi).^2).^0.5)./(sin(phi).^2+0.5*cos(phi).^2);
% Plot results
figure; 
subplot(211), plot(phi,p,'-r.') 
grid on; xlabel('\phi'); ylabel('p'); title('p(\phi)');
subplot(212), plot(p.*cos(phi),p.*sin(phi),'-r.');
grid on; xlabel('x'); ylabel('y'); title('p(x,y)')

The lower plot above looks consistent with the plot of p(x,y) in Fig.1 of your original post.

Your Figure 3 from the original post is below:

I am not sure why use you use subscripts h_i and h_k, above. I think hi=hk when the angle is small, so it would make more sense to label both h's in the figure as . You provide the following equation for the total volume:

We can infer from the equation above that you believe that the volume of the 3-dimensional wedge in the figure above is

I do not think that is correct. IF the thick end of the wedge were perpendicular to the wedge long axis (and it's not perpendicular), then the volume would be

where p is a function of ϕ, and therefore dV is a function of ϕ. That's what I think, but you should check my work.

The thick end of the wedge is not perpendicular to the long axis of the wedge, so the preceding formula needs a correction factor which equals unity when the end is perpendicular, and the correction factor is less than 1, but more than 0, otherwise. The end of the wedge is perpendicular to the long axis only when and when .

You say and . Then . It follows from the equation I provided above for that . Therefore .

[Matt J]

To insert latex, start in text mode anywhere on a line, and either at the beginning of a line or after a space, type a $ and then immediately the latex code, and finish by typing a $ again.

Or enter latex in the formula editor,

This is good when you have multi-line formulas.

[Robert Demyanovich]

@William Rose: Thank you for taking the time to look at this problem. You are quite right that the vertical axis around which the variable area trapezoid rotates is the "z" axis and not the "y" axis (fixed in the original post). I have also fixed the general equation for the trapezoid area and the equation for as recommended. Thank you for pointing out these inconsistencies.

Your eqn. 4 where is set equal to 0.5 is correct. Your second plot is essentially the same as the first plot in my post. So far, so good.

With regards to my 3rd plot, is the height of the trapezoid at ϕ = 0 degrees. The notation is meant to show a height for a trapezoid that is a very short distance away. The area of trapezoid k is slightly different. It is less than the area of trapezoid i, which has the largest area of any trapezoid slice. The figure was meant to show that the differential thickness by which the trapezoid area is multiplied is .

The equation that you show for dV I believe is correct and as you state is inferred from the integral equation in the original post. The thinking is that the volume is simply the integral of the area of the trapezoid (which is a function of ϕ) multiplied by the differential thickness, which I believe is .

I do not understand what you mean by "IF the thick end of the wedge were perpendicular to the wedge long axis (and it's not perpendicular), then the volume would be". I don't know what this means. My understanding is that the volume is the entire rotation of the variable area trapezoid over 2pi. Then you also wrote "where p is a function of ϕ, and therefore dV is a function of ϕ". But the fact that dV is a function of ϕ is already indicated in your Eq. 6, so I guess I am missing the point.

By thick end of the wedge, do you mean the height? At a specific value of ϕ the height is given by:

so the height of the wedge varies as the trapezoid is rotated around the z axis. I'm not sure why a correction factor would be needed. I'm guessing I'm not understanding what you are saying with regards to the thick end of the wedge. My thinking is that for all trapezoid slices where ϕ is less than 90 degrees or greater than 270 degrees, . When ϕ is equal to 90 or 270 degrees, then . For , (see the second figure in my original post which shows the trapezoids at ϕ equal to 0 degrees (towards the right of point "S") and ϕ equal to 180 degrees (towards the left of point "S").

[Robert Demyanovich]

@William Rose: Your derivation for is verified. For the given values in my original post, .

[William Rose]

@Robert Demyanovich, thank you for checking my work.

[Edit: run the code so that it generates a figure. Change my conclusion: Robert D is right, no adjustment is needed.]

I now recognize that you are right and I was thinking incorrectly. As you said earlier, there is no need for an adjustment factor for the differential volume, dV.

The figure below shows a view of the object from above, i.e. a slice through the object in the x-y plane, through z=0. Three triangles are highlighted to show that the edge of the triangle that is opposite point S is generally not perpendicular to the long axis of the triangle. But no adjustment to the volume formula is needed, because the width of the triangle, perpendicular to the long axis, is p(t)*dt.

N=40;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]'-dt/2;    % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
px=p.*cos(t); py=p.*sin(t);
figure
plot(px,py,'--k'); hold on
plot([0;px(1:2);0],[0;py(1:2);0],'-r')
plot([0;px(6:7);0],[0;py(6:7);0],'-g')
plot([0;px(11:12);0],[0;py(11:12);0],'-b')
xlabel('x'); ylabel('y')
text(0,-4e-5,'S')
grid on; axis equal

[Robert Demyanovich]

I could not find an attached figure.

[William Rose]

@Robert Demyanovich, there was no figure because I forgot to run my code before posting the comment. I have fixed that now, and I have adjusted my conclusion - see edited comment above.

[Paul]

I'm confused about the volume integral that's shown above, which is

(unless that's changed in the comments below). I'm trying to relate that back to a simpler 2D case to express the area enclosed by p(phi), which I guess would be the area of solid in the z=0 plane. It seems like the area integral analagous to that volume integral would be

but that wouldn't be correct, would it? It's missing a factor of two, i.e., the correct area integral in polar coordinates is

Similarly, would the volume integral have to account for the change in height going from h(phi) to h(phi + dphi) (I'm not sure about this).

[Robert Demyanovich]

See the third figure. The volume is being estimated based on a variable area trapezoid multiplied by the differential thickness, which I believe is . Also, note that the revolution is around a vertical axis through the point S.

[Paul]

Does the proposed volume integral work for simple case?

For example, instead of variable area trapezoid, let's use a constant area trapezoid. And let's use a simple trapezoid, e.g., a rectangle.

Let hs = 1, h(phi) = 1, and p(phi) =2, so that our trapezoid is really a rectangle with area h*phi = 2. The solid of revolution is a cylinder, and the volume of that cylinder is pi*r^2*h = 4*pi.

If we use the formula for the volume of a solid of revolution(link), we'd have

syms f(x) g(x)
f(x) = 0.5; g(x) = -0.5;
V = 2*sym(pi)*int(x*abs(f(x)-g(x)),x,0,2)

V = 

Using the formula above, we have

syms A(phi) p(phi)
A(phi) = 2; p(phi) = 2;
V = int(A(phi)*phi,phi,0,2*sym(pi))

V = 

[Robert Demyanovich]

@Matt J, William Rose, Paul:

I'm beginning to think that my equation for the differential volume is incorrect and it is more complex than originally thought. What I should be doing is looking at the sector (wedge) as the rotation occurs and multiplying it by the height (like what would be done for a cylinder). It will now take a double integral because I need to integrate over the length "r" from the point S to the edge of the ellipse where r = p. So, maybe something like this:

and the integral then is,

where r at the point S is equal to 0 and at the ellipse edge is equal to p. If I assume that the change in height from r = 0 to r = p is linear along any radial line, then

I think this will result in something similar to what William Rose initially proposed for dV.

Thoughts?

[Matt J]

I agree that,

would be the correct itegral formula.

[Robert Demyanovich]

So, when I work out the double integral, I get

substituting in the following expession for

yields,

which becomes

Evaluating the inner integral yields,

or

or more simply,

This expression looks like what William Rose has been proposing, whereas, I originally and incorrectly proposed:

[William Rose]

@Robert Demyanovich,

What you are proposing is what we already did. That is the integral I solved to get the formula which I posted many comments ago:

The next step is to solve the integral

where

I don't know how to solve that integral, and neither does Wolfram Alpha. I tried.

Therefore you solve it numerically. My code already does this, by two methods, which give the same answers, as shown in one of my other comments. Method 1 (called method 2 in my other comment) is to compute each dV using the formula for dV, above. Method 2 (called method 3 in my other comment) is to compute dV as the volume of the convex hull around the six vertices of the wedge dV. (Method 1, in the other comment, is the tetrahedralization approach. @Matt J says it also gives the same answer, when an earlier typo is fixed.)

[Robert Demyanovich]

@William Rose: Yes, you were correct and I was wrong in my original post. In the post from today I did acknowledge that you originally proposed the equation. I needed to come around to the correctness of that equation by deriving it myself.

I have accepted your answer. Thanks so much for your help.

[William Rose]

@Robert Demyanovich, you're welcome.

[Robert Demyanovich]

Hopefully, one last question. If the height used to calculate the differential volume/volume is linear, then why is the height concave (in comparison to the red trapezoid) as shown in the following image generated for the x-z plane from the MatLab code. The trapezoid shown is for an azimuthal angle equal to 0 degrees.

[Matt J]

The blue area doesn't show us the height of any single one of the wedges. It's the height of the projection of the entire volume on the x-z plane.

The concavity of an object is not preserved under perspective projection. Think of a bowl-shaped solid. In 3D, it is concave, but its projection is convex.

[Robert Demyanovich]

First off, I made a mistake, the azimuthal angle is 0 degrees not 180 degrees. I guess I thought that if I look at the shape directly from the side (i.e. the x-z plane), then what I would see is the largest trapezoid which is at 0 degrees. I think if I'm looking from the side at the x-z plane, that the trapezoid at 0 degrees would have the highest height as viewed by moving my eyes along the x-axis. At other azimuthal angles, the height is decreasing.

[Matt J]

Depends how far away the camera is. If you stood at the foot of the Seattle Space Needle, you wouldn't be able to see the tip of the needle, even though it is the highest point. The observation deck would be in the way.

[Robert Demyanovich]

Well, in the image posted, what is it that is above the red trapezoid line at the top? The image is exactly in the x-z plane, there is no curvature due to the y-plane. So, which heights at other azimuthal angles could be larger than the height at d$\phi$ = 0, since by equation, those height are always less. That radial line is kind of like a backbone and everything drips off the upper red line, for example, in the x-y plane. I guess I just don't see what you are saying regarding the image that I posted (not the Seattle space needle example, which I don't think applies, because in the the image I posted, I can see all of it rather easily. That wouldn't be the case if I was at the foot of the space needle trying to get "perspective" on the entire space needle). It doesn't matter where I am located on the z plane between +7E-04 and -7E-04, I can see all of the image and the greenish blue is pretty much always above the upper red trapezoid line or below the lower red trapezoid line. Of course, if I have my face pressed up on the monitor, then maybe I can't see the entire image. But backing away even a little bit allows me to view the entire image with greenish blue above and below the red lines.

[Matt J]

So, which heights at other azimuthal angles could be larger than the height at d$\phi$ = 0,

All of them, depending again on where the camera is positioned. Even though those parts of the object have a lower height in 3D, they become bigger from the camera's perspective if they are closer to the camera, while the parts farther away get smaller. The bottom line is, we need to know where you put the camera.

doesn't matter where I am located on the z plane

Not sure what the "z-plane" is. Do you mean you've translated the camera up and down in z, but kept it at the same (x,y) coordinates? This might be a good time to mention that we haven't seen the code that generated these perspectives. It limits our ability to assess the rendering precisely.

[Robert Demyanovich]

Okay.

The image that I posted was generated using the view command in the following line,

axis equal; view(0,0); hold on; camlight

Originally, it was view (-45,25) and I thought that changing it to view(0,0) would provide a camera shot producing an exact cross section. Is that not the case?

[Robert Demyanovich]

It looks like William Rose commented on this issue several days ago. What he says make sense.

[Matt J]

It looks like William Rose commented on this issue several days ago

You haven't provided a link to his comment, but you might be refering to this comment, describing that the top surface of your volume is not planar. An implication of this is that the central spine of your volume (at ) is in fact not the highest. For every fixed x, the points at the rim of your volume are higher than at the center.

That does indeed explain what you are seeing, but I still wanted to illustrate how, even when the top surface is planar, you can still see the concave envelope, depending on your display settings. In this first figure below, I plot the cutoff cylinder that I showed you how to generate using AxelRot. The top and bottom surfaces are planes perpendicular to the x-z plane, and therefore the central spine of the object really is maximally high in this case for any fixed x.

 figure;
 plotVolume();
 view(-30,20)

Accordingly at view(0,0), the central trapezoid is snug against the blue envelope of the object, as long as the axes' Projection property is orthographic (the default).

 figure;
 plotVolume();
 view(0,0)

However, if we switch to perspective projection, which is how peoples' eyes more naturally see the world, the central trapezoid will never be flush with the upper an lower envelopes of view(0,0), and you get the concave envelope. This is because the edge closer to the observer will always appear higher than the central spine.

 figure;
 plotVolume();
 view(0,0)
 set(gca,'Projection','perspective');

function plotVolume()


    %Volume parameters (independent)
    N=100; %discretize the perimeter into N points
    diam=6; %tube diameter
    ho=1;
    hi=4;
    xS=-2; %x-coordinate of S
    
    %Volume parameters (derived)
    dh=(hi-ho)/2;
    angle=atand(dh/diam);
    axMajor=hypot(diam,dh);
    hS=ho+2*dh/diam*(xS+diam/2);
    
    %Compute boundary vertices
    t=linspace(0,360,N+1)'; t(end)=[];
    V=0.5*[axMajor*cosd(t),diam*sind(t),0*t];
    
    XYZ1=AxelRot( (V+[0,0,+ho/2])' ,+angle,[0,-1,0], [-diam/2,0,+ho/2])';
    XYZ2=AxelRot( (V+[0,0,-ho/2])' ,-angle,[0,-1,0], [-diam/2,0,-ho/2])';
    XYZ=[XYZ1;XYZ2];

    %Display
    trisurf(convhulln(XYZ), XYZ(:,1), XYZ(:,2), XYZ(:,3), ...
        'FaceColor', 'cyan', 'EdgeColor', 'none','FaceAlpha',0.3);
    axis equal; xlabel x; ylabel y; zlabel z
    camlight

    
    %Movie
    
    i=1;
    A=XYZ1(i,:);
    B=XYZ2(i,:);
    C=[xS,0,-hS/2];
    D=[xS,0,+hS/2];
    
    args=num2cell([A;B;C;D;A],1);
    [x,y,z]=deal(args{:});
    
    line(x,y,z,Color='r',LineWidth=2); 

    view(0,0)
    camdolly(0,+30,0,'fixtarget','data');

end

[Robert Demyanovich]

Thanks for the explanation.

Answer 1

[William Rose]

@Robert Demyanovich

{Edit 10/20/2025: Fix error in code below by changing a 1 to a 3. It was

[hst(1),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

and it should be, and now is

[hst(3),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

}

I like the elegant solution of @Matt J.

You indicate in your comments that you are also interested in solving the problem by explicit rotation of the trapezoid.

Here is code that that generates a wire frame view of the shape, by rotating the trapezoid.

N=40;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]';         % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
h=B*p;  % [length]
hs=B*sqrt(R^2-w^2/2);
% Compute arrays of 3D points
pt=[p.*cos(t),p.*sin(t),h/2];      % x,y,z coords of points along top of shape
pb=[p.*cos(t),p.*sin(t),-h/2];     % x,y,z coords of points along bottom of shape
hst=[0,0,hs/2]; hsb=[0,0,-hs/2];    % x,y,z coords of top, bottom central points
% Plot the trapezoid in 3D at various rotation angles 
figure; hold on
for i=1:N
    plot3([hst(1),pt(i,1),pb(i,1),hsb(1),hst(1)],...
          [hst(2),pt(i,2),pb(i,2),hsb(2),hst(2)],...
          [hst(3),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')
end
plot3(pt(:,1),pt(:,2),pt(:,3),'-g')  % curve connecting top points
plot3(pb(:,1),pb(:,2),pb(:,3),'-b')  % curve connecting bottom points
xlabel('x'); ylabel('y'); zlabel('z'); 
hold off; axis equal; grid on; view(-45,25)

If you run the code locally, you will be able to rotate the 3D plot by clicking on the 3D rotate tool at top of the plot, then click and drag inside the plot.

Wrap a surface around the points, compute the volume, plot the surface:

pAll=[pt;pb;hst;hsb];   % all the 3D points
[k,v]=convhull(pAll);   % wrap the points and compute volume
fprintf('Volume=%.3e\n',v)

Volume=1.323e-09

figure
trisurf(k,pAll(:,1),pAll(:,2),pAll(:,3),'FaceColor','c','EdgeColor','none')
axis equal; camlight

I am attaching a script which has the code above, with more comments and an additional plot.

[Robert Demyanovich]

@William Rose: First: Do you now agree that my equation for the volume is correct or do you think that correction factors are still needed?

Second: Thanks so much for providing code on how to accomplish what I am looking for. On the surface this looks good, but I need some time to check it out in depth. Also, I am going to modify the code a bit to allow me to use different values of θ, which will change the values of w and B, but not R. Then I can check the calculated volumes for different scenarios.

Any chance you could help me create a video showing a trapezoid rotating around the vertical line through "S" or point me in the right direction on how to do this?

[William Rose]

@Robert Demyanovich,

1. You are right that no adjustment to the volume formula is needed.
2. You're welcome.
3. Yes. Later.

[William Rose]

@Robert Demyanovich,

The attached scripts shows the shape growing frame by frame. The scripts do not yet make a video, but I'll add that later.

Script 1 shows the plot of the convex hull surface, as the trapezoid rotates and the shape grows, The partially completed object is non-convex, when the trapezoid is more than halfway around, until it is complete. This causes an undesirable result, because convhull() drapes a convex hull around the points. The resulting surface hides the wedge that has yet to be completed in this squished round of cheese. A complicated solution would be to define and plot the all the polygons on the surface at each iteration, using trapezoidal polygons around the sides and triangles on the top and bottom. That would be doable, but a lot of tedious work.

Script 2 shows the wire frames as the trapezoid rotates.

[William Rose]

@Robert Demyanovich,

This script saves the successive plots of the wire frames as frames in a video file. The video file is also attached.

Sometimes when I make videos like this, I save 1 or 2 seconds worth of frame 1 at the start, and 1 or 2 seconds worth of the last frame, to provide bookends on the video. You could easily adjust the script to do that if you wish.

[Robert Demyanovich]

@William Rose: I am calculating a different volume than what your Matlab code yields. I'm getting 2.087E-09 whereas your code gets 1.347E-09 after increasing N to 100. Increasing N only increases to around 1.352E-09. I used online Wolfram-Alpha to calculate the integral but because of length limits on W-A I had to integrate the terms separately and then add them together. If you have the time, I'm curious what value of the volume you calculate using the equation for the volume. Thanks.

[William Rose]

@Robert Demyanovich,

[edit: correct typo]

I estimate volume two ways: with the convhull command

[k,v]=convhull(pAll);

and by the formula for dV which I suggested.

% Estimate total volume using formula for dV
dV=(B*p(1:N).^3/3+hs*p(1:N).^2/6)*dt;
Vtot=sum(dV);

Results, when I change N=number of steps in one full rotation, then re-run the script:

fprintf('convhull() volume=%.3e, sum(dV) volume=%.3e.  N=%d.\n',v,Vtot,N)

>> rotateTrapezoid

convhull() volume=1.323e-09, sum(dV) volume=1.291e-09. N=40.

>> rotateTrapezoid

convhull() volume=1.345e-09, sum(dV) volume=1.291e-09. N=80.

>> rotateTrapezoid

convhull() volume=1.351e-09, sum(dV) volume=1.291e-09. N=180.

>> rotateTrapezoid

convhull() volume=1.352e-09, sum(dV) volume=1.291e-09. N=360.

>> rotateTrapezoid

convhull() volume=1.352e-09, sum(dV) volume=1.291e-09. N=720.

Discussion

• My sum(dV) volume is a lot less than what you report, and pretty close to the convex hull volume.
• volume=sum(dV) is constant (1.291e-9) as the number of slices (N) grows and each slice gets smaller.
• The volume of the convex hull increases slightly from N=40 to N=180 slices, and after that it is stable to four significant figures.
• There is a persistent difference, even at the highest N values, where the convex hull volume is ~5% more than sum(dV) volume. I don't know why.

[Robert Demyanovich]

Thank for your comment. Are you now of the mind that your original equation for dV is correct? Or are you just using it as a point of reference.

[Matt J]

Now that we know the solid is not convex and convhull is not applicable, you will have to use a tetrahedralization of the solid to compute the volume. Here is a modification of @William Rose's original loop for doing so,

N=360;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]';         % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
h=B*p;  % [length]
hs=B*sqrt(R^2-w^2/2);
% Compute arrays of 3D points
pt=[p.*cos(t),p.*sin(t),h/2];      % x,y,z coords of points along top of shape
pb=[p.*cos(t),p.*sin(t),-h/2];     % x,y,z coords of points along bottom of shape
hst=[0,0,hs/2]; hsb=[0,0,-hs/2];    % x,y,z coords of top, bottom central points
V=0;
wrappedList=[1:N,1];
for k=1:N
    
    i=wrappedList(k);  
    j=wrappedList(k+1);  
    
    x1=[hst(1),pt(i,1),pb(i,1),hsb(1)];
    y1=[hst(2),pt(i,2),pb(i,2),hsb(2)];
    z1=[hst(3),pt(i,3),pb(i,3),hsb(3)];
    x2=[pt(j,1),pb(j,1)];
    y2=[pt(j,2),pb(j,2)];
    z2=[pt(j,3),pb(j,3)];
    
    TR=delaunayTriangulation([x1,x2]',[y1,y2]',[z1,z2]');
    V=V+tetmeshVolume(TR);

end

volume=V

volume = 1.2903e-09

function V_total = tetmeshVolume(TR)
    % TR can be a delaunayTriangulation or struct with .ConnectivityList and .Points

    T = TR.ConnectivityList;
    X = TR.Points;

    % Extract vertex coordinates for all tetrahedra
    v1 = X(T(:,1), :);
    v2 = X(T(:,2), :);
    v3 = X(T(:,3), :);
    v4 = X(T(:,4), :);

    % Compute per-tetra volumes using vectorized cross and dot
    V = abs(dot(v1 - v4, cross(v2 - v4, v3 - v4, 2), 2)) / 6;

    % Total volume
    V_total = sum(V);
end

[William Rose]

@Matt J, @Robert Demyanovich,

Thank you @Matt J for showing how to tetrahedralize a dV, and find its volume as the sum of the volumes of the tetrahedra that make up each dV. The total volume you get by this method is quite a bit less than the volume I got by using a formula for each dV: 1.2144e-09 versus 1.291e-09. So I wonder if I made a mistake when I derived my formula for dV by integrating.

To get more insight, I compute the dVs and total volume three ways:

1. @Matt J's tetrahedralization way
2. Using the formula I derived by integrating, with the average value of p for each wedge.
3. Using convhull() to estimate the volume of each dV

N=360;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]';         % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
h=B*p;  % [length]
hs=B*sqrt(R^2-w^2/2);
% Compute arrays of 3D points
pt=[p.*cos(t),p.*sin(t),h/2];      % x,y,z coords of points along top of shape
pb=[p.*cos(t),p.*sin(t),-h/2];     % x,y,z coords of points along bottom of shape
hst=[0,0,hs/2]; hsb=[0,0,-hs/2];    % x,y,z coords of top, bottom central points

% 1. Volume by Matt J's tetrahedralization method
dV1=zeros(N,1);
wrappedList=[1:N,1];
for k=1:N
    
    i=wrappedList(k);  
    j=wrappedList(k+1);  
    
    x1=[hst(1),pt(i,1),pb(i,1),hsb(1)];
    y1=[hst(2),pt(i,2),pb(i,2),hsb(2)];
    z1=[hst(1),pt(i,3),pb(i,3),hsb(3)];
    x2=[pt(j,1),pb(j,1)];
    y2=[pt(j,2),pb(j,2)];
    z2=[pt(j,3),pb(j,3)];
    TR=delaunayTriangulation([x1,x2]',[y1,y2]',[z1,z2]');
    dV1(i)=tetmeshVolume(TR);
end
volume1=sum(dV1);

% 2. Volume using a formula for dV
pbar=(p(1:N)+p(2:end))/2;           % mean p for each wedge
dV2=(B*pbar.^3/3+hs*pbar.^2/6)*dt;
volume2=sum(dV2);

% 3. Volume by convhull() of each dV
dV3=zeros(N,1);
for i=1:N
    wedgeVerts=[hst;hsb;pt([i,i+1],:);pb([i,i+1],:)];
    [~,dV3(i)]=convhull(wedgeVerts);
end
volume3=sum(dV3);
fprintf('N=%d: volume 1=%.4e, volume 2=%.4e volume 3=%.4e\n',N,volume1,volume2,volume3)

N=360: volume 1=1.2401e-09, volume 2=1.2905e-09 volume 3=1.2903e-09

function V_total = tetmeshVolume(TR)
% TR can be a delaunayTriangulation or struct with .ConnectivityList and .Points
T = TR.ConnectivityList;
X = TR.Points;
% Extract vertex coordinates for all tetrahedra
v1 = X(T(:,1), :);
v2 = X(T(:,2), :);
v3 = X(T(:,3), :);
v4 = X(T(:,4), :);
% Compute per-tetra volumes using vectorized cross and dot
V = abs(dot(v1 - v4, cross(v2 - v4, v3 - v4, 2), 2)) / 6;
% Total volume
V_total = sum(V);
end

Here are the total volume results, for different values of N (multiply each volume by 1e-9):

N volume 1 volume 2 volume 3

40 1.2144 1.2747 1.2638

80 1.2338 1.2866 1.2839

360 1.2401 1.2905 1.2903

720 1.2403 1.2906 1.2906

Remember that vol.1=Matt J's tetra method; vol.2=formula; vol.3=convhull of each dV.

The tetrahedron method gives a smaller volume by about 4%. The other two methods converge to one another as N gets large. I don;t know why.

Does the difference between the methods vary as a function of rotation angle? To find out, plot dV2/dV1 and dV3/dV1, versus wedge number, when N=360. Therefore wedge number=degrees.

plot(1:N,dV2./dV1,'rx',1:N,dV3./dV1,'b+');
xlabel('Degrees'); ylabel('Ratio'); legend('dV2/dV1','dV3/dV1')
grid on; xlim([0 360]); xticks([0 90 180 270 360])

The plot shows that methods 2 and 3 give a 50% larger value than the tetrahedron method, at 180 degrees. At this angle, the wedges are very small. I'm not sure which method is right, and why the other(s) are wrong.

[Matt J]

I had a typo in my code for z1. Should have been,

 z1=[hst(3),pt(i,3),pb(i,3),hsb(3)];

at which point all 3 methods agree. And, convhull does seem to make more sense.

[Sam Chak]

Hi @William Rose

+1 👍 to you and @Matt J. Loved all three solutions as they are practical and well-explained. Initially, I was confused because I thought the OP wanted to estimate the air volume of a section of the circular ventilation duct with a 45° bend, rather than the elliptical duct.

[William Rose]

@Sam Chak, nice photo! THank you! Your comme nt means a lot, considering how many nice solutions and explanations you have posted.

Answer 2

[Matt J]

Along the lines of what @Mathieu NOE commented, I think it does make more sense to start with an elliptic cylinder and clip off the top and bottom with cutting planes. That will be much easier than approaching it as a solid of revolution:

%Volume parameters
N=1000;
ellipse=scale(nsidedpoly(N),[1.5,1]);
planeSlope=0.1; 
planeHeight=0.2;
cylLength=0.5;

%Get the outside vertices
X=ellipse.Vertices(:,1);
Y=ellipse.Vertices(:,2);
Z=linspace(-cylLength,+cylLength,N);

[X,Z]=ndgrid(X,Z);
Y=repmat(Y,1,N);

keep=Z<=planeSlope*X+planeHeight  & Z>=-planeSlope*X-planeHeight;

X=X(keep); Y=Y(keep); Z=Z(keep);

%Triangulate and compute volume
[K,volume]=convhulln([X,Y,Z]);
volume,

volume = 1.8849

%Display
trisurf(K, X, Y, Z, ...
    'FaceColor', 'cyan', 'EdgeColor', 'none');
axis equal
camlight
view(-15,15)

●

[Robert Demyanovich]

@Matt J: There are two goals that I have. The main goal is to verify that my equations will yield the correct pictorial representation of the volume and correctly calculate the actual volume itself. Secondly, I would eventually like to have a "video" that shows a trapezoid slice revolving around the vertical line through "S" so that I can visually demonstrate my thinking of how the volume should be determined.

Visually, your chart appears to be correct. I see that you were able to determine the volume of the solid produced from your code. Can MatLab also measure the volume of the solid produced by rotating the trapezoid around the z axis (i.e. by not requiring the volume equation shown in my original post)? That would be a great help in verifying that the equations are, in fact, calculating the correct volume (or perhaps more accurately, that there is a consistency between the two). I would check this further by varying the angle θ, which then varies the constants that I listed in my original post.

[Matt J]

Can MatLab also measure the volume of the solid produced by rotating the trapezoid around the z axis (i.e. by not requiring the volume equation shown in my original post)?

That is what my code has already done (in numerical approximation), assuming you agree that the solid I have graphed matches what you are looking for.

As for the rest of what you are trying to do, I still think you are barking up the wrong tree by analyzing this as a solid of revolution, including your analytical volume computation. The volume can be decomposed into 3 parts: a central elliptic cylinder section, whose volume is trivially and two identical upper and lower sections with triangular x-z profiles (see below). The upper and lower sections should be fairly easy triple integrals in Cartesian coordinates. Note from symmetry that you really only need to integrate over one half (illustrated below) of one section

%keep=Z<=planeSlope*X+planeHeight  & Z>=planeSlope*min(X)+planeHeight & Y>=0;

[Robert Demyanovich]

Although I don't know how to create a surface/volume of revolution in MatLab, I'm surprised that this seems to be so difficult. There are many videos explaing how a volume can be determined by revolving an area around an axis; but these videos do not reference Matlab. Still, since this is a tried and true technique for creating volumes, I guess I'm surprised that this hasn't been implemented or is very difficult to implement in MatLab.

Further, the creation of a video showing the trapezoid slice revolving around the axis is rather important. It doesn't seem to me that the method you propose would provide for that. However, your method is certainly viable if all I want is to produce a solid.

[Matt J]

There are many videos explaing how a volume can be determined by revolving an area around an axis;

But when the area is changing with angle? That seems non-standard.

[Matt J]

To make the movie, you could use this File Exchange tool,

https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection?s_tid=srchtitle

which computes the intersection of two triangulated surfaces. To triangulate your solid you would give inputs from my code,

TR = triangulation(K,X,Y,Z);

while the second triangulated surface would just be the cut plane you want for a particular frame of your movie. You would represent the plane as a sufficiently big rectangle (two triangles) passing through S.

Answer 3

[Matt J]

Here's a little bit more of a polished version of my original answer, which also generates a movie of the trapezoidal cross-sections. You can adjust the VideoWriter properties to your liking.

This code requires the download of AxelRot,

https://www.mathworks.com/matlabcentral/fileexchange/30864-3d-rotation-about-shifted-axis

%Volume parameters (independent)
N=100; %discretize the perimeter into N points
diam=8; %tube diameter
ho=1;
hi=4;
xS=-3; %x-coordinate of S

%Volume parameters (derived)
dh=(hi-ho)/2;
angle=atand(dh/diam);
axMajor=hypot(diam,dh);
hS=ho+2*dh/diam*(xS+diam/2);

%Compute boundary vertices
t=linspace(0,360,N+1)'; t(end)=[];
V=0.5*[axMajor*cosd(t),diam*sind(t),0*t];

XYZ1=AxelRot( (V+[0,0,+ho/2])' ,+angle,[0,-1,0], [-diam/2,0,+ho/2])';
XYZ2=AxelRot( (V+[0,0,-ho/2])' ,-angle,[0,-1,0], [-diam/2,0,-ho/2])';
XYZ=[XYZ1;XYZ2];

%Triangulate and compute volume
[K,volume]=convhulln(XYZ);
volume,

volume = 122.1617

%Display
trisurf(K, XYZ(:,1), XYZ(:,2), XYZ(:,3), ...
    'FaceColor', 'cyan', 'EdgeColor', 'none','FaceAlpha',0.5);
axis equal; xlabel x; ylabel y; zlabel z
camlight
view(-30,20)

%Movie
shg
v=VideoWriter('trapezoidMovie.mp4','MPEG-4'); 
v.Quality=95;
v.FrameRate=30; 
v.open;
for i=1:N
    
    A=XYZ1(i,:);
    B=XYZ2(i,:);
    C=[xS,0,-hS/2];
    D=[xS,0,+hS/2];

    args=num2cell([A;B;C;D;A],1);
    [x,y,z]=deal(args{:});

    h=line(x,y,z,Color='r',LineWidth=2); 

    drawnow;
    v.writeVideo(getframe(gcf));

    if i<N, delete(h); end

end

v.close

[Robert Demyanovich]

@Matt J: This looks very interesting. I copied your code and put it into Matlab. Then I downloaded AxelRot and added that code at the end, which is what I think you are supposed to do with Matlab script (version 2021). But I'm getting the error "Function definitions are not supported in this context. Functions can only created as local or nested function in code files."

[Matt J]

which is what I think you are supposed to do with Matlab script (version 2021).

No, you were meant to leave AxelRot in its own file, just as you downloaded it. You need to put the file in a folder that is on your Matlab path, and then you can call it from anywhere.

[Matt J]

I've modified the example to include an export of the movie frames to an .mp4 video file.

[Robert Demyanovich]

I opened AxelRot.m and your code that contains the call to AxelRot. Matlab asked me if I wanted to add both to the path and I answered Yes. Still get an error when your code calls function varargout=AxelRot(varargin)

[Matt J]

I don't know what you did. What are the folder paths where you have put AxelRot.m and the other code? Did you look at your Matlab search path to see that those folders are on it? Does executing this at the command line give you anything?

which AxelRot

If all else fails, putting all code files that you want to use in the current Matlab folder is enough for Matlab to see everything.

[Robert Demyanovich]

I used the Open menu item and Open submenu command to open the Matlab code, let's call it Solution1, from let's say Folder1. I can clearly see that Axelrot is in Folder1 along with Solution1. But when I run Solution1, I still get the following error:

"Function definition are not supported in this context. Functions can only be created as local or nested functions in code files."

I am using MATLAB R2021b.

[Robert Demyanovich]

@Matt J: Got your code to work. Stupid error on my part. The only thing is I had to comment out a single line which was generating an error. The line is

v.Quality = 95

1)How do I adjust the speed with which the red trapezoid rotates? I though it would be on the line

v.FrameRate=30

but changing the value doesn't seem to do anything.

2) Is the video stored somewhere on my computer? I thought it would be by default in the Folder that contains the codes, but it isn't there.

[Matt J]

I though it would be on the line

Yes, that sets the playback speed for the video file (not the speed at which Matlab updates the figure). If you want Matlab to loop more slowly, you can stick a pause command at the bottom of the loop.

I thought it would be by default in the Folder that contains the codes, but it isn't there.

It should be.

[Matt J]

Sorry, it may be because you are writing to a JPEG 2000 video. The settings should have been,

v=VideoWriter('trapezoidMovie.mp4','MPEG-4'); 
v.Quality=95;

[Robert Demyanovich]

Thanks for that, but I already changed it to 'MPEG-4'. v.Qaulity=95 does not now generate an error.

[Robert Demyanovich]

@Matt J: The quality of the video is poor despite my setting v.Quality=100. I have been looking online on how to increase the quality and the consensus seems to be that I need to create the figure (which ultimately becomes a frame) to 720p. But your code doesn't look anything like the examples so I am loathe to mess with your code (mostly because I don't fully understand it). Any suggestions on how I can increase the resolution of the video besides v.Quality?

Also, I would like to accept both your answer and William Rose's answer since each better addresses certain parts of my problem. Is that possible, or do I have to choose one or the other?

Thanks.

[Matt J]

You can't accept two, but don't worry. Just accept one and up-vote the other.

To improve video quality, use a maximized figure window when recording the video, and increase the LineWidth, e.g.,

h=line(x,y,z,Color='r',LineWidth=3); 

There are programmatic ways to maximize the figure window, but it might be easiest if you just do it manually.

[Robert Demyanovich]

To improve the video quality, I added pause(5) just ahead of the loop that adds the frames. This gives me 5 seconds to click the maximize button on the figure that is shown before the "animation" starts. I think that the frames that are added, subsequently, are at the higher resolution (at least it looks that way).

[William Rose]

I always learn a lot from @Matt J's answers.

I'm curious: do you have the same video quality issue with the code I posted? I realize it makes a different movie.

Regarding speed of rotation in @Matt J's code: His code saves one angular step per frame. So you can increase N in his code: then it will take a smaller angular step each time. Or you can decrease v.FrameRate from the usual 30 (fps) to something slower. Adding a pause statement will make is slower when you run the code and see it in Matlab, but it will not change how the video plays back, because the video is always 1 angular step per frame.

[Robert Demyanovich]

@William Rose: Matt J's answer is excellent for creating the video that I was looking for. However, the actual image generated is not quite correct (but it doesn't really matter for the video). The equations that I presented indicate that to the rear of the stagnation point (S), the height does not seem to indicate a linear change (if viewed in the x-z plane). There is a steeper dropoff near the edge of the ellipse at 180 degrees for example. This is not meant to be a criticism of Matt J's work on this problem, as he was working from the original blue figures that I posted (Figs 1 and 2). I really appreciate the help that both of you have given me.

I put pause(5) before the loop in your movie and then maiximized the figure. The resulting video has a size that is almost 3 times larger and appears to be a bit crisper even though the font sizes are smaller.

The only issue I still have is resolving why the equations seem to be calculating a larger volume than Matlab calculates.

[Matt J]

The equations that I presented indicate that to the rear of the stagnation point (S), the height does not seem to indicate a linear change (if viewed in the x-z plane).

So, the volume we are talking about is no longer meant to have a purely trapezoidal cross-section in the x-z plane, like you had here?

I don't see how the change can be "nonlinear" in the x-z plane to the rear of S, if every cut through S is supposed to be a trapezoid.

As for the video generation, I don't see anything in my video-generation code that wasn't also in @William Rose's.

[Robert Demyanovich]

If you look at the x-z plot of William Rose's code shown below you can see that where the arrows point it looks like the height is falling faster. Actually looking at the entire green or blue line, it's likely that the slope is not linear throughout in this view. For a video I think it is fine to assume linear. The derivation of my equations simply say that a trapezoid can be assume at any specific azimuthal angle. The area of that trapezoid varies based on values of and . I don't think that my equation forces a view in the x-z plane to show a linear reduction in the top and bottom planes.

[William Rose]

@Robert Demyanovich,

[Edit: add a screenshot which I meant to include earlier.]

I see what you mean about the shape in the x-z plane. When I use the 3D rotate tool to rotate the 3D figure that my code makes, to view it looking down the y axis (az=0, el=0), I get this view:

This view makes clear something I had not realized, but which you obviously did realize: the "top" and "bottom" sides of the volume swept out by the trapezoid are not planes.

Non-planarity of the top and bottom also explains the persistent difference of about 5% between the volume of the convex hull and the volume obtained as the sum of the skinny slices. I should have realized this before. The figure below shows the trapezoids at 0, 45, and -45 degrees of rotation. The view is almost exactly down the y- axis, but not quite, because if it were exactly down the y-axis, the +45 and -45 trapezoids would exactly overlap. Notice how the top and bottom edges of the 0 degree trapezoid below the edges of the +,-45 degree trapezoids, in this view. That means the top surface dips down slightly, i.e. it is concave.

If you play around with the wireframe diagram, rotating with your mouse or touchpad in 3D, you can also appreciate the concave nature of the top and bottom surfaces. Here's a screenshot - but its easier to see the concavity when it's moving around, so try it on your machine, if you're skeptical.

The script I used to draw the overall wire frame and the convex hull and the figure with three trapezoids is attached: rotateTrapezoid.m. It does not save a video file.

Since the shape is concave in places, convhull() overestimates the volume. It also means that the volume obtained by taking planar slices through the top and bottom of the elliptical cylinder is not exactly the volume swept out by the rotating trapezoid.

[Matt J]

I don't think that my equation forces a view in the x-z plane to show a linear reduction in the top and bottom planes.

It does if that view is cross-sectional. Since the solid is created by revolving trapezoids, then the x-z plane cross section through the solid must be comprised of two trapezoids stitched together, one at and one at . I acknowledge that the two trapezoids do not have to combine to form a single larger trapezoid, but they must combine to form something that is at least piecewise linear.

[Robert Demyanovich]

@William Rose: Yes. I didn't know how convhull() works, but agree that it will overestimate. I attempted to estimate the volume by assuming planar slices through an elliptical cylinder (not sure I can call it that). Using this method, I estimate the volume as 1.18E-09, which is significantly different than what I have been calculating from my equation which is 2.09E-09. There may be issues with the numerical calculation of "p" or it's possible that my volume calculation using the equation is simply a factor of 2 too high (assuming that 1.18E-09 is a bit too high for reasons you state). I'll continue to work on this and report back.

[Robert Demyanovich]

@Matt J: Thank you for this comment. It would seem that you are correct in that if I am rotating trapezoids, then there can't be convex shapes in a cross sectional view despite being piecewise linear. I think my technique of using trapezoids will provide a good estimate, though, but I need to be clear that the trapezoids are estimated cross sections.

The equation for "p" is quite complicated. There is a unique condition where the equation for p is much simpler. This is the case where . I'm going to play around with this to see what the cross section looks like for this case and also what volume is calculated.

[Robert Demyanovich]

@Matt J and William Rose: As I suspected, the theoretical case where does in fact yield a cross section that is trapezoids throughout as shown in the figure below.

As such this was the basis of my initial drawings shown as Fig. 1 and Fig.2 in the original post. Now, in the non-theoretical case, w and B may not be exactly correct so this could potentially contribute to the reason why trapezoids are not seen throughout, or maybe it is really that way - I don't know. But I think I can use the idea of trapezoids to estimate the volume for the non-theoretical case.

[William Rose]

@Robert Demyanovich,

[Edit: add .mpeg file made by the script]

This script makes a series of figures similar to the nice ones by @John D'Errico, but the surrounding volume is the convex hull around all the trapezoids, instead of being a cylinder with an elliptical cross section, sliced by two planes. I realize that this convex hull is not perfect, as discussed in an earlier comment. But this convex hull is exactly correct along the non-planar top and bottom edges, which you can see by rotating it to az=0, el=0. Screenshots from the .mpeg file during playback are shown, at the 8th and final frames.

[William Rose]

@Robert Demyanovich, you wrote:

"Are you now of the mind that your original equation for dV is correct? Or are you just using it as a point of reference."

I think my equation for dV is correct,

,

to the extent that each dV has the same length, p, on each long sides. This is a pretty good approximation, when the wedges are skinny. In my code, I use the length of the trailing edge of each wedge to calculate dV, all the way around. This means I overestimate the volume of each dV, for phi=0 to 180, because the leading edge is shorter than the trailing edge. But from phi=180 to 360, I underestimate the volume of each wedge, because the leading edge is longer than the trailing edge. The errors in my dVs from phi=0 to 180 pretty much cancel out the errors from 180 to 360, and that is why vol=sum(dV) does not change (at 4 significant digits), when I vary N from 40 wedges (9 deg each) to 720 (0.5 deg each).

If you need each dV to be more accurate, use p=mean(leading & trailing edge lengths) in the formula for each dV, instead of using p=trailing edge length.

For one wedge, comprised of 6 vertices, I expect the volume from convhull() will be very close to the volume from the formula, if we use the mean value of p. Check this with the following commands, after running rotateTrapezoid, with N=80:

wedge1verts=[hst;hsb;pt([1,2],:);pb([1,2],:)];  % array with 6 vertices of wedge 1
[~,wedge1vol]=convhull(wedge1verts);            % wedge volume from convhull
p1=mean(p(1:2));                                % mean p for wedge 1
dV1=(B*p1^3/3+hs*p1^2/6)*dt;                    % wedge volume from formula
fprintf('convhull volume=%.3e, formula volume=%.3e.\n',wedge1vol,dV1)

convhull volume=1.052e-10, formula volume=1.053e-10.

The two volume estimates are very close, as expected.