Yukterez

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;
  140. setj={"GlobalAdaptive","MaxErrorIncreases"->100,Method->"GaussKronrodRule"};
  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",
  168. PlotRange->{{-PR,PR},{-PR,PR},{-PR,PR}},ImageSize->imgsize,ImagePadding->1,ViewPoint->{Infinity,0,0}],
  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