trapezoidal_rule_double_integral(x, y, mat)

function out = trapezoidal_rule_double_integral(x, y, mat)
%  
% Calculates a double integral using the trapezoidal rule.
% 
%
% By: Mohammed S. Al-Rawi
% IEETA, University of Aveiro
% Portugal
% March, 4, 2013 
%
% Please send feedback if you have any issue with this function to
%  [e-mail]
%
% -- mat is a 2D matrix that contains the function values, e.g. f(x,y),
% size of f is [M N]
% -- x and --y are the vectors to be used with the function to generate mat
% size of x is M, and size of y is N
%
%
% see the examples below for more illustration
%
% Example:

%  double_integral(e^(y-x)) x=0 to .5, y= 0 to .5
% % 
%  x=0:.01:.5; 
%  y=0:.02:.5;
%  [xx,yy] = meshgrid(x,y);
%  z=yy-xx;
%  mat = exp(z)'; % the transpose is needed since meshgrid returns the
%  values transposed
%  out = trapezoidal_rule_double_integral(x, y, mat)
% 
% %  out =
% 
%     0.255262565920498

%   Comparing to the function quad2d:
%  fun = @(x,y) exp(y-x); Q = quad2d(fun, 0,.5,0,.5)
% 
% Q =
% 
%    0.255251930412767
% 
%
% here is another more complex Example 
% fun = @(x,y) exp(sin(y)-x.^2); Q = quad2d(fun, 0,.5,0,.5)
% 
% Q =
% 
%    0.297471075439585
% 
% x=sort(rand(100,1)/2); %non-uniform sampling
% y=0:.01:.5;
% [xx,yy] = meshgrid(x,y);
%  z=yy-xx;
% mat= exp(sin(yy)-xx.^2)';
% out = trapezoidal_rule_double_integral(x, y, mat)
% 
% out =
% 
%    0.292243052795743
%
%
%
% But, why would one wanna use this function, having the accurate quad2d?  
% The reason is that quad2d only works with functions, with an upper and a  
%  lower limit while this function accepts as input the vectors x and y. 
% BTW, x and/or y could be non-uniformly spaced. So, this function could be 
% used to implement fast integration in some applications where the
% computation of f is extremely expensive. Of course the accuracy is  much
% lower than that of quad2d, but as I said, you may find it useful in some
% situations that require high speed, and you already have the vectors
%  x, y, and mat. As a better example, consider integrating empirical data.
%  To apply quad2d for empirical data, you will have to write a function
%  that interpoltes from mat(x,y) and pass it to quad2d, but, for
%  trpz...2d, this is not needed.

% N=length(y);
% tmp=zeros(N,1);
% for i=1:N
%     tmp(i)= trapz(x,mat(:,i));
% end
% out =trapz(y,tmp);
% 

% simple, ain't it? Enjoy khafla.


% Well, there is a simpler statement that uses only one line, suggested by 
% Richard Crozier http://www.mathworks.com/matlabcentral/fileexchange/authors/34660
% 
 out = trapz(y,trapz(x,mat,1),2);
 
 % Now it is even simpler thanks to Richard :)