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



14 Oct 2002 (Updated )

Companion Software (amamhlib)

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;
xd=real(Zd)+4; zd=imag(Zd);
th=[-pi/2,pi]; nth=31;

% Compute the geometrical properties
disp('  ')                   
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('      axx       axz       azz')

% 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;
xxd=real(Zd)+4; zzd=imag(Zd);
th=[-pi/2,pi]; Nth=121;

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]=...
% [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;
  xd=real(Zd)+4; zd=imag(Zd);
  th=[-pi/2,pi]; nth=31;

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

% 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; 

% Evaluate complex points and derivative values
% on the spline curve which is rotated to form
% the volume of revolution
% 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
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])];


function f=anglefun(t)
% f=anglefun computes multipliers involving 
% t, the rotation angle of the volume.
c=cos(t); s=sin(t);


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;
  xd=real(Zd)+4; zd=imag(Zd);
  th=[-pi/2,pi]; nth=31;
xd=xd(:); zd=zd(:); nd=length(xd);
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
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