Yukterez

lichtablenkung krümmungsradius

Sep 20th, 2018
55
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* |||| Mathematica Syntax | yukterez.net | tinyurl.com/ycxlabz3 |||||||||||||||||||||||| *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. ClearAll["Global`*"]
  6.  
  7. mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  8. mt2 = {"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
  9. mt3 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  10. mt4 = {"EquationSimplification"-> "Residual"};
  11. mt0 = Automatic;
  12. mta = mt3;
  13. wp = 32;
  14.  
  15. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  16. (* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
  17. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  18.  
  19. A = a; (* pseudosphärisch: A=0, kartesisch: A=a *)
  20. crd = 1; (* Boyer-Lindquist: crd=1, Kerr-Schild: crd=2 *)
  21. dsp = 1; (* Display Modus *)
  22.  
  23. tmax = 3; (* Eigenzeit *)
  24. Tmax = 3; (* Koordinatenzeit *)
  25. TMax = Min[Tmax, т[plunge-1*^-3]]; tMax = Min[tmax, plunge]; (* Integrationsende *)
  26.  
  27. r0 = 10^6; (* Startradius *)
  28. θ0 = π/2; (* Breitengrad *)
  29. φ0 = 0; (* Längengrad *)
  30. a = 78/10; (* Spinparameter *)
  31.  
  32. v0 = 1; (* Anfangsgeschwindigkeit *)
  33. α0 = 0; (* vertikaler Abschusswinkel *)
  34. ψ0 = 0; (* Bahninklinationswinkel *)
  35.  
  36. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  37. (* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
  38. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  39.  
  40. vr0 = v0 Sin[α0]; (* radiale Geschwindigkeitskomponente *)
  41. vθ0 = v0 Cos[α0] Sin[ψ0]; (* longitudinale Geschwindigkeitskomponente *)
  42. vφ0 = v0 Cos[α0] Cos[ψ0]; (* latitudinale Geschwindigkeitskomponente *)
  43.  
  44. x0BL[A_] := Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0]; (* Anfangskoordinaten *)
  45. y0BL[A_] := Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
  46. z0[A_] := r0 Cos[θ0];
  47.  
  48. x0KS[A_] := ((r0 Cos[φ0]+A Sin[φ0]) Sin[θ0]);
  49. y0KS[A_] := ((r0 Sin[φ0]-A Cos[φ0]) Sin[θ0]);
  50.  
  51. x0[A_] := If[crd==1, x0BL[A], x0KS[A]];
  52. y0[A_] := If[crd==1, y0BL[A], y0KS[A]];
  53.  
  54. ε = Sqrt[δ Ξ/χ]/j[v0]+Lz ω0; (* Energie und Drehimpulskomponenten *)
  55. Lz = vφ0 Ы/j[v0];
  56. pθ0 = vθ0 Sqrt[Ξ]/j[v0];
  57. pr0 = vr0 Sqrt[(Ξ/δ)/j[v0]^2];
  58. Q = Simplify[Limit[pθ0^2+(Lz^2 Csc[ϑ]^2-a^2 (ε^2+μ)) Cos[ϑ]^2, ϑ->θ0]]; (* Carter Q *)
  59. k = Q+Lz^2+a^2 (ε^2+μ); (* Carter k *)
  60. μ = If[Abs[v0]==1, 0, -1]; (* Baryon: μ=-1, Photon: μ=0 *)
  61.  
  62. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  63. (* |||||||| 3) RADIUS NACH GESCHWINDIGKEIT UND VICE VERSA ||||||||||||||||||||||||||||||| *)
  64. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  65.  
  66. rPro = 2 (1+Cos[2/3 ArcCos[-a]]); (* prograder Photonenorbitradius *)
  67. rRet = 2 (1+Cos[2/3 ArcCos[+a]]); (* retrograder Photonenorbitradius *)
  68. rTeo = 1+2 Sqrt[1-a^3/3] Cos[ArcCos[(1-a^2)/(1-a^2/3)^(3/2)]/3];
  69.  
  70. δp[r_,a_] := Quiet[δi/.NSolve[(a^4(-1+r)+2(-3+r)r^4+a^2r(6+r(-5+3 r))+ (* Eq. Ink. Winkel *)
  71. 4a Sqrt[a^2+(-2+r)r](a^2+3 r^2)Cos[δi]-a^2(3+r)(a^2+(-2+r)r)Cos[2δi])/(2r(a^2+
  72. (-2+r)r)(r^3+a^2(2+r)))==0&&δi<=π&&δi>=0,δi][[1]]];
  73.  
  74. vPro = (a^2-2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a+r0^(3/2)));(* Kreisgeschwindigkeit + *)
  75. vRet = (a^2+2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a-r0^(3/2)));(* Kreisgeschwindigkeit - *)
  76.  
  77. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  78. (* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
  79. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  80.  
  81. rE = 1+Sqrt[1-a^2 Cos[θ]^2]; (* äußere Ergosphäre *)
  82. rG = 1-Sqrt[1-a^2 Cos[θ]^2]; (* innere Ergosphäre *)
  83. rA = 1+Sqrt[1-a^2]; (* äußerer Horizont *)
  84. rI = 1-Sqrt[1-a^2]; (* innerer Horizont *)
  85.  
  86. RE[A_, w1_, w2_] := Xyz[xyZ[
  87. {Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2];
  88. RG[A_, w1_, w2_] := Xyz[xyZ[
  89. {Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
  90. RA[A_, w1_, w2_] := Xyz[xyZ[
  91. {Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2];
  92. RI[A_, w1_, w2_] := Xyz[xyZ[
  93. {Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2];
  94.  
  95. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  96. (* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *)
  97. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  98.  
  99. horizons[A_, mesh_, w1_, w2_] := Show[
  100. ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  101. Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
  102. ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  103. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
  104. ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  105. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
  106. ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  107. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
  108. BLKS := Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];
  109.  
  110. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  111. (* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  112. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  113.  
  114. j[v_] := Sqrt[1+μ v^2]; (* Lorentzfaktor *)
  115. mirr = Sqrt[(Sqrt[1-a^2]+1)/2]; (* irreduzible Masse *)
  116. я = Sqrt[Χ/Σ]Sin[θ[τ]]; (* axialer Umfangsradius *)
  117. яi[τ_] := Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
  118. Ы = Sqrt[χ/Ξ]Sin[θ0];
  119. ц = 2r[τ]/Σ; ц0=2r0/Ξ;
  120.  
  121. Σ = r[τ]^2+a^2 Cos[θ[τ]]^2; (* poloidialer Umfangsradius *)
  122. Σi[τ_] := R[τ]^2+a^2 Cos[Θ[τ]]^2;
  123. Ξ = r0^2+a^2 Cos[θ0]^2;
  124.  
  125. Δ = r[τ]^2-2r[τ]+a^2;
  126. Δi[τ_] := R[τ]^2-2R[τ]+a^2;
  127. δ = r0^2-2r0+a^2;
  128.  
  129. Χ = (r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
  130. Χi[τ_] := (R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
  131. χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
  132.  
  133. т[τ_] := Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit nach Eigenzeit *)
  134. д[ξ_] := Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 0}]]; (* Eigenzeit nach Koordinatenzeit *)
  135. T := Quiet[д[tk]];
  136.  
  137. ю[τ_] := Evaluate[t'[τ]/.sol][[1]];
  138. γ[τ_] := If[μ==0, "Infinity", ю[τ]]; (* totale ZD *)
  139. R[τ_] := Evaluate[r[τ]/.sol][[1]]; (* Boyer-Lindquist Radius *)
  140. Φ[τ_] := Evaluate[φ[τ]/.sol][[1]];
  141. Θ[τ_] := Evaluate[θ[τ]/.sol][[1]];
  142. ß[τ_] := Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ]];
  143.  
  144. gs[τ_] := (2 (R[τ]^2-a^2 Cos[Θ[τ]]^2) Sqrt[((a^2+R[τ]^2)^2- (* Gravitation *)
  145. a^2 (a^2+(R[τ]-2) R[τ]) Sin[Θ[τ]]^2)/(a^2+2 R[τ]^2+a^2 Cos[2 Θ[τ]])^2])/(R[τ]^2+a^2 Cos[Θ[τ]]^2)^2;
  146.  
  147. ς[τ_] := Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0 = Sqrt[χ/δ/Ξ]; (* gravitative ZD *)
  148. ω[τ_] := 2R[τ] a/Χi[τ]; ω0 = 2r0 a/χ; ωd=2r[τ] a/Χ; (* Frame Dragging *)
  149. Ω[τ_] := ω[τ] Sqrt[X[τ]^2+Y[τ]^2]; (* Frame Dragging beobachtete Geschwindigkeit *)
  150. й[τ_] := ω[τ] яi[τ] ς[τ]; й0 = ω0 Ы ς0; (* Frame Dragging lokale Geschwindigkeit *)
  151.  
  152. ж[τ_] := Sqrt[ς[τ]^2-1]/ς[τ]; ж0 = Sqrt[ς0^2-1]/ς0; (* Fluchtgeschwindigkeit *)
  153.  
  154. vd[τ_] := Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2-
  155. r[τ]^4(ε-Lz ωd)^2+Δ(Σ+a^2 Sin[θ[τ]]^2 (ε-
  156. Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+
  157. a^2 Sin[θ[τ]]^2 Δ](ε - Lz ωd)))];
  158.  
  159. v[τ_] := If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]]; (* lokale Dreiergeschwindigkeit *)
  160. dst[τ_] := Evaluate[str[τ]/.sol][[1]]; (* Strecke *)
  161.  
  162. pΘ[τ_] := Evaluate[pθ[τ] /. sol][[1]]; pΘks[τ_] := Σi[τ] Θ'[τ]; (* Impuls *)
  163. pR[τ_] := Evaluate[pr[τ] /. sol][[1]]; pRks[τ_] := Σi[τ]/Δi[τ] R'[τ];
  164. sh[τ_] := Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
  165. epot[τ_] := ε+μ-ekin[τ]; (* potentielle Energie *)
  166. ekin[τ_] := If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1]; (* kinetische Energie *)
  167.  
  168. (* beobachtete Inklination *)
  169. ink0 := б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]];
  170.  
  171. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  172. (* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  173. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  174.  
  175. dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_] := Chop[N[Simplify[z]]];
  176.  
  177. (* Boyer-Lindquist-Koordinaten *)
  178.  
  179. pr2τ[τ_] := 1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ;
  180. pθ2τ[τ_] := (Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ));
  181.  
  182. DG1={
  183. t'[τ] == ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),
  184. t[0] == 0,
  185.  
  186. r'[τ] == (pr[τ] Δ)/Σ,
  187. r[0] == r0,
  188.  
  189. θ'[τ] == pθ[τ]/Σ,
  190. θ[0] == θ0,
  191.  
  192. φ'[τ] == (2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),
  193. φ[0] == φ0,
  194.  
  195. pr'[τ] == pr2τ[τ],
  196. pr[0] == pr0,
  197.  
  198. pθ'[τ] == pθ2τ[τ],
  199. pθ[0] == pθ0,
  200.  
  201. str'[τ]== vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
  202. str[0] == 0,
  203.  
  204. vlt'[τ]== vd[τ],
  205. vlt[0] == 0
  206. };
  207.  
  208. (* Kerr-Schild-Koordinaten *)
  209.  
  210. dr = (pr0 δ)/Ξ; dθ=pθ0/Ξ;
  211. dφ = (2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ;
  212. dΦ = If[θ0==0, 0, If[θ0==π, 0, dφ]];
  213. φk = φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}];
  214.  
  215. DG2={
  216. t''[τ] == (2 ((a^4 Cos[θ[τ]]^4+a^2 Cos[θ[τ]]^2 r[τ]-r[τ]^3-r[τ]^4) r'[τ]^2+r[τ] ((a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-2 a^3 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^3 θ'[τ] φ'[τ]+Sin[θ[τ]]^2 (r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]^2+a t'[τ] (a (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+2 (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]))+r'[τ] ((a^4 Cos[θ[τ]]^4+2 a^2 Cos[θ[τ]]^2 r[τ]-2 r[τ]^3-r[τ]^4) t'[τ]+a (a r[τ] (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+(-a^4 Cos[θ[τ]]^4-2 a^2 Cos[θ[τ]]^2 r[τ]+2 r[τ]^3+r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^3,
  217. t'[0] == Limit[(ц0 (-dr+a Sin[θ1]^2 dΦ))/(-1+ц0)+\[Sqrt]((1/((-1+ц0)^2 Ξ))(Ξ (dr^2+(-1+ц0) (μ-Ξ dθ^2))+Sin[θ1]^2 dΦ (-2a Ξ dr-(-1+ц0) χ dΦ+ц0^2 a^2 Ξ Sin[θ1]^2 dΦ))), θ1->θ0],
  218. t[0] == 0,
  219.  
  220. r''[τ] == (-8 (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2 Cos[2 θ[τ]]+r[τ] (2+r[τ])) r'[τ]^2+4 r'[τ] (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) (-2 r[τ]+a^2 Sin[θ[τ]]^2) t'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]-a Sin[θ[τ]]^2 (2 r[τ] (a^2 Cos[θ[τ]]^2 (-4+a^2+a^2 Cos[2 θ[τ]])+2 r[τ] ((2+a^2+a^2 Cos[2 θ[τ]]) r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+a^4 Sin[2 θ[τ]]^2) φ'[τ])+2 (a^2+(-2+r[τ]) r[τ]) (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+4 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+8 a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2))/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
  221. r'[0] == dr,
  222. r[0] == r0,
  223.  
  224. θ''[τ] == (4 a^2 r[τ] Sin[2 θ[τ]] (r'[τ]+t'[τ])^2-8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]^2-8 a ((Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]]) r'[τ]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ]) φ'[τ]+(2 a^6 Cos[θ[τ]]^4+r[τ] (a^4 Cos[θ[τ]]^2 (5+Cos[2 θ[τ]]) r[τ]+2 a^2 (2+Cos[2 θ[τ]]) r[τ]^3+2 r[τ]^5+2 a^2 (a^2 (3+Cos[2 θ[τ]])+4 r[τ]^2) Sin[θ[τ]]^2)) Sin[2 θ[τ]] φ'[τ]^2)/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
  225. θ'[0] == dθ,
  226. θ[0] == θ0,
  227.  
  228. φ''[τ] == If[θ[τ]==0, 0, (4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) r'[τ]^2+4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]^2+4 a r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-8 (a^2 Cos[θ[τ]]^2+r[τ]^2) (Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 r[τ] Sin[2 θ[τ]]) θ'[τ] φ'[τ]+a Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2+8 a t'[τ] (2 Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ]+a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ])+8 r'[τ] ((a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]+a Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ] (2+r[τ])) θ'[τ]-(r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]))/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3)],
  229. φ'[0] == dΦ,
  230. φ[0] == φk,
  231.  
  232. str'[τ]== vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
  233. str[0] == 0,
  234.  
  235. vlt'[τ]== vd[τ],
  236. vlt[0] == 0
  237. };
  238.  
  239. DGL = If[crd==1, DG1, DG2];
  240.  
  241. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  242. (* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  243. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  244.  
  245. sol=
  246.  
  247. NDSolve[DGL, {t, r, θ, φ, str, vlt, pr, pθ}, {τ, 0, tmax},
  248. WorkingPrecision-> wp,
  249. MaxSteps-> Infinity,
  250. Method-> mta,
  251. InterpolationOrder-> All,
  252. StepMonitor :> (laststep=plunge; plunge=τ;
  253. stepsize=plunge-laststep;), Method->{"EventLocator",
  254. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  255.  
  256. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  257. (* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  258. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  259.  
  260. XBL[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
  261. YBL[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  262. Z[τ_] := Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
  263.  
  264. XKS[τ_] := Evaluate[((r[τ] Cos[φ[τ]]+a Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  265. YKS[τ_] := Evaluate[((r[τ] Sin[φ[τ]]-a Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  266.  
  267. X[τ_] := If[crd==1, XBL[τ], XKS[τ]];
  268. Y[τ_] := If[crd==1, YBL[τ], YKS[τ]];
  269.  
  270. xBL[τ_] := Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
  271. yBL[τ_] := Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  272. z[τ_] := Z[τ];
  273.  
  274. xKS[τ_] := Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  275. yKS[τ_] := Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  276.  
  277. x[τ_] := If[crd==1, xBL[τ], xKS[τ]];
  278. y[τ_] := If[crd==1, yBL[τ], yKS[τ]];
  279.  
  280. XYZ[τ_] := Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_] := Sqrt[X[τ]^2+Y[τ]^2]; (* kartesisches R *)
  281.  
  282. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *)
  283. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  284. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  285.  
  286. α[τ_] := ArcTan[y'[τ], x'[τ]]; (* Winkel *)
  287. ř[τ_] := Sqrt[y'[τ]^2+x'[τ]^2]/α'[τ]; (* Krümmungsradius *)
  288.  
  289. "r"->r0
  290. "ř"->ř[0]
RAW Paste Data