3
$\begingroup$

I want to try plotting this: enter image description here

As you can see the first axis is v and the second is just dependent of the angle of $\Phi$ .

The function is

plot = ParametricPlot[
  141.20^2/(16 π) (Abs[Subscript[I, 0]]^2 + 
 4 Abs[Subscript[I, 1]]^2 *Cos[Φ]^2 + 
 Abs[Subscript[I, 2]]^2 + 
 2*Cos[2 Φ]*
  Re[Subscript[I, 0] Conjugate[Subscript[I, 2]]]),
 {Φ,0, π/2}, {v, 0, 10}]

and looks like this $w=141,20^2(\lvert{I_0}\rvert+4 \lvert{I_1}\rvert*cos^2(\Phi)+\lvert{I_2}^2\lvert+2 cos(2\Phi)*\Re{(I_0 I_2^*)})$

Before that I have defined I_0, I_1 and I_2. They are dependent from the variable v. and

α = 45/180 π

Subscript[I, 0] = 
NIntegrate[
Sqrt[Cos[θ]] Sin[θ] (1 + Cos[θ]) BesselJ[
0, (v Sin[θ])/Sin[α]], {θ, 0, α}, 
PrecisionGoal -> 2, MaxRecursion -> 20]


Subscript[I, 1] = 
NIntegrate[
Sqrt[Cos[θ]] Sin[θ]^2 BesselJ[
1, (v Sin[θ])/Sin[α]] , {θ, 0, α}, 
PrecisionGoal -> 2, MaxRecursion -> 20]

Subscript[I, 2] = NIntegrate[
Sqrt[Cos[θ]] Sin[θ] (1 - Cos[θ]) BesselJ[
2, (v Sin[θ])/Sin[α]] , {θ, 0, α}, 
PrecisionGoal -> 2, MaxRecursion -> 20]

What is wrong with my code?

Can I use Contours -> {} and ContourLabels -> All with ParametricPlot aswell, or is it just for ContourPlot?

There is a paper focusing on this problem:http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/29/032/29032386.pdf The Intensities I_0,1,2are on page 211 of the paper and is page 35 of the pdf, the plots are on page 218 of this paper and are pdf-page 42.

$\endgroup$
  • $\begingroup$ I think you want ContourPlot rather than ParametricPlot. What are the expressions for I_0, I_1, and I_2? What is the y-axis (is it Phi?). Also, don't use subscripts. $\endgroup$ – march Jul 9 '15 at 17:49
  • $\begingroup$ I used ParametricPlot because I want something like a unit circle. There isnt really a y-axis, just the angles. $\endgroup$ – donut Jul 9 '15 at 17:53
  • $\begingroup$ So is w the variable corresponding to the y-axis then? In other words, is $\tan(\Phi) = w/v$? Or is $w$ the "$z$" value, i.e. the values of the contours shown? Also, in order to be able to reproduce what you are doing, we need the expressions for the three I's. Also note that I is a protected symbol (the complex unit) in Mathematica. It is a good idea to use lower-case symbols for variable names in Mathematica since built-in functions are all capitalized. $\endgroup$ – march Jul 9 '15 at 17:55
  • $\begingroup$ Better yet: can you link to the book/article from where you obtained this figure? $\endgroup$ – J. M.'s technical difficulties Jul 9 '15 at 17:59
  • $\begingroup$ no w is not corresponding to y. Thats why the first circle in the middle has the value w=90 and the ''value of the y'' and v are almost 0 $\endgroup$ – donut Jul 9 '15 at 17:59
6
$\begingroup$

It seems to me that what you want is a ContourPlot, with a strong caveat (see below). Your function is defined in polar coordinates, so in order to use ContourPlot, we need to transform to polar coordinates, as:

v -> Sqrt[x^2 + y^2]
phi -> ArcTan[x, y]

Thus, replace any instance of v or phi in w with the corresponding expression in terms of of x and y, and call ContourPlot as

