[V,Rc,VRR,I]=volprop(pdat,sdat,t1,t2,titl,doplot)

function [V,Rc,VRR,I]=volprop(pdat,sdat,t1,t2,titl,doplot)
% [V,Rc,VRR,I]=volprop(pdat,sdat,t1,t2,titl,doplot)
% This function has a structure similar to function areaprop. 
% It analyzes a volume produced by rotating an area in the xz
% plane about the z-axis through through angles t1 to t2 
% (measured in degrees). See areaprop for a more detailed
% description of the input variables pdat and sdat. Results
% computed are the volume, centroidal radius, second moments 
% of volume, and the inertia tensor corresponding to a cross-
% section defined using polyline and spline segments. 
% pdat,sdat - cell variables defining the polyline and
%             spline boundary curves
% t1,t2     - the volume is generated by rotating an area
%             in the xz-plane from angle t1 to angle t2
%             counterclockwise about the z axis. t1 and 
%             t2 are measured in degrees
% titl      - a character variable used as a problem title.
%             Input [] for titl if no title is desired.
% doplot    - doplot = -1 for no print of results or plots. 
%             doplot =  0 to print results without a plot.
%             doplot =  1 prints both output and a surface
%                       for the volume of revolution.
%             doplot =  2 prints the computed output, the
%                         planar cross-section plot and 
%                         the volume of revolution.
% V,Rc      - the volume and centroidal radius  
% VRR       - moment integral([x;y;z]*[x,y,z]*dxdydz)
% I         - inertia tensor for a body of unit mass 
%             density defined by
%             I = eye(3,3)*sum(diag(VRR))-VRR

%             HBW, 9/16/09
if ~exist('doplot','var'), doplot=0; end  
if nargin==0
  pdat=@volbook; sdat={}; t1=-90; t2=180; doplot=1;
  titl='Example from Adv. Math. & Mech. Appl. page 164'; 
end
if isa(pdat,'function_handle'), [pdat,sdat]=feval(pdat); end    
if ~isempty(pdat), pdat=chkdat(pdat); end
if ~isempty(sdat), sdat=chkdat(sdat); end
if doplot==2, plotdat(pdat,sdat,titl), pause, end
if doplot>0,  plotvol(pdat,sdat,t1,t2,titl), shg, end 
t1=pi/180*t1; s1=sin(t1); c1=cos(t1);
t2=pi/180*t2; s2=sin(t2); c2=cos(t2);
f2=[(t2+s2*c2)/2,  s2^2/2,      s2;
     s2^2/2,    (t2-s2*c2)/2,  -c2;
     s2,           -c2,         t2];
f1=[(t1+s1*c1)/2,  s1^2/2,       s1;
    s1^2/2,     (t1-s1*c1)/2,   -c1;
    s1,            -c1,          t1]; 
V10=0; V20=0; V11=0; V30=0; V21=0; V12=0;
np=length(pdat); ns=length(sdat);
if ~isempty(pdat)
  for k=1:np
    u=pdat{k}; xd=u(1,:); yd=u(2,:); 
    [v10,v20,v11,v30,v21,v12]=splpol(2,xd,yd); 
    V10=V10+v10; V20=V20+v20; V11=V11+v11;
    V30=V30+v30; V21=V21+v21; V12=V12+v12;
  end
end
if ~isempty(sdat)
  for k=1:ns
    u=sdat{k}; xd=u(1,:); yd=u(2,:);  
    [v10,v20,v11,v30,v21,v12]=splpol(2,xd,yd);
    V10=V10+v10; V20=V20+v20; V11=V11+v11;
    V30=V30+v30; V21=V21+v21; V12=V12+v12;
  end
end
V=(t2-t1)*V10;
VR=[s2-s1;-c2+c1;t2-t1].*[V20;V20;V11]; Rc=VR/V;
VRR=(f2-f1).*[V30,V30,V21;V30,V30,V21;V21,V21,V12];
I=eye(3,3)*sum(diag(VRR))-VRR;
if doplot>-1
  disp(' '), disp(titl)   
  outputv(V,Rc,VRR,I)
