clear
echo on
%Solution to Problem 27
%
%In this problem we are asked to determine
%the moment of inertia of a prolate spheroid
%about the three pricipal axes of revolution.
%This integral is actually quite easy to evaluate
%analytically.  In this solution I will demonstrate
%several approaches to obtaining the solution.
pause
%First we do the analytical solution.  The interior
%of the spheroid is described by the inequality:
%
%  (z/a)^2 + (x/b)^2 + (y/b)^2 < 1
%
%This immediately suggests a mapping onto a sphere.
%This is easily accomplished by defining z*=z/a,
%x*=x/b, and y*=y/b.  Note that the volume differential
%must also be mapped, e.g., dV* = dV / (a b^2).
pause
%Next we define spherical coordinates in this mapped
%domain, which is itself another mapping.  Thus, we
%let z*=r* cos(theta), x*=r* sin(theta) cos(phi), and
%y*=r* sin(theta) sin(phi).  The limits of integration
%are now 0< r* <1, 0< theta< pi, 0< phi< 2pi, so the sphere
%has been mapped onto a rectangular prism!  The volume
%element dV* in this new system of variables becomes:
%dV* = r*^2 sin(theta) dtheta dphi dr*.
pause
%It turns out that if we plug everything back in, the 
%triple integral -separates- into the product of three 
%single integrals.  This results in a tremendous simplification
%because three single integrals are much faster to evaluate
%accurately than one triple integral!  If you plug it
%all back in and evaluate the resulting simple integrals
%analytically, you get that the inertia about the z-axis
%is given by:
%
%  Iz = a b^4 dens (8 pi)/15
%
%and that about the other two axes is given by:
%
%  Ix = Iy = a b^4 dens (4 pi)/15 (1 + a^2/b^2)
%
pause
%Now we turn to the numerical solution.  There are many
%ways we can do this.  The easiest way to code it up is
%to use Monte Carlo integration of the original function
%imbedded in a rectangular prism.  It is also the slowest!
%Let's try this approach first, but only to a limited 
%accuracy so that you won't spend forever waiting for it
%to finish running!  We shall put the imbedded function in
%the file fimbed.m
pause
n=10^4;
%Note that the limits of integration of the
%rectangular prism the function is imbedded in
%are -1< x<1, -1< y<1, -2< z<2.  Thus:
sum=zeros(1,3);
sumsq=sum;
for i=1:n
  x=rand*2-1;
  y=rand*2-1;
  z=2*(rand*2-1);
    f=fimbed(x,y,z);
    sum=sum+f;
    sumsq=sumsq+f.^2.;
    echo off
  end
end
echo on
integ=sum/n*16;
error=16*(sumsq/n-(sum/n).^2).^.5/n^.5;
pause
%Note that we have set up the problem so that
%we calculate all of the moments simultaneously.
%Let's compare the results to the exact results:
integexact=2*2.5*4*pi/15*[5,5,2];
[integ;error;integexact]'
%From which it is clear that the x and y moments
%are the same, and that the error in the Monte
%Carlo integration routine is around 1%.  To get
%the desired 0.1% error we require around 10^6
%points!  It is hard to get very high accuracy
%from Monte Carlo integration!
pause
%OK, now let's try doing this via three nested
%quadrature routines.  We will use the matlab routine
%myquad (from an earlier example) and we will only do 
%the z-moment part - it is simple to extend it to the other 
%two moments as well.  Also, rather than do it using an 
%imbedded function (which would lead to singularities) we shall 
%do it using the mapped function in spherical coordinates.
pause
%In the innermost loop we will integrate over phi, in
%the next we integrate over theta, and in the outer loop
%we integrate over r.  This will require the attached
%programs rfun, thetafun, and phifun.  Anyway, we have:
tolr=0.0005;
zmoment=myquad('rfun',0,1,tolr)
%which has an error of only
zerror=zmoment-integexact(3)
%which is much less than the desired tolerance.
pause
%It makes a whole lot of sense to do the above
%nested integration in the mapped domain, because that
%eliminates all of the singularities in the derivatives.
%If you used the imbedded domain approach, the integration
%takes practically forever.
pause
%As a final example, we will integrate the expression
%by integrating over chords in the spheroid.  We first
%integrate over x, then over y, and finally over z.  As
%in the example, the limits of integration over x are
%from -xlim to xlim, for y from -ylim to ylim, and for
%z from -a to a.  The quantities xlim and ylim can be
%obtained from:
%
%  ylim = b*(1 - (z/a)^2)^.5
%
%  xlim = b*(1 - (z/a)^2 - (y/b)^2)^.5
pause
%In this example we will use the attached routines to
%calculate the nested integrals.  We will again use myquad.m as
%the integration routine, and make use of symmetry to reduce 
%our computations by an order of magnitude.  We can do all three
%integrals with the same routines by setting a flag in the
%main script and in the innermost function 'xfun' too.
global a b flag
a=2;
b=1;
pause
moments=zeros(3,1);
for flag=1:3
 moments(flag)=8*myquad('zfun',0,a,.001);
end
pause
%Which yields the moments and the errors:
format short e
result = [moments,moments-integexact']
format short
%which is well within the desired tolerance, although not nearly as
%accurate or as fast as we had with the mapped domain results.  If 
%we had used the matlab routine 'quad' it would have taken longer
%still, and we would have gotten pages of recursion warnings!  In 
%part this is because of the difference between a closed (quad) and open 
%(myquad) quadrature rule:  The closed rule evaluates the function
%on the boundary, which in a nested integral such as this can cause 
%alot of trouble.
echo off



_______________________


function rf=rfun(r)
toltheta=0.00025;
global rt
rf=zeros(size(r));
for i=1:length(r)
 rt=r(i);
 rf(i)=myquad('thetafun',0,pi,toltheta);
end



_______________________

function thetaf=thetafun(theta)
tolphi=0.0001;
global thetat
thetaf=zeros(size(theta));
for i=1:length(theta)
 thetat=theta(i);
 thetaf(i)=myquad('phifun',0,2*pi,tolphi);
end



_______________________


function phif=phifun(phi)
global rt thetat
phif=2*2.5*rt^4*sin(thetat)^2*(cos(phi).^2+sin(phi).^2)*sin(thetat);


_______________________


function f=fimbed(x,y,z)
f=[0,0,0];
if (x^2+y^2+(z/2)^2)<1,
 f(1)=(y^2+z^2)*2.5;
 f(2)=(x^2+z^2)*2.5;
 f(3)=(x^2+y^2)*2.5;
end

_______________________


function out=zfun(z)
%This function returns the integral over y
%which is in turn the integrand for the
%integral over z.
global a b zpass
out=zeros(size(z));
for i=1:length(z)
 zpass=z(i);
 ylim=b*(1-(zpass/a)^2)^.5;
 out(i)=myquad('yfun',0,ylim,.0005);
end



_______________________


function out=yfun(y)
%This function returns the integral over x
%which is in turn the integrand for the
%integral over y.
global a b zpass ypass
out=zeros(size(y));
for i=1:length(y)
 ypass=y(i);
 xlim=b*(1 - (zpass/a)^2 - (ypass/b)^2)^.5;
 out(i)=myquad('xfun',0,xlim,.0002);
end



_____________________


function out=xfun(x)
%This function returns the integrand we
%are integrating.
global a b zpass ypass flag

if flag==1
 out=(zpass^2+ypass^2)*2.5*ones(size(x));
end
if flag==2
 out=(x.^2+zpass^2)*2.5;
end
if flag==3
 out=(x.^2+ypass^2)*2.5;
end


__________________________