4
$\begingroup$

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:

enter image description here

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – MarcoB Aug 14 '15 at 20:00
  • $\begingroup$ 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. $\endgroup$ – george2079 Aug 14 '15 at 20:17
  • $\begingroup$ I tried it first without Assumptions and without the a. But it was the same. I tried it without Compile, too. $\endgroup$ – donut Aug 14 '15 at 20:26
  • $\begingroup$ on further thought I don't see any point to using Compile here in the first place. $\endgroup$ – george2079 Aug 14 '15 at 20:26
  • $\begingroup$ 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. $\endgroup$ – donut Aug 14 '15 at 20:29
6
$\begingroup$

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]"}]

enter image description here

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]

enter image description here

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – donut Aug 15 '15 at 10:26
  • $\begingroup$ @donut Please see material added to answer. I hope that this is helpful. $\endgroup$ – bbgodfrey Aug 15 '15 at 11:49
  • $\begingroup$ this helps! One last question, now I want to put ex and ey into my s0,1,2,3. I have to transform the e back right? $\endgroup$ – donut Aug 15 '15 at 14:10
  • $\begingroup$ 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}] $\endgroup$ – donut Aug 15 '15 at 14:12
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 15 '15 at 14:47

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.