%Partial-sphere volume calculation
% -Andy Ruina 24 September 2009
% To understand this whole file you have to understand
% several matlab commands that output cute things:
% sprintf, disp, datestr,num2st and now.
% The basic calculation uses:
% linspace, zeros, ones, and
% end (as an index, and also at end of for loop)
% Then the final graph uses
% plot, title, xlabel and ylabel.
%It helps debugging if you know your program started
% and stopped. Lines like the following, and the
% similar ones at the end, helps with this.
disp(sprintf('\n************************'))
disp(['Time at start of calculation: ' datestr(now)])
clear; %Clear all variables
% Set some basic parameters
R = 1; % Sphere radius
n = 1000001; % Number of cylindrical shells +1
%Parameterize the cylidrical shells by the
% angle between the the sphere symmetry plane
% and the top of the cylinder.
theta= linspace(0, pi/2, n);
%Arc-increments on the sphere:
deltatheta = (pi/2)/(n-1);
%The set of heights of the cylindrical shells:
h = R * sin(theta);
%The set of radii of the cylindrical shells:
r = R * cos(theta);
%The set of thicknesses of the cylindical shells:
deltar = sin(theta)*deltatheta;
%The set of volumes of the cylindrical shells:
deltaV = 2 * pi * r .* h .* deltar;
%The set of volumes of the solid cylinders
%inside:
Vcyl = pi * r.^2.*h;
%Place holders for the volumes
V=zeros(1,n);
Vtot = zeros(1,n);
Vhalfsphere = (2 * pi * R^3) / 3;
%This next is the `Reimann sum`, saving the partial sums
for i = 2:n
V(i) = V(i-1) + deltaV(i);
end
%Add the three parts
% 1) The half of a wedding ring calculated with the sum,
% 2) The solid cylinder in the middle of the ring
% 3) The bottom half of the sphere.
Vtot = V + Vhalfsphere + Vcyl;
%Total volume of sphere calculated numerically above
spheretotalnumerical = Vtot(end);
disp(['The volume is: ' num2str(spheretotalnumerical)])
%Error in total volume of sphere
error = spheretotalnumerical/(4*pi/3) - 1; %C
disp(['The fractional error is: ' num2str(error)])
%Back to the bouyancy problem where the volume of
% the sphere in oil is supposed to be 2/3 of the total
% volume of the sphere.
Voil = (2/3)* (4/3) * pi*R^3 * ones(n,1);
plot(h,Vtot,h,Voil)
title(['Volume in oil vs h: crossing pt is solution -' ...
datestr(now)])
xlabel('h = height from center that is in oil')
ylabel('V(h) = volume submerged')
axis([0 1 0 5])
disp(['Time at end of calculation: ' datestr(now)] )