# 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