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