ContourPlot[ w
  , {x, -10, 10}, {y, -10, 10}
  , PlotRange -> All
  , ContourShading -> False
  , Contours -> {0.25, 0.5, 1.0, 2.0, 3.0, 10.0, 30.0, 50.0, 70.0, 90.0}
]

where

w = 141.20^2/(16*Pi)*(
  Abs[i0[Sqrt[x^2 + y^2]]]^2 +
  Abs[i2[Sqrt[x^2 + y^2]]]^2 + 
  4 Abs[i1[Sqrt[x^2 + y^2]]]^2 Cos[ArcTan[y, x]] + 
  2 Cos[2 ArcTan[y, x]] Re[i0[Sqrt[x^2 + y^2]] Conjugate[i2[Sqrt[x^2 + y^2]]]]
 )

The next issue is that your expressions for I_0, I_1, and I_2 are functions defined as numerical integrals, but the argument depend on v. So, we define these functions with SetDelay (:=) as

a = 45*Pi/180.;
Clear[i0, i1, i2]

i0[v_] :=  NIntegrate[
  Sqrt[Cos[t]] Sin[t] (1 + Cos[t]) BesselJ[0, (v Sin[t])/Sin[a]]
  , {t, 0, a}
  , PrecisionGoal -> 2, MaxRecursion -> 20]

i1[v_] :=  NIntegrate[
  Sqrt[Cos[t]] Sin[t]^2 BesselJ[1, (v Sin[t])/Sin[a]]
  , {t, 0, a}
  , PrecisionGoal -> 2, MaxRecursion -> 20]

i2[v_] :=  NIntegrate[
  Sqrt[Cos[t]] Sin[t] (1 - Cos[t]) BesselJ[2, (v Sin[t])/Sin[a]]
  , {t, 0, a}
  , PrecisionGoal -> 2, MaxRecursion -> 20]

Now, the problem here is that since these functions are defined in terms of integrals, the numerics of this problem cause it to be extremely slow. To see how long it will take, just try to find one contour with a restricted range via

ContourPlot[
  w
 , {y, -2, 2}, {x, -2, 2}
 , PlotRange -> All
 , ContourShading -> False
 , Contours -> {50.0}
]

This results in the image shown below, but it takes (on my machine) 45 seconds. Alternatively, you can evaluate the function on a fine-enough grid and use ListContourPlot. I did

points = Flatten[
   Table[{x, y, w}, {x, 0.000001, 10.000001, 0.1}, {y, 0.000001, 10.000001, 0.1}]
   , 1]

