How to plot a 2 variable function of an integral?

Question

I have this function:

$$E_p(x_0,y_0,t)=E_0e^{-i\omega t}\int^{\frac{a}2}_{\frac{-a}2}\int^{\frac{b}2}_{\frac{-b}2}\frac{e^{ik\sqrt{(x-x_0)^2+(y-y_0)^2+z^2}}}{\sqrt{(x-x_0)^2+(y-y_0)^2+z^2}}dxdy$$

Where $E_0,\omega,a,b,z$ are known constants.

Now I want to plot a new function defined as:

$$I_p(x_0,y_0)=c\epsilon_0\langle E_pE_p^*\rangle$$ Where $c,\epsilon_0$ are given constants and $\langle \rangle$ is the time average. So the term $e^{i\omega t}$ cancels and we can write the function as:

$$I_p(x_0,y_0)=c\epsilon_0E_0^2\int^{\frac{a}2}_{\frac{-a}2}\int^{\frac{b}2}_{\frac{-b}2}\frac{1}{[(x-x_0)^2+(y-y_0)^2+z^2]}dxdy$$

In an attempt to plot this function I wrote the following code in Mathematica:

z=10;
r[x0_,y0_,x_,y_]:=Sqrt[(x-x0)^2 + (y-y0)^2 + z^2];
E0=1; k=1.15*10^7;
a=10^-3; b=10^-3;
c=10^8; e0=10^-12;
Ip[x0_, y0_] := c*e0*E0^2*Integrate[1/(r[x0,y0,x,y])^2,{x,-a/2,a/2},{y,-b/2,b/2}];

But when I add the following line to get the value of $I_p$ in $(0,0)$

Ip[0,0]

I get an unsolved integral of x. I think also because of this the ContourPlot or Plot3D commands don't work properly too.

Answer 1

A few points:

1. in your definition of r you missed a _ after y0;
2. use NIntegrate rather than Integrate for numerical integration, such as in the definition of Ip;
3. protect the input to Ip so it is only evaluated with explicitly numerical input, especially if you intend to use it in plotting; use NumericQ to do so; see also this FAQ: User-defined functions, numerical approximation, and NumericQ.

Here is the modified code (the rest of the definitions can be taken from your own original code):

r[x0_, y0_, x_, y_] := Sqrt[(x - x0)^2 + (y - y0)^2 + z^2];

Ip[x0_?NumericQ, y0_?NumericQ] := 
  c*e0*E0^2*
   NIntegrate[1/(r[x0, y0, x, y])^2, {x, -a/2, a/2}, {y, -b/2, b/2}];

Plot3D[Ip[x, y], {x, -10, 10}, {y, -10, 10}]