No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

volrevol
function volrevol
%
% volrevol
% ~~~~~~~~
% This program determines geometrical properties
% for a solid generated by rotating a closed spline
% curve through an arbitrary angle about the z-axis.
% A detailed description of the geometry is given in
% function volrev.
%
% User m functions called: volrev srfv
%----------------------------------------------
%
% Data for a cross section consisting of the lower
% half of a circle plus a square capped by the 
% upper half of a smaller semicircle. The geometry
% is rotated through 270 degrees about the z-axis.

n1=9; t1=-pi:pi/n1:0; n2=6; t2=0:pi/n2:pi;
Zd=[0,exp(i*t1),1/2+i+exp(i*t2)/2,0];
xd=real(Zd)+4; zd=imag(Zd);
th=[-pi/2,pi]; nth=31;
kn=[1,2,n1+2,n1+3,n1+n2+3,n1+n2+4];


% Compute the geometrical properties
[v,rg,Irr,x,y,z,aprop]=volrev(...
                       xd,zd,kn,th,nth);
disp('  ')                   
disp('PROPERTIES OF A VOLUME OF REVOLUTION')
disp(' ')
disp('Results Using Function VOLREV')
disp(['Volume = ',num2str(v)]), %disp(' ')
disp(['Rg = [',num2str(rg(:)'),']']), %disp(' ')
disp('Inertia Tensor ='), disp(Irr), %disp(' ')
disp('Area Properties'), %disp(' ')
disp('     area     xcentr    zcentr.')
disp(aprop(1:3))
disp('      axx       axz       azz')
disp(aprop(4:6))

% Run a second case to generate a dense set of
% surface coordinates to check results using
% function srfv.

N1=61; T1=-pi:pi/N1:0; N2=41; T2=0:pi/N2:pi;
Zd=[0,exp(i*T1),1/2+i+exp(i*T2)/2,0];
xxd=real(Zd)+4; zzd=imag(Zd);
th=[-pi/2,pi]; Nth=121;
Kn=[1,2,N1+2,N1+3,N1+N2+3,N1+N2+4];  

[V,Rg,IRR,X,Y,Z]=volrev(...
                   xxd,zzd,Kn,th,Nth,1);
[vt,rct,vrrt]=srfv(X,Y,Z);
disp('Results Using Function SRFV')
disp(['Volume = ',num2str(vt)])
disp(['Rg = [',num2str(rct(:)'),']'])
disp('Inertia Tensor ='), disp(vrrt)

%=========================================

function [v,rg,Irr,X,Y,Z,aprop,xd,zd,kn]=...
                   volrev(xd,zd,kn,th,nth,noplot)
%
% [v,rg,Irr,X,Y,Z,aprop,xd,zd,kn]=...
%                  volrev(xd,zd,kn,th,nth,noplot)
%~~~~~~~~~~~~~~~~~~~~~~~~~

% This function computes geometrical properties
% for a volume of revolution resulting when a
% closed curve in the (x,z) plane is rotated,
% through given angular limits, about the z axis.
% The cross section of the volume is defined by
% a spline curve passed through data points 
% (xd,zd) in the same manner as was done in
% function areaprop for plane areas.

% xd,zd - data vectors defining the spline 
%         interpolated boundary, which is 
%         traversed in a counterclockwise
%         direction
% kn    - indices of any points where slope
%         discontinuity is allowed to turn
%         sharp corners
% p     - vector of volume properties containing
%         [v, xcg, ycg, zcg, vxx, vyy, vzz,...
%         vxy, vyz, vzx] where v is the volume,
%         (xcg,ycg,zcg) are coordinates of the
%         centroid, and the remaining properties
%         are volume integrals of the following
%         integrand:
%         [x.^, y.^2, z.^2, xy, yz, zx]*dxdyxz
% X,Y,Z - data arrays containing points on the
%         surface of revolution. Plotting these
%         points shows the solid volume with 
%         the ends left open. Function fill3 
%         is used to plot the surface with ends
%         closed
% aprop - a vector containing properties of the
%         area in the (x,z) plane which was used
%         to generate the volume. aprop=[area,...
%         xcentroidal, ycentroidal, axx, axz, azz]. 

% User m functions called: rotasurf, gcquad,
%         curve2d, anglefun, splined
%----------------------------------------------
if nargin==0 
  t1=-pi:pi/6:0; t2=0:pi/6:pi;
  Zd=[0,exp(i*t1),1/2+i+exp(i*t2)/2,0,-1];
  xd=real(Zd)+4; zd=imag(Zd);
  kn=[1,2,8,9,15,16];
  th=[-pi/2,pi]; nth=31;
end

% Plot a surface of revolution based on the 
% input data points
if nargin==6
  [X,Y,Z]=rotasurf(xd,zd,th,nth,1);
else
  [X,Y,Z]=rotasurf(xd,zd,th,nth); pause
end

% Obtain base points and weight factors for the
% composite Gauss formula of order seven used in 
% the numerical integration
nd=length(xd); nseg=nd-1; 
[dum,bp,wf]=gcquad([],1,nd,7,nseg);

% Evaluate complex points and derivative values
% on the spline curve which is rotated to form
% the volume of revolution
[u,uplot,up]=curve2d(xd,zd,kn,bp);
% plot(real(uplot),imag(uplot)), axis equal,shg
u=u(:); up=up(:); n=length(bp);
x=real(u); dx=real(up); z=imag(u);
dz=imag(up); da=x.*dz-z.*dx;

% Evaluate line integrals for area properties
p=[ones(n,1), x, z, x.^2, x.*z, z.^2, x.^3,...
        (x.^2).*z, x.*(z.^2)].*repmat(da,1,9);
p=(wf(:)'*p)./[2 3 3 4 4 4 5 5 5];

% Scale area properties by multipliers involving
% the rotation angle for the volume
f=anglefun(th(2))-anglefun(th(1));
v=f(1)*p(2); rg=f([2 3 1]).*p([4 4 5])/v;
vrr=[f([4 5 2]); f([5 6 3]); f([2 3 1])].*...
    [p([7 7 8]); p([7 7 8]); p([8 8 9])];
Irr=eye(3)*sum(diag(vrr))-vrr;
aprop=[p(1),p(2:3)/p(1),p(4:6)];

%=========================================

function f=anglefun(t)
% f=anglefun computes multipliers involving 
% t, the rotation angle of the volume.
c=cos(t); s=sin(t);
f=[t,s,-c,(t+c*s)/2,s*s/2,(t-c*s)/2];

%=========================================

function [x,y,z,xd,zd]=rotasurf(xd,zd,th,nth,noplot)
% [x,y,z,xd,zd]=rotasurf(xd,zd,th,nth,noplot)
% This function generates points on a surface of
% revolution generated by rotating an area in
% the (x,z) plane about the z-axis
%
% xd,yz  - coordinate data for the curve in the
%          (x,y) which forms the cross section
% th     - [ThetaMin,ThetaMax] defining limits of
%          rotation angle about the z-axis
% nth    - number of theta values used to generate
%          surface values
% noplot - option given any value if no plot is
%          desired. Otherwise omit this value.
% x,y,z  - arrays of points on the surface
%
% User m functions called: none
%----------------------------------------------

if nargin==0
  n1=9; t1=-pi:pi/n1:0; n2=6; t2=0:pi/n2:pi;
  Zd=[0,exp(i*t1),1/2+i+exp(i*t2)/2,0];
  xd=real(Zd)+4; zd=imag(Zd);
  th=[-pi/2,pi]; nth=31;
end
xd=xd(:); zd=zd(:); nd=length(xd);
t=linspace(th(1),th(2),nth);
x=xd*cos(t); y=xd*sin(t); z=repmat(zd,1,nth);
if nargin==5, return; end
close; surf(x,y,z), title('VOLUME OF REVOLUTION')
xlabel('x axis'), ylabel('y axis')
zlabel('z axis'), colormap([1 1 1]); hold on
fill3(x(:,1),y(:,1),z(:,1),'w')
fill3(x(:,end),y(:,end),z(:,end),'w')
axis equal, grid on, hold off, shg

%=========================================

% function [z,zplot,zp]=curve2d(xd,yd,kn,t)
% See Appendix B

%=========================================

% function [val,bp,wf]=gcquad(func,xlow,...
%                    xhigh,nquad,mparts,varargin)
% See Appendix B

%=========================================

% function range=cubrange(xyz,ovrsiz)
% See Appendix B

%=========================================

% function val=splined(xd,yd,x,if2)
% See Appendix B

%=========================================

% function [v,rc,vrr]=srfv(x,y,z)
% See Appendix B

Contact us