Break this file up into the appropriate bits....

----------------------------------

clear
echo on
%In this example we demonstrate the use of the matlab routine quad.m to
%perform a 2-D integral.  We shall integrate the function 1/exp(x^2+y^2)
%over the elliptical disk described by:

%  (x/a)^2 + (y/b)^2 < 1

%where a and b are the major and minor axes of the bounding ellipse,
%respectively.
pause
%First, let's plot up this function over its domain.  We shall take a to be
%2 and b to be unity:
global apass bpass
a=2;
b=1;
apass=a;
bpass=b;
x=[-a:.05:a];
y=[-b:.05:b];
z=zeros(length(y),length(x));
for j=1:length(x)
 for i=1:length(y)
  if (x(j)^2/a^2+y(i)^2/b^2)<1
   z(i,j)=1/exp(x(j)^2+y(i)^2);
  end
  echo off
 end
end
echo on
figure(1)
mesh(x,y,z)
xlabel('x')
ylabel('y')
set(gca,'Fontsize',18)
pause
%As you can see from the plot, the function is a maximum in the center and
%decreases toward the edges of the ellipse.  This is clearer in the contour
%plot:
figure(2)
contour(x,y,z)
xlabel('x')
ylabel('y')
hold on
plot(x,(1-x.^2/a^2).^.5,'k')
plot(x,-(1-x.^2/a^2).^.5,'k')
%The roughness around the edges is due to the finite discretization of the
%bounding ellipse.
pause
%Now for the integral itself.  We shall do this as two nested integrals.
%The inner integral integrates over:

% -b*(1-(x/a)^2)^.5 < y < b*(1-(x/a)^2)^.5

%and the outer integral integrates over -a < x < a. Actually, because of
%symmetry, we can just integrate over the positive quadrant and multiply
%the result by four, this being four times faster!
pause
%In the main program we just integrate over the outer loop:
integ=4*quad('inner',0,a)
%where the inner loop is given in the attached function inner.m, which
%returns the integral over y for an array of input x's.
pause
%Note that the points are unevenly spaced:  most of the points are located
%near the end of the ellipsoidal region, where the function value is small
%and isn't changing much!  This is a consequence of the mapping we are
%doing using the above integration method!  An optimized 2-D gaussian
%quadrature routine would probably give us just as accurate result with
%only a few function evaluations!
hold off
pause
%We can also evaluate the integral using matlab's two-dimensional
%quadrature routine quad2d.m.  Matlab also has routines dblquad.m and
%triplequad.m, but these are limited only to rectangles and rectangular
%prisms.  quad2d is a bit more sophisticated, in that it will allow y to be
%a function of x, such as is required here.  Both the 

ymax = @(x) b*(1-x.^2/a^2).^.5;

%We also generate figure so we can see where quad2d is putting the points:
figure(3)
contour(x,y,z)
xlabel('x')
ylabel('y')
hold on
plot(x,(1-x.^2/a^2).^.5,'k')
plot(x,-(1-x.^2/a^2).^.5,'k')

pause
%OK, now to calculate the integral:
fun = @(x,y) funin2(x,y); %This is the same function, with graphics added

integ2=4*quad2d(fun,0,a,0,ymax)
%Which is the same value...

%Note that the points are placed a bit differently, and a bit more
%uniformly...
hold off

echo off


------------------------------------------------------------

function out=inner(x)
% This function evaluates the integral over y for an array of values of x.
% Because this is the inner integral, we set the tolerance to be tighter
% than the default tolerance for the outer integral.

global apass bpass xpass
%We use the variable xpass to pass the values of x to the function we are
%integrating.

n=length(x);
out=zeros(size(x));
a=apass;
b=bpass;

for i=1:n
 xpass=x(i); %We evaluate the inner integral one x at a time.
 ylim=b*(1-xpass^2/a^2)^.5;
 out(i)=quad('funin',0,ylim,1e-7);
end
drawnow %We draw the figure showing where the function calls are.

end


--------------------------------------------------------------

function out=funin(y)
global xpass
out=1./exp(y.^2+xpass^2);
%Comment this next bit of graphics out
%if you want the program to run faster.
for i=1:length(y)
 plot(xpass,y(i),'o')
end

end


--------------------------------------------------------------

function out=funin2(x,y)
out=1./exp(y.^2+x.^2);
plot(x,y,'o')
end