end

function outputv(V,Rc,VRR,I)
%disp(' '), disp('VOLPROP RESULTS'), disp(' ')
disp('R  = [x, y, z], dV = dxdydz'), disp(' ')
disp(['V   = Integral(dV) = ',num2str(V)]), disp(' ')
disp('Rc  = Integral(R*dV) / V ='), disp(Rc')
disp('VRR = Integral( R''*R*dV) ='), disp(VRR)
disp('I = eye(3,3)*sum(diag(VRR))-VRR ='), disp(I)

function [V10,V20,V11,V30,V21,V12]=splpol(type,xd,yd,doplot)
% [V10,V20,V11,V30,V21,V12]=splpol(type,xd,yd,doplot)
% Inertial properties contributions of a cubic spline or 
% polyline boundary part are computed. If the first and 
% last data points match,then a closed boundary is defined.
if nargin==0 % example ellipse geometry
  th=linspace(0,2*pi,41); xd=1+cos(th); yd=1+.5*sin(th); 
  type=1; doplot=1;
end 
if exist('doplot','var')
  if type ==1
    nd=length(xd); zd=xd(:)+i*yd(:); 
    z=spline(1:nd,xd(:)+i*yd(:),linspace(1,nd,81));
    plot(real(z),imag(z),'k',xd,yd,'.r')
  else
    plot(xd,yd,'k')
  end
  axis equal, xlabel('x axis'), ylabel('y axis')
  title('PLOT OF THE BOUNDARY DATA'), shg, pause
end
V10=xymom(type,xd,yd,1,0); V20=xymom(type,xd,yd,2,0);
V11=xymom(type,xd,yd,1,1); V30=xymom(type,xd,yd,3,0);
V21=xymom(type,xd,yd,2,1); V12=xymom(type,xd,yd,1,2); 

function v=xymom(type,xd,yd,n,m)
% v=xymom(type,xd,yd,n,m) integrates x^n*y^m*dx*dy 
% for an area bounded by either a cubic spline or a 
% polygon defined by data points xd,yd traversing
% the boundary in the a counter-clockwise sense.
% xd,yd - data points on the boundary
% n,m   - exponents for x^n*y^m
% type  - 1 for a spline or 2 for a polyline
% v     - value of the integral
zd=xd(:)+i*yd(:); nd=length(zd); td=1:nd; 
if type==1 | nargin< 5 % use a spline curve 
  v=quadgk(@(t)intgrn(td,zd,n,m,t),td(1),td(end));   
else                   % use a polyline
  f=@(z)real(z).^n.*imag(z).^m.*conj(z);
  v=quadgk(f,zd(1),zd(end),'waypoints',zd(2:end-1)); 
  v=imag(v);
end
v=v/(n+m+2); 

function v=intgrn(td,zd,n,m,t)
% v=intgrn(td,zd,n,m,t)
z=spline(td,zd,t); zp=splder(td,zd,t);
v=real(z).^n.*imag(z).^m.*imag(conj(z).*zp); 

function v=splder(td,zd,t)
% v=splder(td,zd,t) cubic spline differentiation
%    HBW, 9/6/09
n=length(td);
if n>3
  [b,c]=unmkpp(spline(td,zd));
  c=[3*c(:,1),2*c(:,2),c(:,3)];
  v=ppval(mkpp(b,c),t);
elseif n==3
  c=[[1;1;1],td(:),td(:).^2]\zd(:);
  v=c(2)+2*c(3)*t;    
elseif n==2   
  v=(zd(2)-zd(1))/(td(2)-td(1))*ones(size(t));   
else
 v=nan;  
 error('More than 1 data point is required.')
end

function u=chkdat(v)
u=v; f=@(z)[real(z(:))';imag(z(:))'];
if iscell(v)
  n=length(v); 
  for k=1:n
    if ~isreal(u{k}), u(k)={f(u{k})}; end
  end
else
  if all(imag(u)==0)
    u={[u(1,:);u(2,:)]};
  else
    u={f(v)};
  end
end

function plotdat(pdat,sdat,titl)
% plotdat(pdat,sdat,titl) plots a geometry
% composed of splines and polygon lines
if nargin<3, titl='BOUNDARY DATA PLOT'; end
ns=length(sdat); np=length(pdat); x=[]; y=[];
hold off
close
if ~isempty(sdat)
  for j=1:ns
    u=sdat{j};
    n=size(u,2); tp=linspace(1,n,80);
    zp=spline(1:n,u(1,:)+i*u(2,:),tp);
    xp=real(zp); yp=imag(zp);
    x=[x,xp]; y=[y,yp];
    plot(xp,yp,'k'); hold on  
  end
end
if ~isempty(pdat)
  for j=1:np
    u=pdat{j}; x=[x,u(1,:)]; y=[y,u(2,:)];
    plot(u(1,:),u(2,:),'k'); 
    hold on 
  end
end
xlabel('x axis'), ylabel('y axis')
title(titl)
xmin=min(x); xmax=max(x);
ymin=min(y); ymax=max(y);
dx=0.1*(xmax-xmin); dy=0.1*(ymax-ymin);
range=[xmin-dx,xmax+dx,ymin-dy,ymax+dy];
axis(range), axis equal, hold off
shg

function plotvol(pdat,sdat,t1,t2,titl,nrot,nspl)
% plotvol(pdat,sdat,t1,t2,titl,nrot,nspl)
% plots a volume of revolution, with the cross-
% section composed of splines and polygon lines
if nargin<7, nspl=21; end, if nargin<6, nrot=30; end
if nargin<5, titl='VOLUME OF REVOLUTION'; end
if nargin<4, t1=-90; t2=270; end
if nargin==0, [pdat,sdat]=feval(@volbook); end
pdat=chkdat(pdat); sdat=chkdat(sdat);
t1=pi/180*t1; t2=pi/180*t2;
nrot=max(nrot,ceil(20*abs(t2-t1)/(2*pi)));
thr=linspace(t1,t2,nrot+1); cr=cos(thr); sr=sin(thr);
ns=length(sdat); np=length(pdat); x=[]; y=[];
hold off
close
big=1e300; xmin=big; ymin=big; zmin=big;
xmax=-big; ymax=-big; zmax=-big;
if ~isempty(sdat)
  for j=1:ns
    u=sdat{j};
    n=size(u,2); tp=linspace(1,n,nspl)';
    zp=spline(1:n,u(1,:)+i*u(2,:),tp); 
    X=real(zp)*cr; Y=real(zp)*sr; Z=imag(zp)*ones(1,nrot+1);
    xmin=min([xmin;X(:)]); xmax=max([xmax;X(:)]);
    ymin=min([ymin;Y(:)]); ymax=max([ymax;Y(:)]);
    zmin=min([zmin;Z(:)]); zmax=max([zmax;Z(:)]);  
    surf(X,Y,Z), colormap([1 1 0]), hold on    
  end
end
if ~isempty(pdat)
  for j=1:np
    u=pdat{j};  
    X=u(1,:)'*cr; Y=u(1,:)'*sr; Z=u(2,:)'*ones(1,nrot+1);
    xmin=min([xmin;X(:)]); xmax=max([xmax;X(:)]);
    ymin=min([ymin;Y(:)]); ymax=max([ymax;Y(:)]);
    zmin=min([zmin;Z(:)]); zmax=max([zmax;Z(:)]);
    surf(X,Y,Z), colormap([1 1 0]), hold on     
  end
end
xlabel('x axis'), ylabel('y axis'), zlabel('z axis')
title(titl) 
dx=0.1*(xmax-xmin); dy=0.1*(ymax-ymin); dz=0.1*(zmax-zmin); 
range=[xmin-dx,xmax+dx,ymin-dy,ymax+dy,zmin-dz,zmax+dz];
axis equal, axis(range), hold off
shg