(Flatten is necessary in order to turn make sure that points is a List of ordered triples.) This takes a good amount of time (3 minutes on my machine), but you only need to do it once. Then you can make the contour plot and play around with the formatting all you want without having to re-evaluate the function. (Note that we need to avoid x = 0 in our parameterization of the function; hence the 0.000001's. You can extend this to negative values of x and y in obvious ways.) Then, we call

ListContourPlot[points
 , PlotRange -> All
 , ContourShading -> False
 , Contours -> {0.25, 0.5, 1.0, 2.0, 3.0, 10.0, 30.0, 50.0, 70.0, 90.0}]

resulting in the second figure below.

Another possible solution is to numerically integrate the I functions on a grid first, then create InterpolatingFunctions that can then be used in ContourPlot. The next step is to see if that works.

enter image description here

Pedagogical Note

The reason ParametricPlot is not the right answer here is that the correct format for that type of plot is

ParametricPlot[ {x[t], y[t]}, {t, 0, 1} ]

or, if you are in polar coordinates

ParametricPlot[ v[t]*{Cos[phi[t]], Sin[phi[t]]}, {t, 0, 1} ]

So, your idea of putting w in as the argument for ParametricPlot and trying to use two independent variables v and phi doesn't work for this type of plot.

Here, you could parameterize, say, v in terms of phi by solving

0.25 = w

for v to get the 0.25 contour of w (although there are multiple contours, so you would need to be careful), but you don't have any way of doing this analytically, and if you are doing things numerically, you should just ContourPlot anyway.

| improve this answer | |
$\endgroup$
  • $\begingroup$ no this isnt it, because my u isnt v*tan[phi]. Its not Cartesian coordinates $\endgroup$ – donut Jul 9 '15 at 18:31
  • $\begingroup$ thank you so much, this helps me alot. now I understand how to apply it here $\endgroup$ – donut Jul 10 '15 at 9:40
  • $\begingroup$ sorry march I still have one question. I tied to reproduce it and everything worked fine till ListContourPlot. I'm getting a nothing there. just a plane koordinate. Then I tried it with ContourPlot. That solution was not quite right. The contours in the middle are not shown.Then I tried Contourplotwith only the inner Contours. It gave me just one circle. $\endgroup$ – donut Jul 10 '15 at 10:23
  • $\begingroup$ @donut. See the updated answer. $\endgroup$ – march Jul 10 '15 at 17:00
6
$\begingroup$

Continuing with the InterpolatingFunctions suggested my March.

I was getting errors unless I set a finite AccuracyGoal.

I created data in the range -2 to 20 in steps of 0.1 (the range of interest is 0 to 10 so wanted data beyond this range) and then used that as input to the Interpolation.

i0[v_] := 
  NIntegrate[
   Sqrt[Cos[t]] Sin[t] (1 + Cos[t]) BesselJ[0, (v Sin[t])/Sin[a]], {t,
     0, a}, PrecisionGoal -> 2, MaxRecursion -> 20, 
   AccuracyGoal -> 10^-8];

d0 = Map[{#, i0[#]} &, Range[-2, 20, 0.1]];

e0 = Interpolation[d0];

To test it I plotted the original function with the InterpolationFunction (here I used the Range -20 to 20).

Show[
 Plot[i0[v], {v, -10, 10}, PlotStyle -> {Thick, Blue}],
 Plot[e0[v], {v, -10, 10}, PlotStyle -> {Thick, White, Dashed}]
 ]

Mathematica graphics

I repeated these steps for i1 and i2 resulting in e1 and e2 InterpolationFunctions.

The contour plot using the numerical integration is painfully slow.

It took about 5 minutes to produce on my machine.

w = 141.20^2/(16*Pi)*(Abs[i0[Sqrt[x^2 + y^2]]]^2 + 
    Abs[i2[Sqrt[x^2 + y^2]]]^2 + 
    4 Abs[i1[Sqrt[x^2 + y^2]]]^2 Cos[ArcTan[y, x]] + 
    2 Cos[2 ArcTan[y, x]] Re[
      i0[Sqrt[x^2 + y^2]] Conjugate[i2[Sqrt[x^2 + y^2]]]])

and then

ContourPlot[w,
 {x, -10, 10},
 {y, -10, 10},
 ContourShading -> False,
 Contours -> {0.25, 0.5, 1.0, 2.0, 3.0, 10.0, 30.0, 50.0, 70.0, 90.0}
 ]

Mathematica graphics

I then tried essentially the same function replacing it with the interpolated functions.

we = 141.20^2/(16*Pi)*(Abs[e0[Sqrt[x^2 + y^2]]]^2 + 
     Abs[e2[Sqrt[x^2 + y^2]]]^2 + 
     4 Abs[e1[Sqrt[x^2 + y^2]]]^2 Cos[ArcTan[y, x]] + 
     2 Cos[2 ArcTan[y, x]] Re[
       e0[Sqrt[x^2 + y^2]] Conjugate[e2[Sqrt[x^2 + y^2]]]]);

ContourPlot[we,
 {x, -10, 10},
 {y, -10, 10},
 ContourShading -> False,
 Contours -> {0.25, 0.5, 1.0, 2.0, 3.0, 10.0, 30.0, 50.0, 70.0, 90.0}
 ]

I was surprised to see that when I used the InterpolationFunction I got some errors.

Mathematica graphics

The CountourPlot did a good job in the middle but had some countours not present in the original plot.

Mathematica graphics

This is where I have to stop. I judge the interpolation functions to be 80 percent successful.

| improve this answer | |
$\endgroup$
  • $\begingroup$ This is very helpful, thank you $\endgroup$ – donut Jul 10 '15 at 9:41

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.