Hi everyone,

I am trying to calculate a 2d Integral which has another 2d integral in the integrands. here is the setup:

alpha = 1; sigma1 = 4; sigma2 = 1; phi = 0.5;

a = 0.5*(Cos[phi]/sigma1)^2 + 0.5*(Sin[phi]/sigma2)^2

b = 0.25*Sin[2*phi]*(-1*sigma1^-2 + sigma2^-2)

c = 0.5*(Sin[phi]/sigma1)^2 + 0.5*(Cos[phi]/sigma2)^2

A = alpha/(2*pi*sigma1*sigma2)

t[x1_, y1_, x_, y_] := A*Exp[-1*(a*(x - x1)^2 + b*(x - x1)*(y - y1) + c*(y - y1)^2)]

B2[x1_, y1_, x_, y_] := HeavisideLambda[x - x1, y - y1]

trans:= NIntegrate[B2[x1, y1, xz, yz]*NIntegrate[t[xz, yz, xp, yp]*(B2[x2, y2, xz, yz] - B2[x2, y2, xp, yp]), {xp, x2 - 10*sigma1, x2 + 10*sigma1}, {yp, y2 - 10*sigma2, y2 + 10*sigma2}], {xz, x1 - 1, x1 + 1}, {yz, y1 - 1, y1 + 1}]

trans[0,0,1,0]

trans[0,0,10,10]

trans[0,0,50,50]

I am trying to get numbers for trans[0,0,x,y]for different x,y, but I couldn't get any results yet. I tried analytic integration too (Integrate instead of NIntegrate) but it didn't help either.

I am quite new to Mathematica so I don't know what options in NIntegrate can help me here.

I have tried the same integral in Maple, there it works for some values of phi, and doesn't for others (so i'm sure the thing should some how work here too). for example, it works fine for phi=0, but for phi=0.5 it fails. For phi=Pi/4, it works but the results don't go to zero as x,y get far from zero (e.g. trans[0,0,50,50]), which is the expected behavior for this integral.

any ideas?

Thanks in advance