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