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:


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

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)$


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

  • 1
    $\begingroup$ 1) in your definition of r you missed a _ after y0; 2) use NIntegrate rather than Integrate for the numerical integral in Ip`. $\endgroup$ – MarcoB Jun 18 '18 at 22:39
  • $\begingroup$ @MarcoB the first mistake just happened here. When using symbols from the basic math assistant how can I indicate NIntegrate? $\endgroup$ – Alireza Jun 18 '18 at 22:46
  • 1
    $\begingroup$ I am not sure that you can from the palettes. It is best to enter it manually. $\endgroup$ – MarcoB Jun 18 '18 at 22:48

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] := 
   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}]

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.