# Plot of massless particles on unstable null circular geodesic : Finding an alternative for GraphicsMeshFindIntersections and reducing code length

My friend and I were trying to reproduce Figure 4(left) from this paper(page 7). We succeeded in reproducing the result after writing a lengthy and complicated code. Here I dissect the code in steps:

Step 1:
Find the parametric plot for Free Energy v/s Temperature and the Transition Temperature.

l = 1;
T[rp_, Q_] := 1/(4 \[Pi] rp) (1 - Q^2/rp^2 + (3 rp^2)/l^2)
M[rp_, Q_] := rp/2 (1 + Q^2/rp^2 + rp^2/l^2)
S[rp_] := \[Pi] rp^2
\[Phi][rp_, Q_] := Q/rp
F[rp_, Q_] := Simplify[M[rp, Q] - (T[rp, Q] S[rp])]
Q = Qt l;
rp = rpt l;
Tt[rpt_, Qt_] = T[rp, Q] l;
Ft[rpt_, Qt_] = F[rp, Q]/l;
Mt[rpt_, Qt_] = M[rp, Q]/l;
r = rt l;
pplt[Qt_] :=
ParametricPlot[Evaluate@{Tt[rpt, Qt], Ft[rpt, Qt]}, {rpt, 0.01, 3},
PlotRange -> {{0.2, 0.6}, {-0.05, 0.15}}, AspectRatio -> 1]
ip[Qt_] := GraphicsMeshFindIntersections[pplt[Qt]][[1]]
Tp = Table[ip[Qt][[1]], {Qt, 0.01, 0.166, 0.001}];


Step 2:
Find the critical values and generate a Table for Qtt(Charge).

rptc = SolveValues[D[Tt[rpt, Qt], {rpt, 2}] == 0, rpt][[2]];
Qtc = SolveValues[D[Tt[rpt, Qt], rpt] == 0, Qt][[2]];
rpt1 = rptc /. Qt -> Qtc;
rptc = Last @@ SolveValues[rpt1 == rpt, rpt];
Qtc = Qtc /. rpt -> rptc;
Ttc = Tt[rpt, Qt] /. rpt -> rptc /. Qt -> Qtc;
\[Phi]c = \[Phi][rp, Q] /. Qt -> Qtc /. rpt -> rptc;
f[r_, rpt_, Qt_] := 1 - (2 Mt[rpt, Qt])/r + Q^2/r^2 + r^2/l^2
Veff[r_] := f[r, rpt, Qt] (L^2/r^2 + \[Delta]1)
r0 = l/2 (3/2 rpt (Qt^2/rpt^2 + rpt^2 + 1) +
Sqrt[9/4 rpt^2 (Qt^2/rpt^2 + rpt^2 + 1)^2 - 8 Qt^2]);
l = 1; \[Delta]1 = 0;
\[Lambda][rpt_, Qt_] =
Sqrt[-((r0^2 f[r0, rpt, Qt])/L^2) Veff''[r0]] // Simplify;
Qtt = Table[Qt, {Qt, 0.01, 0.166, 0.001}];


Step 3:
Find the parametric plot for Lyapunov exponent v/s Temperature and find the intersection points.

pplt1[Qt_] :=
ParametricPlot[{Tt[rpt, Qt], \[Lambda][rpt, Qt]}, {rpt, 0.01, 3},
PlotRange -> {{0.0, 0.6}, {0.75, 26.00}}, AspectRatio -> 1,
AxesLabel -> {"T", "\[Lambda]"}, PlotPoints -> 500]
pts = Table[
GraphicsMeshFindIntersections[
Show[pplt1[Qtt[[t]]],
Graphics[{Red, Line[{{Tp[[t]], 1}, {Tp[[t]], 26}}]}]]], {t, 1,
Length[Tp]}];


which gives the following result after using animate:
Step 4:
Find the minimum and maximum intersection points and find their difference to get the raw plot:

minmax = ReplacePart[pts, {{_, _, 1} -> Nothing, {_, 2 | 2} -> Nothing}];
\[CapitalDelta]\[Lambda] = Flatten[Differences[#] & /@ minmax];
data = Transpose@{Tp, \[CapitalDelta]\[Lambda]};
ListLinePlot[data, PlotRange -> All]


Step 5
Find the rescaled(normalized) plot using the critical values.

\[Lambda]c = \[Lambda][rptc, Qtc] // N;
diff\[Lambda] = \[CapitalDelta]\[Lambda]/\[Lambda]c;
t1 = Tp/Ttc;
data2 = Transpose@{t1, diff\[Lambda]};
ListLinePlot[data2, PlotRange -> All]


I wish to use functional programming techniques to reduce the code length to get the desired result. Is there a way to shorten the code? Any help in this regard would be truly beneficial!

• Why the downvote?
Sep 10, 2022 at 22:49
• Other than the downvote, how is this a trivial syntax error? I do not understand
– bmf
Sep 11, 2022 at 13:25
• I'm not the downvoter, but I'd say the question is too localized and unlikely to help future visitors as it currently written. Your code is already Mathematica-style, I don't see any bad coding practice. I'm afraid it's hard to give advice without understanding the purpose of the code, but the code is lenghy and it's not immediately clear to me what problem is solved here. Can you summarize the problem with a few lines (probably with some traditional math formulas)? That'll make the question on-topic. Sep 12, 2022 at 5:43
• @xzczd basically I am trying to find a parametric for $\frac{\Delta \lambda }{\lambda_c}$ v/s $t$ where $\lambda_c$ is the critical Lyapunov exponent and t is the rescaled Temperature given by $\frac{T_p}{T_c}$. I used the command "GraphicsMeshFindIntersections" to generate a table from the previous plots, which makes the code lengthy as I use the data further in my calculations. I wished to find a novel approach to get the desired result, without using some intermediate steps. Hope it helps. Sep 12, 2022 at 7:20
• @GustavoDelfino thank you for the feedback. Actually, I have to repeat the same process for massive particles as well, where I have a different equation for $\lambda$. I wish there was an alternative for 'GraphicsMeshFindIntersections command, i.e. I can get the desired result without using the intermediate plots. Sep 12, 2022 at 15:05

This is a more functional way to do step 1:

In[1]:= step1 = RightComposition[
Range,
Map[Qt|->{(1-Qt^2/rpt^2+3 rpt^2)/(4 \[Pi] rpt),(3 Qt^2+rpt^2-rpt^4)/(4 rpt)}],
Map[ParametricPlot[#,{rpt,0.01,3},PlotRange->{{0.2,0.6},{-0.05,0.15}},AspectRatio->1]&],
Map[GraphicsMeshFindIntersections],
Map[First],
Map[First]
];

In[2]:= step1[0.01,0.166,0.001] == Tp
Out[2]= True
`

I have no time now to help with the other steps. But this is certainly more in the style of functional programming.

• Thank you for the informative answer. It gave me a head start. I will try to build upon your suggestion! Sep 12, 2022 at 15:37