%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)] )