1
$\begingroup$

I'm trying to get the following intensity pattern from this paper

enter image description here

In the Newtonian case, the equations for the intensity (eq. 9 and 10 in the article), I put them in Mathematica with the change $b^2 = x^2 + y^2$ because the intensity is circularly symmetric with the parameter impact $b$ as radius (page 8 in this article)

H[x_, y_] := 10^3 1/(32 \[Pi] (Sqrt[x^2 + y^2])^3);
G[x_, y_] := 10^3 (1/(32 \[Pi]^2 (Sqrt[x^2 + y^2])^3)) ArcTan[Sqrt[x^2 + y^2]/Sqrt[4 - (Sqrt[x^2 + y^2])^2]] -10^3 Sqrt[4 - (Sqrt[x^2 + y^2])^2]/(128 \[Pi]^2 (Sqrt[x^2 + y^2])^2);
combinedFunction[x_, y_] := Piecewise[{{G[x, y], x^2 + y^2 < 2}, {H[x, y], x^2 + y^2 > 2}}] 
DensityPlot[combinedFunction[x, y], {x, -10, 10}, {y, -10, 10}, ColorFunction -> "SunsetColors", Exclusions -> All,PlotPoints -> 100, PlotLegends -> Automatic]

and I get

enter image description here

However, it does not look like the image in the article. Could anyone advise me or give me a suggestion? You should obtain a circumference of radius $2$ as in the paper. I'm also trying to get the intensity pattern for a more general case where numerical integration was performed (I actually asked the question here) but I can't find how to make the change $b^2 = x^2 + y^2$. I appreciate suggestions or observations

$\endgroup$
3
  • 1
    $\begingroup$ Your plot is only showing a range up to ~0.2, and showing everything above that as white. Add PlotRange -> {0, 1.2}. $\endgroup$
    – MelaGo
    Nov 3, 2023 at 22:02
  • 1
    $\begingroup$ Also, shouldn't the combinedFunction be Piecewise[{{G[x, y], Sqrt[x^2 + y^2] < 2}, {H[x, y], Sqrt[x^2 + y^2] > 2}}]? $\endgroup$
    – MelaGo
    Nov 3, 2023 at 22:03
  • $\begingroup$ @MelaGo You're right. I honestly didn't notice that detail. Thank you so much! $\endgroup$ Nov 3, 2023 at 22:26

0