In[]:=
f[r_]:=1-
2M
2
r
3
r
+
3
g
+
8
3
πP
2
r
;
In[]:=
Solve[f[r]==0,M]
Out[]=
M
(3+8Pπ
2
r
)(
3
g
+
3
r
)
6
2
r

In[]:=
mass=
(3+8Pπ
2
r
)(
3
g
+
3
r
)
6
2
r
;
In[]:=
D[f[r],r]
4π
/.Mmass
Out[]=
16Pπr
3
-
23+8Pπ
2
r

3r
+
2
r
3+8Pπ
2
r

3
g
+
3
r
4π
In[]:=
FullSimplify[%]
Out[]=
-2
3
g
+
3
r
+8Pπ
5
r
4πr(
3
g
+
3
r
)
In[]:=
temp=
-2
3
g
+
3
r
+8Pπ
5
r
4πr(
3
g
+
3
r
)
;
In[]:=
ξ=
3
(3+8π
2
r
)
3
(
3
g
+
3
r
)
216
6
r
+
6
g
-
3
g
3
(3+8π
2
r
)
3
(
3
g
+
3
r
)
108
6
r
-
3
g
;
In[]:=
rp=
3
g
2
2
r
+
r
2
+
4
3
(
3
g
+
3
r
)Pπ+
2
(3+8π
2
r
)
2
(
3
g
+
3
r
)
36
4
r
1
3
ξ
+
1
3
ξ
;
In[]:=
rs=
rp
f[rp]
;
In[]:=
shadowcontour[met_,g_,p_,T1_,T2_]:=Module{T,r,θ,shadowrad,photonrad,xi,rh},​​rh=r/.LastNSolveT-
-2
3
g
+
3
r
+8pπ
5
r
4πr(
3
g
+
3
r
)
==0,r;​​xi=
3
(3+8π
2
r
)
3
(
3
g
+
3
r
)
216
6
r
+
3
g
-
3
g
3
(3+8π
2
r
)
3
(
3
g
+
3
r
)
108
6
r
-
3
g
;​​photonrad=
3
g
2
2
r
+
r
2
+
4
3
(
3
g
+
3
r
)pπ+
2
(3+8π
2
r
)
2
(
3
g
+
3
r
)
36
4
r
1
3
xi
+
1
3
xi
;​​shadowrad=
photonrad
met[photonrad]
/.rrh;​​plot=Region[ParametricRegion[{shadowradCos[θ],shadowradSin[θ]},{{θ,0,2π},{T,T1,T2}}]];​​plot​​​​
In[]:=
shadowcontourf,0.3,
0.002418
0.09
,0,0.21
Out[]=
Region
Embedding dimension: 2

In[]:=
g=0.3;​​P=
0.002418
0.09
;​​ParametricPlot[{{rsCos[θ],rsSin[θ]},{{θ,0,2π},{T,0,0.21}}}]
In[]:=
ParametricPlot[{rsCos[θ],rsSin[θ]},{θ,0,2π},{T,0,0.21}]
Out[]=