Plotting with MATLAB

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:
[tex]
\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
[/tex]
Where [itex]\mu=\sqrt{k^{2}+l^{2}}[/itex]. 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
[tex]
\begin{array}{rcl}
\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 \\
& - & \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 -{} \\
& - & \frac{1}{4}(f(a,c)+f(a,d)+f(b,c)+f(b,d))\delta x\delta y
\end{array}
[/tex]
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;