Kerr photon regions by He-Xu Zhang

--------------------务必使用12.1及以上版本打开

Δ[M_,a_]:=((((Sqrt[(x^2)+(z^2)])-M))^2)-2 M ((Sqrt[(x^2)+(z^2)])-M)+(a^2);
gtt[M_,a_]:=((((Sqrt[(x^2)+(z^2)])-M))^2)-2 M ((Sqrt[(x^2)+(z^2)])-M)+(a^2) ((z^2)/((x^2)+(z^2)));

PhReigon[M_,a_]:=(((4 ((Sqrt[(x^2)+(z^2)])-M) (((((Sqrt[(x^2)+(z^2)])-M))^2)-2 M ((Sqrt[(x^2)+(z^2)])-M)+(a^2))-(((((Sqrt[(x^2)+(z^2)])-M))^2)+(a^2) ((z^2)/((x^2)+(z^2)))) 2 (((Sqrt[(x^2)+(z^2)])-M)-M)))^2)-16 (a^2) ((((Sqrt[(x^2)+(z^2)])-M))^2) (((((Sqrt[(x^2)+(z^2)])-M))^2)-2 M ((Sqrt[(x^2)+(z^2)])-M)+(a^2)) ((x^2)/((x^2)+(z^2)));

(*RegionPlot[PhReigon[1,0.8]<0&&(x^2)+(z^2)>1,{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[0.788235,0.309803,0],PlotPoints->50,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]*)

(*Δr[M_,a_,r_]:=(r^2)-2 M r+(a^2);
NSolve[2 r Δr[1,0.8,r] D[Δr[1,0.8,r],r]+(r^2) (((D[Δr[1,0.8,r],r]))^2)-2 (r^2) Δr[1,0.8,r] D[Δr[1,0.8,r],{r,2}]==0,r]*)

aa=RegionPlot[PhReigon[1,0.8]<0&&(x^2)+(z^2)>(((1.288621))^2),{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[0.788235,0.309803,0],PlotPoints->200,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

bb=RegionPlot[PhReigon[1,0.8]<0&&(x^2)+(z^2)<(((1.288621))^2)&&(x^2)+(z^2)>1,{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[1,1,0.501960],PlotPoints->100,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

PhReigon1[M_,a_]:=(((4 M Log[((Sqrt[(x^2)+(z^2)])/M)] ((((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2)-2 M (M Log[((Sqrt[(x^2)+(z^2)])/M)])+(a^2))-((((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2)+(a^2) ((z^2)/((x^2)+(z^2)))) 2 ((M Log[((Sqrt[(x^2)+(z^2)])/M)])-M)))^2)-16 (a^2) (((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2) ((((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2)-2 M (M Log[((Sqrt[(x^2)+(z^2)])/M)])+(a^2)) ((x^2)/((x^2)+(z^2)));
(*cc=RegionPlot[PhReigon1[1,0.8]<0&&(x^2)+(z^2)>1&&(x^2)+(z^2)<(((1.288621))^2),{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[1,1,0.501960],PlotPoints->200,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]*)

dd=RegionPlot[PhReigon1[1,0.8]<0&&(x^2)+(z^2)<1,{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[0.788235,0.309803,0],PlotPoints->200,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

ee=RegionPlot[Δ[1,0.8]<0,{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->RGBColor[0,0.121568,0.501960],PlotPoints->100,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

ff=ContourPlot[gtt[1,0.8]==0,{x,-5,5},{z,-5,5},ContourStyle->RGBColor[0.501960,0.501960,0.501960],PlotPoints->50,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

gϕϕ[M_,a_]:=((((((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2)+(a^2)))^2)-(a^2) ((((M Log[((Sqrt[(x^2)+(z^2)])/M)]))^2)-2 M (M Log[((Sqrt[(x^2)+(z^2)])/M)])+(a^2)) ((x^2)/((x^2)+(z^2)));
gg=RegionPlot[gϕϕ[1,0.8]<0,{x,-5,5},{z,-5,5},BoundaryStyle->RGBColor[0.6,0.69019,1],PlotStyle->{RGBColor[0.6,0.69019,1],HatchFilling[((-π)/4)]},PlotPoints->50,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

(*gg0=ContourPlot[gϕϕ[1,0.8]==0,{x,-5,5},{z,-5,5},ContourStyle->RGBColor[0.6,0.69019,1],PlotPoints->50,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]*)

hh=ContourPlot[(x^2)+(z^2)==1,{x,-2,2},{z,-2,2},ContourStyle->{Dashed,Black},PlotPoints->50,Axes->True,AxesLabel->{x,z},AxesStyle->Thickness[0.004],Ticks->{{-4.8,4.8},None},TicksStyle->Black,Frame->None]

(*Show[aa,bb,dd,ee,ff,gg,hh]*)

(*reg1=ImplicitRegion[Δ[1,0.8]>0&&gtt[1,0.8]<0&&(x^2)+(z^2)>(2^2),{x,z}];
reg11=DiscretizeRegion[reg1,{{-5,5},{-5,5}}];
ii0=ContourPlot[2.5 x-z,{x,z}∈reg11,ContourShading->None,Contours->50,ContourStyle->RGBColor[0.501960,0.501960,0.501960]]*)

ii=RegionPlot[Δ[1,0.8]>0&&gtt[1,0.8]<0&&(x^2)+(z^2)>(2^2),{x,-5,5},{z,-5,5},BoundaryStyle->None,PlotStyle->{RGBColor[0.501960,0.501960,0.501960],HatchFilling[(π/4)]}]

(*reg2=ImplicitRegion[gϕϕ[1,0.8]<0,{x,z}];
reg22=DiscretizeRegion[reg2,{{-5,5},{-5,5}}];
jj=ContourPlot[1 x+z,{x,z}∈reg22,ContourShading->None,Contours->25,ContourStyle->RGBColor[0.6,0.69019,1]]*)

kk=Graphics[{{PointSize[Large],Point[{-1,0}]},{PointSize[Large],Point[{1,0}]}}](*画龙点睛，最后画奇点并把它放在Show的最后面也是防止它被其他图形盖住*)

Show[aa,bb,dd,ee,ff,gg,hh,ii,kk]