Yukterez

Schwarzschild deflection of light

Oct 12th, 2017
178
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* schwarzschild.yukerez.net *) (* gravitylense.yukterez.net *) (* Simon Tyran, Vienna *)
  2. (* |||||||||||||||||||||||||||||||| Syntax: Mathematica |||||||||||||||||||||||||||||||| *)
  3.  
  4. ClearAll["Global`*"]
  5. G=1; M=1; c=1; rs=2 G M/c^2;
  6. wp=MachinePrecision;
  7. para=20; pstep=1/2;
  8. j[v_]:=Sqrt[1-v^2]; J=j[v0];
  9. k[r_]:=Sqrt[1-rs/r];
  10. к=k[r0];
  11. r0=Sqrt[2] 10;
  12. θ0=0;
  13. тmax=200;
  14. Ф=β Pi/180;
  15. vr0=v0 Sin[Ф] к/J;
  16. vθ0=v0/r0 Cos[Ф]/J;
  17. v0=9999/10000;
  18. Table[Subscript[sol, β]=
  19. NDSolve[{r''[t]==-((G M)/r[t]^2)+r[t] θ'[t]^2-(3 G M)/c^2 θ'[t]^2,
  20. r'[0]==vr0,
  21. r[0]==r0,
  22. θ''[t]==-((2 r'[t] θ'[t])/r[t]),
  23. θ'[0]==vθ0,
  24. θ[0]==θ0,
  25. τ'[t]==Sqrt[c^2 r[t]+r[t] r'[t]^2-c^2 rs+r[t]^3 θ'[t]^2-r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t]-rs] Sqrt[1-rs/r[t]]),
  26. τ[0]==0,
  27. cl'[t]==((r'[t]/k[r[t]])^2+(θ'[t] r[t])^2)/c^2,
  28. cl[0]==0}, {r, θ, τ, cl}, {t, 0, тmax},
  29. MaxSteps-> Infinity, Method-> Automatic, WorkingPrecision-> wp,
  30. InterpolationOrder-> All], {β, 170, 370, 1}];
  31. t[Χ_, β_]:=Quiet[ξ /. FindRoot[Evaluate[τ[ξ] /. Subscript[sol, β]][[1]]-Χ, {ξ, 0},
  32. WorkingPrecision-> wp, Method-> Automatic]];
  33. Τ[β_]:=Quiet[t[time, β]];
  34. x[t_, β_]:=(Sin[Evaluate[θ[t] /. Subscript[sol, β]]] Evaluate[r[t] /. Subscript[sol, β]])[[1]]
  35. y[t_, β_]:=(Cos[Evaluate[θ[t] /. Subscript[sol, β]]] Evaluate[r[t] /. Subscript[sol, β]])[[1]]
  36. R[t_, β_]:=Evaluate[r[t] /. Subscript[sol, β]][[1]];
  37. γ[t_, β_]:=Evaluate[τ'[t] /. Subscript[sol, β]][[1]];
  38. и[t_, β_]:=Evaluate[τ[t] /. Subscript[sol, β]][[1]];
  39. crθ[t_, β_]:=Evaluate[cl'[t] /. Subscript[sol, β]][[1]];
  40. vrθ[t_, β_]:=crθ[t]/Sqrt[1+crθ[t]^2];
  41. clr[t_, β_]:=Evaluate[r'[t] /. Subscript[sol, β]][[1]];
  42. clθ[t_, β_]:=R[t] Evaluate[θ'[t] /. Subscript[sol, β]][[1]];
  43. s[text_]:=Style[text, FontSize-> font]; font=11;
  44. PR=Sqrt[2] 10;
  45. plot=Do[Print[Rasterize[Grid[{{Show[
  46. Graphics[{{LightGray, Disk[{0, 0}, rs]}}, Frame-> True,
  47. ImageSize-> 400, PlotRange-> PR, ImagePadding-> 1],
  48. Table[{Graphics[{Red, Point[{x[Τ[β], β], y[Τ[β], β]}]}],
  49. ParametricPlot[{x[ε, β], y[ε, β]}, {ε, 0, Τ[β]},
  50. PlotStyle-> {Thickness[0.001], Red}]}, {β, 170, 370, 1}]]},
  51. {Grid[{{s["t"], "=", s[N[time]], s[" GM/c³"]}},
  52. Alignment-> Left, Spacings-> {0, 1/2}]}}, Alignment-> Left]]],
  53. {time, 10, 50, 10}]
RAW Paste Data