# Flamm's Paraboloid for the Kerr Metric

Mar 1st, 2018
118
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. (* Yukterez, version 02.03.2018 || syntax: Wolfram Mathematica || kerr.yukterez.net || flamm.yukterez.net *)
2.
3. ClearAll["Global`*"]
4. Quiet[Do[
5.
6. mt1={"StiffnessSwitching",Method->{"ExplicitRungeKutta",Automatic}};
7. mt2={"EventLocator","Event"->(r[t]-1000001/1000000 rA)};
8. mt3={"ImplicitRungeKutta","DifferenceOrder"->20};
9. mt4={"EquationSimplification"->"Residual"};
10. mt0=Automatic;
11. mta=mt0;
12. wp=MachinePrecision;
13.
14. a; A=a;
15. crd=1; dsp=1;
16. tmax=If[a>99/100, 500, 120]; Tmax=tmax;
17. TMax=Min[Tmax,т[plunge-1*^-4]];tMax=Min[tmax,plunge];
18.
19. r0=Sqrt[2]10 mirr;
20. θ0=Pi/2;
21. φ0=0;
22. μ=-0;
23. v0=1;
24. α0=-Pi/2;
25. ψ0=0;
26.
27. vr0=v0 Sin[α0];
28. vθ0=v0 Cos[α0] Sin[ψ0];
29. vφ0=v0 Cos[α0] Cos[ψ0];
30. x0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];
31. y0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
32. z0[A_]:=r0 Cos[θ0];
33.
34. x0KS[A_]:=((r0 Cos[φ0]+A Sin[φ0]) Sin[θ0]);
35. y0KS[A_]:=((r0 Sin[φ0]-A Cos[φ0]) Sin[θ0]);
36.
37. x0[A_]:=If[crd==1,x0BL[A],x0KS[A]];
38. y0[A_]:=If[crd==1,y0BL[A],y0KS[A]];
39.
40. ε=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0;
41. Lz=vφ0 Ы/j[v0];
42. pθ0=vθ0 Sqrt[Ξ]/j[v0];
43. pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
44. Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2,θ1->θ0];
45. k=Q+Lz^2+a^2 (ε^2+μ);
46.
47. rE=1+Sqrt[1-a^2 Cos[θ]^2]; RE[A_,w1_,w2_]:=Xyz[xyZ[{Sqrt[rE^2+A^2] Sin[θ] Cos[φ],Sqrt[rE^2+A^2] Sin[θ] Sin[φ],rE Cos[θ]},w1],w2];
48. rG=1-Sqrt[1-a^2 Cos[θ]^2]; RG[A_,w1_,w2_]:=Xyz[xyZ[{Sqrt[rG^2+A^2] Sin[θ] Cos[φ],Sqrt[rG^2+A^2] Sin[θ] Sin[φ],rG Cos[θ]},w1],w2];
49. rA=1+Sqrt[1-a^2]; RA[A_,w1_,w2_]:=Xyz[xyZ[{Sqrt[rA^2+A^2] Sin[θ] Cos[φ],Sqrt[rA^2+A^2] Sin[θ] Sin[φ],rA Cos[θ]},w1],w2];
50. rI=1-Sqrt[1-a^2]; RI[A_,w1_,w2_]:=Xyz[xyZ[{Sqrt[rI^2+A^2] Sin[θ] Cos[φ],Sqrt[rI^2+A^2] Sin[θ] Sin[φ],rI Cos[θ]},w1],w2];
51.
52. j[v_]:=Sqrt[1+μ v^2];
53. mirr=Sqrt[(Sqrt[1-a^2]+1)/2];
54. я=Sqrt[Χ/Σ]Sin[θ[τ]];
55. яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
56. Ы=Sqrt[χ/Ξ]Sin[θ0];
57. Σ=r[τ]^2+a^2 Cos[θ[τ]]^2;
58. Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
59. Ξ=r0^2+a^2 Cos[θ0]^2;
60. Δ=r[τ]^2-2r[τ]+a^2;
61. Δi[τ_]:=R[τ]^2-2R[τ]+a^2;
62. δ=r0^2-2r0+a^2;
63. Χ=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
64. Χi[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
65. χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
66. ц=2r[τ]/Σ;ц0=2r0/Ξ;
67.
68. т[τ_]:=Evaluate[t[τ]/.sol][[1]];
69. д[ξ_]:=Quiet[Ξ/.FindRoot[т[Ξ]-ξ,{Ξ,0}]];
70. T:=Quiet[д[tk]];
71.
72. ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
73. γ[τ_]:=ю[τ];
74. R[τ_]:=Evaluate[r[τ]/.sol][[1]];
75. Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];
76. Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
77. ß[τ_]:=Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2]/ю[τ]];
78.
79. ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]];ς0=Sqrt[χ/δ/Ξ];
80. ω[τ_]:=2R[τ] a/Χi[τ];ω0=2r0 a/χ;ωd=2r[τ] a/Χ;
81. Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];
82. й[τ_]:=ω[τ] яi[τ] ς[τ];й0=ω0 Ы ς0;
83. ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ];ж0=Sqrt[ς0^2-1]/ς0;
84.
85. vd[τ_]:=Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2-r[τ]^4(ε-Lz ωd)^2+Δ(Σ+a^2 Sin[θ[τ]]^2 (ε-Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+a^2 Sin[θ[τ]]^2 Δ](ε-Lz ωd)))];
86.
87. v[τ_]:=If[μ==0,1,Evaluate[vlt'[τ]/.sol][[1]]];
88. dst[τ_]:=Evaluate[str[τ]/.sol][[1]];
89. pΘ[τ_]:=Evaluate[pθ[τ]/.sol][[1]];pΘks[τ_]:=Σi[τ] Θ'[τ];
90. pR[τ_]:=Evaluate[pr[τ]/.sol][[1]];pRks[τ_]:=Σi[τ]/Δi[τ] R'[τ];
91. sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
92. epot[τ_]:=ε+μ-ekin[τ];
93. ekin[τ_]:=If[μ==0,ς[τ],1/Sqrt[1-v[τ]^2]-1];
94. dp=Y^Y;n0[z_]:=Chop[N[z]];
95.
96. pr2τ[τ_]:=1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ;
97. pθ2τ[τ_]:=(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ));
98.
99. DGL={t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),t[0]==0,r'[τ]==(pr[τ] Δ)/Σ,r[0]==r0,θ'[τ]==0,θ[0]==θ0,φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),φ[0]==φ0,pr'[τ]==pr2τ[τ],pr[0]==pr0,pθ'[τ]==0,pθ[0]==pθ0,str'[τ]==vd[τ]/Max[1*^-16,Abs[Sqrt[1-vd[τ]^2]]],str[0]==0,vlt'[τ]==vd[τ],vlt[0]==0};
100.
101. dr=(pr0 δ)/Ξ;dθ=pθ0/Ξ;dφ=(2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ;dΦ=If[θ0==0,0,If[θ0==\[Pi],0,dφ]];
102. φk=φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}];
103.
104. sol=NDSolve[DGL,{t,r,θ,φ,str,vlt,pr,pθ},{τ,0,tmax},WorkingPrecision->wp,MaxSteps->Infinity,Method->mta,InterpolationOrder->All,StepMonitor:>(laststep=plunge;plunge=τ;
105. stepsize=plunge-laststep;),Method->{"EventLocator","Event":>(If[stepsize<1*^-4,0,1])}];
106.
107. x[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]+n Pi/180]/.sol][[1]];
108. y[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]+n Pi/180]/.sol][[1]];
109. z[τ_]:=0;
110.
111. XKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+a Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
112. YKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-a Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
113.
114. X[τ_]:=If[crd==1,XBL[τ],XKS[τ]];
115. Y[τ_]:=If[crd==1,YBL[τ],YKS[τ]];
116.
117. xBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
118. yBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
119.
120. xKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
121. yKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
122.
123. XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2];XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];
124. Xyz[{x_,y_,z_},α_]:={x Cos[α]-y Sin[α],x Sin[α]+y Cos[α],z};
125. xYz[{x_,y_,z_},β_]:={x Cos[β]+z Sin[β],y,z Cos[β]-x Sin[β]};
126. xyZ[{x_,y_,z_},ψ_]:={x,y Cos[ψ]-z Sin[ψ],y Sin[ψ]+z Cos[ψ]};
127.
128. plp=50;
129. w1l=0;
130. w2l=0;
131. w1r=0;
132. w2r=0;
133. Mrec=300;
134. mrec=10;
135. imgsize=482;
136.
137. s[text_]:=Style[text,FontSize->font];font=11;
138.
139. tk=TMax;
141.
142. Σj[r_]:=r^2+a^2 Cos[θ0]^2;
143. Δj[r_]:=r^2-2 r+a^2;
144. Χj[r_]:=(r^2+a^2)^2-a^2 Sin[θ0]^2 Δj[r];
145.
146. grr[r_]:=Σj[r]/Δj[r];
147.
148. Ȓj[я_]:=Sqrt[NIntegrate[Sqrt[Abs[grr[r]]],{r,0,Abs[я]},Method->setj,MaxRecursion->100]^2+a^2];
149. гj[я_]:=Re[x/.FindRoot[Ȓj[x]-я,{x,3}]];
150.
151. Ȓ[я_]:=NIntegrate[Sqrt[Abs[grr[r]]],{r,0,Abs[я]},Method->setj,MaxRecursion->100];
152. Ř[я_]:=NIntegrate[r/Sqrt[r^2+a^2] Sqrt[Abs[(grr[r]-1)]],{r,Abs[я],r0},Method->setj,MaxRecursion->100];
153. ř[я_]:=Ř[0]-Ř[я];
154. г[я_]:=x/.FindRoot[Ȓ[x]-я,{x,1}];
155. z[я_]:=ř[я]-ř[r0];
156.
157. rint=Evaluate[Interpolation[Table[{tint,R[д[tint]]},{tint,0,tk,1}]]];
158. pint=Evaluate[Interpolation[Table[{tint,Φ[д[tint]]},{tint,0,tk,1}]]];
159. zint=Evaluate[Interpolation[Table[{tint,z[rint[tint]]},{tint,0,tk,1}]]];
160. xint[τ_,n_]:=Sqrt[rint[τ]^2+a^2] Sin[θ0] Cos[pint[τ]+n Pi/180];
161. yint[τ_,n_]:=Sqrt[rint[τ]^2+a^2] Sin[θ0] Sin[pint[τ]+n Pi/180];
162.
163. plot3D[PR_,w2_]:=
164. Rasterize[Grid[{{
165. Show[
166. ParametricPlot3D[Table[xYz[{xint[tt,n]/mirr, yint[tt,n]/mirr,zint[tt]/mirr},w2],{n,10,360,10}],{tt,0,tk},
167. PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic",
169. ParametricPlot3D[xYz[{2 Sin[uuu], 2 Cos[uuu],zint[tk]/mirr-0.005},w2],{uuu,0,2Pi},
170. PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic"],
171. If[a == 1,ParametricPlot3D[xYz[{2 Sin[uuu], 2 Cos[uuu], zint[tk]/mirr - uuu/100}, w2],
172. {uuu, 0, 2000}, PlotStyle -> {Thickness[0.0025]}, PlotPoints -> Automatic, PlotTheme -> "Classic",
173. ColorFunction -> Function[{x, y, z}, Darker[Hue[2/3, 1, 1, z]]]], {}],
174. ParametricPlot3D[xYz[{r0/mirr Sin[uuu], r0 /mirr Cos[uuu],zint[0]/mirr},w2],{uuu,0,2Pi},
175. PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic"]
176. ]},{s[" a"->N[a]]}},Alignment->Left]];
177.
178. Export["A-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,0]];
179. Export["B-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/8]];
180. Export["C-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/4]];
181. Export["D-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/2]];
182.
183. ClearAll["Global`*"];
184. ClearAll["Local`*"],
185.
186. {a,1,0,-1/5}]];
187.
188. (* Bildsequenz wird ins Hauptverzeichnis exportiert *)
RAW Paste Data