Plotting Free Surface Using MATLAB | Get Results Quickly

[hunt_mat]

Hi,

I like to think that I can generally get MATLAB to do what I want and when I want but this has got me stumped. I want to plot a free surface defined by a double integral:
$$\eta =\frac{1}{4\pi^{2}}\int_{\mathbb{R}^{2}}\frac{\mu e^{-\mu^{2}/4}e^{i(kx+ly)}\tanh\mu}{U^{2}k^{2}-\mu (B-E_{b}\mu+\mu^{2})\tanh\mu}dkdl$$
Where $\mu=\sqrt{k^{2}+l^{2}}$. I wrote a routine that does a double integral trapezium rule reasonably well but I need to get it working for the integrand above. Is there a quick method I can use to do this?

I should add that U is chosen such that the denominator has no zeros.

 

[hunt_mat]

I should add, the algorithm I used to compute the 2D trapezium rule is
$$\begin{array}{rcl}<br /> \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy & = & \left(\sum_{i=1}^{N}\sum_{j=1}^{M}f(x_{i},y_{j})-\frac{1}{2}\sum_{i=1}^{N}f(x_{i},c)-\frac{1}{2}\sum_{i=1}^{N}f(x_{i},d)\right)\delta x\delta y \\<br /> & - & \left(\frac{1}{2}\sum_{j=1}^{M}f(a,y_{j})+\frac{1}{2}\sum_{j=1}^{M}f(b,y_{j})\right)\delta x\delta y -{} \\<br /> & - & \frac{1}{4}(f(a,c)+f(a,d)+f(b,c)+f(b,d))\delta x\delta y<br /> \end{array}$$
The piece of Matlab code I used to compute the integral is given in here. I tested it out by computing the inverve 2D Fourier transform of a Gaussian and that seemed to work fine.
function y=trap_2d(A,dx,dy)
N=length(A(:,1));
M=length(A(1,;

a=A(1,1)+A(1,M)+A(N,1)+A(N,M);
b=sum(A(1,+sum(A(N,;
c=sum(A(:,1))+sum(A(:,M));
u=zeros(1,N);
for i=1:N
u(i)=sum(A(i,;
end

d=sum(u);

y=(d-0.5*c-0.5*b-0.25*a)*dx*dy;