题目还是这个：
平面内两质点O和P，距离为r0，O固定，质量为M，P初速度v0，方向垂直两点连线。P仅受万有引力作用。仅考虑牛顿力学（相对论一边去），求P轨迹
转换成数学语言
P''=-GM*P/|P|3
P(0)=(r0,0)
P'(0)=(0,v0)
求P(t)

设P=(rcosθ,rsinθ)
则|P|=r
P''=((rcosθ)'',(rsinθ)'')
=(r''cosθ-2r'θ'sinθ-rθ''sinθ-rθ'2cosθ,
r''sinθ+2r'θ'cosθ+rθ''cosθ-rθ'2sinθ)
原方程变为：
r''cosθ-2r'θ'sinθ-rθ''sinθ-rθ'2cosθ=-GMcosθ/r2
r''sinθ+2r'θ'cosθ+rθ''cosθ-rθ'2sinθ=-GMsinθ/r2
↓
r''-2r'θ'tanθ-rθ''tanθ-rθ'2=-GM/r2……①
r''+2r'θ'cotθ+rθ''cotθ-rθ'2=-GM/r2……②
②-①：
(2r'θ'+rθ'')(tanθ+cotθ)=0
∵tanθ+cotθ≠0
∴2r'θ'+rθ''=0……③
③代入①：
r''-rθ'2=-GM/r2

r''-rθ'2=-GM/r2……④
④整理可得
θ'=±√(GM/r3+r''/r)⑤
θ''=±[-3GMr'/r4+(r'''r-r''r')/r2]/√(GM/r3+r''/r)⑥
⑤⑥代入③
±r'√(GM/r3+r''/r)±r[-3GMr'/r4+(r'''r-r''r')/r2]/√(GM/r3+r''/r)=0
4r'(GM/r3+r''/r)=3GMr'/r3-r'''+r''r'/r
GMr'/r3+r'''+3r''r'/r=0
GMr'+r3r'''+3r2r'r''=0
r3r'''+3r2r'r''=-GMr'
(r3r'')'=(-GMr)'
r3r''=-GMr + constA

r3r''=-GMr + constA……⑦
考虑constA：
constA=r(0)3r''(0)+GMr(0)=r03r''(0)+GMr0
至于r''(0)，考虑④式
r''(0)=-GM/r(0)2+r(0)θ'(0)2=-GM/r02+v02/r0
∴constA=-GMr0+v02r02+GMr0=v02r02
因此⑦式为
r3r''=-GMr+v02r02

r3r''=-GMr+v02r02
r''=-GM/r2+v02r02/r3
dr/dt d(dr/dt) = dr(-GM/r2+v02r02/r3)
r'2/2=GM/r-v02r02/(2r2)+ constB

san值：99→0

r'2/2=GM/r-v02r02/(2r2)+ constB ……⑧
考虑constB：
constB=r'(0)2/2-GM/r(0)+v02r02/(2r(0)2)=-GM/r0+v02/2
∴⑧式为
r'2/2=GM/r-v02r02/(2r2)-GM/r0+v02/2
r'2=2GM(1/r-1/r0)+v02(1-r02/r2)

lz作死

好吧，好像写成这样更好：
r'2=2GM/r-v02r02/r2-2GM/r0+v02
这个方程的解：
Ⅰ.当v02<2GM/r0时
引进变参φ∈[0,∞）
r与t的关系满足参数方程
t=[GMφ+(GM-v02r0)sinφ]/(2GM/r0-v02)^(3/2)
r=[GM+(GM-v02r0)cosφ]/(2GM/r0-v02)
Ⅱ.当v02>2GM/r0时
引进变参ψ∈[0,∞）
r与t的关系满足参数方程
t=[GMψ+(GM-v02r0)sinhψ]/(v02-2GM/r0)^(3/2)
r=[GM+(GM-v02r0)coshψ]/(2GM/r0-v02)
事实上ⅠⅡ两种情况的结果可用φ=iψ相互转化
Ⅲ.当v02=2GM/r0时
r与t的关系满足方程
2/3 * (r-r0)^(3/2) + 2r0*(r-r0)^(1/2) = t√(2GM)

在ⅠⅡ两种情况下继续推导
（之所以一块写是因为两个很相似，且对比起来很有意思）
Ⅰ的用黑色表示
Ⅱ的用红色表示
r'
=(dr/dφ)/(dt/dφ)
=(dr/dψ)/(dt/dψ)
={-(GM-v02r0)sinφ/(2GM/r0-v02)} / {[GM+(GM-v02r0)cosφ]/(2GM/r0-v02)^(3/2)}
={(GM-v02r0)sinhψ/(2GM/r0-v02)} / [GM+(GM-v02r0)coshψ]/(v02-2GM/r0)^(3/2)}
=-(GM-v02r0)sinφ/[GM+(GM-v02r0)cosφ]√(2GM/r0-v02)
=-(GM-v02r0)sinhψ/[GM+(GM-v02r0)coshψ]√(v02-2GM/r0)

r''
=(d(r')/dφ)/(dt/dφ)
=(d(r')/dψ)/(dt/dψ)
={-(GM-v02r0)cosφ*[GM+(GM-v02r0)cosφ]-(GM-v02r0)sinφ*(GM-v02r0)sinφ}/[GM+(GM-v02r0)cosφ]^2*√(2GM/r0-v02) / {[GM+(GM-v02r0)cosφ]/(2GM/r0-v02)^(3/2)}
={-(GM-v02r0)coshψ*[GM+(GM-v02r0)coshψ]+(GM-v02r0)sinhψ*(GM-v02r0)sinhψ}/[GM+(GM-v02r0)coshψ]^2*√(v02-2GM/r0) / {[GM+(GM-v02r0)coshψ]/(v02-2GM/r0)^(3/2)}
=(2GM/r0-v02)^2[-GM(GM-v02r0)cosφ-(GM-v02r0)^2]/[GM+(GM-v02r0)cosφ]^3
=(2GM/r0-v02)^2[-GM(GM-v02r0)coshψ-(GM-v02r0)^2]/[GM+(GM-v02r0)coshψ]^3

根据前面的⑤式
θ'2=GM/r3+r''/r
=GM(2GM/r0-v02)^3/[GM+(GM-v02r0)cosφ]^3 + (2GM/r0-v02)^3[-GM(GM-v02r0)cosφ-(GM-v02r0)^2]/[GM+(GM-v02r0)cosφ]^4
=GM(2GM/r0-v02)^3/[GM+(GM-v02r0)coshψ]^3 + (2GM/r0-v02)^3[-GM(GM-v02r0)cohψ-(GM-v02r0)^2]/[GM+(GM-v02r0)cohψ]^4
=(2GM/r0-v02)^3/[GM+(GM-v02r0)cosφ]^4 *{GM[GM+(GM-v02r0)cosφ] -GM(GM-v02r0)cosφ-(GM-v02r0)^2}
=(2GM/r0-v02)^3/[GM+(GM-v02r0)coshψ]^4 *{GM[GM+(GM-v02r0)coshψ] -GM(GM-v02r0)coshψ-(GM-v02r0)^2}
=(2GM/r0-v02)^3/[GM+(GM-v02r0)cosφ]^4 *(2GMv02r0-v04r02)
=(2GM/r0-v02)^3/[GM+(GM-v02r0)coshψ]^4 *(2GMv02r0-v04r02)
=(2GM/r0-v02)^4/[GM+(GM-v02r0)cosφ]^4 * v02r02
=(2GM/r0-v02)^4/[GM+(GM-v02r0)coshψ]^4 * v02r02
∴θ'
=v0r0(2GM/r0-v02)^2/[GM+(GM-v02r0)cosφ]^2
=v0r0(2GM/r0-v02)^2/[GM+(GM-v02r0)coshψ]^2

θ=∫θ'dt
=∫θ'dφ*dt/dφ
=∫θ'dψ*dt/dψ
=∫v0r0(2GM/r0-v02)^2/[GM+(GM-v02r0)cosφ]^2 * {[GM+(GM-v02r0)cosφ]/(2GM/r0-v02)^(3/2)} *dφ
=∫v0r0(2GM/r0-v02)^2/[GM+(GM-v02r0)coshψ]^2 * {[GM+(GM-v02r0)coshψ]/(v02-2GM/r0)^(3/2)} *dψ
=∫√(2GMv02r0-v04r02)/[GM+(GM-v02r0)cosφ] *dφ
=∫√(v04r02-2GMv02r0)/[GM+(GM-v02r0)coshψ] *dψ

θ
=∫√(2GMv02r0-v04r02)/[GM+(GM-v02r0)cosφ] *dφ
=∫√(v04r02-2GMv02r0)/[GM+(GM-v02r0)coshψ] *dψ
=√(2GMv02r0-v04r02)∫dφ[cos²(φ/2)+sin²(φ/2)]/[(2GM-v02r0)cos²(φ/2)+v02r0sin²(φ/2)]
=√(v04r02-2GMv02r0)∫dψ[cosh²(ψ/2)-sinh²(ψ/2)]/[(2GM-v02r0)cosh²(ψ/2)-v02r0sinh²(ψ/2)]
=2√(2GMv02r0-v04r02)∫d(φ/2)[1+tan²(φ/2)]/[2GM-v02r0+v02r0tan²(φ/2)]
=2√(v04r02-2GMv02r0)∫d(ψ/2)[1-tanh²(ψ/2)]/[2GM-v02r0-v02r0tanh²(ψ/2)]
=2/(2GM-v02r0)*√(2GMv02r0-v04r02)∫d(φ/2)sec²(φ/2)/[1+v02r0/(2GM-v02r0)*tan²(φ/2)]
=2/(2GM-v02r0)*√(v04r02-2GMv02r0)∫d(ψ/2)sech²(ψ/2)/[1-v02r0/(2GM-v02r0)*tanh²(ψ/2)]
=2/√(2GMv02r0-v04r2)*√(2GMv02r0-v04r02)∫d(φ/2)√[v02r0/(2GM-v02r0)]sec²(φ/2)/[1+v02r0/(2GM-v02r0)*tan²(φ/2)]
=2/√(v04r2-2GMv02r0)*√(v04r02-2GMv02r0)∫d(ψ/2)√[v02r0/(v02r0-2GM)]sech²(ψ/2)/[1+v02r0/(v02r0-2GM)*tanh²(ψ/2)]
=2∫d{√[v02r0/(2GM-v02r0)]tan(φ/2)}/{1+[√[v02r0/(2GM-v02r0)]tan(φ/2)]^2}
=2∫d{√[v02r0/(v02r0-2GM)]tanh(ψ/2)}/{1+[√[v02r0/(v02r0-2GM)]tanh(ψ/2)]^2}
=2arctan{√[v02r0/(2GM-v02r0)]tanφ/2} + constC
=2arctan{√[v02r0/(v02r0-2GM)]tanhψ/2} + constC
∵θ(0)=0,φ(0)=0,ψ(0)=0,constC=constD=0
tan(θ/2)=
√[v02r0/(2GM-v02r0)]tan(φ/2)
√[v02r0/(v02r0-2GM)]tanh(ψ/2)