# Does not evaluate the Integral and a plane plot

I am trying to compute a ContourPlot and a VectorPlot using the following functions:

k = 2 π/(632.8*10^(-9))
α = ArcSin[0.95]

a = Exp[-(Tan[θ]^2/Tan[α]^2)] BesselJ[1, Tan[θ]/ Tan[α]]


a is to apply into the functions ex and ey. These are the electric field components.

ex = Compile[{{x, _Real}, {y, _Real}},
NIntegrate[
NIntegrate[
a Sqrt[Cos[θ]]
Sin[θ] Exp[
I k ( x Sin[θ] Cos[ϕ] +
y Sin[θ] Sin[ϕ] +
z Cos[θ])] (-Sin[ϕ]), {ϕ, 0, 2 π},
PrecisionGoal -> 5, MaxRecursion -> 20], {θ,
0, α}, PrecisionGoal -> 5,
MaxRecursion -> 20] Assumptions -> {y > 0}]

ey = Compile[{{x, _Real}, {y, _Real}},
NIntegrate[
NIntegrate[
a Sqrt[Cos[θ]]
Sin[θ] Exp[
I k ( x Sin[θ] Cos[ϕ] +
y Sin[θ] Sin[ϕ] +
z Cos[θ])] (Cos[ϕ]), {ϕ, 0, 2 π},
PrecisionGoal -> 5, MaxRecursion -> 20], {θ,
0, α}, PrecisionGoal -> 5,
MaxRecursion -> 20] Assumptions :> Element[x, Reals]]


If I try to set values into the function it only gives me back the function itself.Then I tried to evaluate the integrals with Assumptions. It does not work and I don't know how to apply this to 2 variables.

{ex, ey} /. {x -> 2, y -> 1}


Plotting ex and ey with ContourPlot does not give me anything either. (if it is done right, there are round contours.)

I used bbgodfreys suggestions and plotted a ContourPlot. The Outcome wasn't expected.

ContourPlot[
Re[ex[r, -10, 10] + Re[ey[r, -10, 10]]], {r, 0, 10}, {\[Phi]0, 0,
10}]


and got:

• 1) As far as I know, NIntegrate does not take Assumptions options. 2) you have syntax errors in the positioning of your Assumptions statements as well. 3) There is a lot going on in your code: have you tried constructing a much simpler example and getting that to work first? – MarcoB Aug 14 '15 at 20:00
• Your usage of Compile is wrong. Work thorough some of the examples in the docs to see how it works. (usage will be like ex[1,2] not ex/. x->... ). There is of course no point to specifying Assumptions on x,y since they always will have distinct values inside the function. – george2079 Aug 14 '15 at 20:17
• I tried it first without Assumptions and without the a. But it was the same. I tried it without Compile, too. – donut Aug 14 '15 at 20:26
• on further thought I don't see any point to using Compile here in the first place. – george2079 Aug 14 '15 at 20:26
• I tried ey[0, 0] but it took so long so I canceled it. Plus it says: NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. – donut Aug 14 '15 at 20:29

This problem can be simplified enormously by performing the integral over ϕ symbolically. To do so, first consider the ϕ-dependent exponential in the integrand.

Exp[I k (x Sin[θ] Cos[ϕ] + y Sin[θ] Sin[ϕ] + z Cos[θ])]


It can be transformed and simplified as follows.

expt = Exp[I k ( x Sin[θ] Cos[ϕ] + y Sin[θ] Sin[ϕ] + z Cos[θ])] /.
{x -> r Cos[ϕ0], y -> r Sin[ϕ0]} // Simplify
(* E^(I k (z Cos[θ] + r Cos[ϕ - ϕ0] Sin[θ])) *)


Now, perform the ϕ integrals

exϕ = Integrate[expt (-Sin[ϕ]), {ϕ, 0, 2 π}][[1]]
(* -2 I E^(I k z Cos[θ]) π BesselJ[1, k r Sin[θ]] Sin[ϕ0] *)

eyϕ = Integrate[expt Cos[ϕ], {ϕ, 0, 2 π}][[1]]
(* 2 I E^(I k z Cos[θ]) π BesselJ[1, k r Sin[θ]] Cos[ϕ0] *)


These results next are copied into the θ integrals

ex[r_, ϕ0_, z_] := NIntegrate[a Sqrt[Cos[θ]] Sin[θ]
(-2 I E^(I k z Cos[θ]) π BesselJ[1, k r Sin[θ]] Sin[ϕ0]), {θ,  0, α}]
ey[r_, ϕ0_, z_] := NIntegrate[a Sqrt[Cos[θ]] Sin[θ]
(2 I E^(I k z Cos[θ]) π BesselJ[1, k r Sin[θ]] Cos[ϕ0]), {θ,  0, α}]


A sample result is

Plot[Im[ey[r, π/4, 0]], {r, 0, 20/k}, AxesLabel -> {r, "Im[ey]"}]


Edit: Contour Plot

ContourPlots of the functions ex and ey can be generated in a straightforward manner, although the process is a bit slow. To modestly improve speed and simplify coding, rescale r and z by k, and factor -2 π I Sin[ϕ0] from ex and 2 π I Cos[ϕ0] from ey to obtain

e[r_, z_] := NIntegrate[a Sqrt[Cos[θ]] Sin[θ] Exp[I z Cos[θ]] BesselJ[1, r Sin[θ]],
{θ,  0, α}, PrecisionGoal -> ∞, AccuracyGoal -> 5]


Note that AccuracyGoal is reduced as well. I have not experimented with Method, which may further improve speed. Finally, if only Re[e] is desired, use

eRe[r_, z_] := NIntegrate[a Sqrt[Cos[θ]] Sin[θ] Cos[z Cos[θ]] BesselJ[1, r Sin[θ]],
{θ,  0, α}, PrecisionGoal -> ∞, AccuracyGoal -> 5]
ContourPlot[eRe[r, z], {r, 0, 10}, {z, 0, 10},
FrameLabel -> {"k r", "k z"}, ContourLabels -> True]


The imaginary part can be obtained by replacing Cos[z Cos[θ]] by Sin[z Cos[θ]].

• Thank you very much, it works much more quicker now. What do I do if I want to do a ContourPlot with ex+ey? I can not plot it over x and y axis anymore right? But over the r and ϕ0? – donut Aug 15 '15 at 10:26
• s0 = ex Conjugate[ex] + ey Conjugate[ey]; s1 = ex Conjugate[ex] - ey Conjugate[ey]; s2 = ex Conjugate[ey] + ey Conjugate[ex]; s3 = I (ex Conjugate[ey] + ey Conjugate[ex]); like this? VectorPlot[{s1[r, z], s2[r, z]}, {r, -10, 10}, {z, -10, 10}] – donut Aug 15 '15 at 14:12
• @donut Yes, although it is easier to work with -2 π I Sin[ϕ0] e[r, z] and 2 π I Cos[ϕ0] e[r, z] instead of ex[r, ϕ0, z] and ey[r, ϕ0, z]. Also, I would use {r, 0, 10}, {z, 0, 10}. because r has a minimum value of zero, and the s are symmetric or antisymmetric about z = 0. If my answer meets your need, please do not forget to accept it. Best wishes. – bbgodfrey Aug 15 '15 at 14:47