
eqr3=r'[tau]^2 + (L^2/r[tau]^2 +1) (1-2/r[tau]) -e^2
eqp3=p'[tau]-L/r[tau]^2
eqct3=ct'[tau]-e/(1-2/r[tau])

Mercury:
Eccentricity	0.205630
Semi-major axis 57,909,050 km -> 39217076.1
Aphelion 69,816,900 km        -> 47281291.6   
Perihelion 46,001,200 km      -> 31152860.6
L=sqrt(GM a (1-e^2))          -> 6128.5
e=Sqrt[(L^2/r0^2 +1)(1-2/r0)] -> .9999999872504522
T=0.240846 year =             ->  1543063000236.21
* ? gn*Msun/c^2
    1476.628442581272    
* set rb=69816900e3 ra=46001200e3 a=57909050e3
* ? rb/(gn*Msun/c^2)
    47281291.61453379    
* ? ra/(gn*Msun/c^2)
    31152860.57986665    
? 0.240846*LY/1476.628442581272
    1543063000236.216    
* ? Sqrt((L^2/ra^2 +1)*(1-2/ra))
    .9999999872504522    
* ? Sqrt((L^2/rb^2 +1)*(1-2/rb))
    .9999999872504521    
? sqrt(a/(gn*Msun/c^2)*(1-.205630^2))
    6128.526179939347    

L=6128.5
e=.9999999872504522
T= 1543063000236.21
rA=31152860.6
rB=47281291.6
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,rA*.95,rB*1.05}]
r0=39217076.1
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,10 T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,.99*T},AspectRatio->Automatic]

finite proper time; infinite coordinate time
L=0
Integrate[1/Sqrt[1/r-1/rmax],{r,0,rmax} ]/Sqrt[2]
Simplify[% , rmax>0]        
               3/2
         Pi rmax
Out[54]= ----------
             2

whip around:
r0=25.6
L=3.75
e=.9704
T=1000
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]

Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,3,30}]

solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,10 T}]

ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,T},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,9 T,10 T},AspectRatio->Automatic]

sucked in:
r0=25.6
L=3.3
e=.9704
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,2,30}]
dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,150}]
Plot[Evaluate[r[tau]  /. Last[solution]],{tau,0,141.249}]
Plot[Evaluate[ct[tau]  /. Last[solution]],{tau,0,141.249}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,141.249},AspectRatio->Automatic]

spiral takeoff:
L=5
r0=L^2/2(1-Sqrt[1-12/L^2])
e=Sqrt[(L^2/r0^2 +1) (1-2/r0)]
r0=r0+.00001
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,r0,30},PlotRange->All]

solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,1000}]


ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,100},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,1000},AspectRatio->Automatic]

Plot[Evaluate[r[tau] /. Last[solution]],{tau,0,100},AspectRatio->Automatic]
Plot[Evaluate[p[tau] /. Last[solution]],{tau,0,100},AspectRatio->Automatic]

large scatter:
L=5
r0=L^2/2(1-Sqrt[1-12/L^2])
e=Sqrt[(L^2/r0^2 +1) (1-2/r0)]
r0=30
dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,500}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,200},AspectRatio->Automatic]
Plot[Evaluate[r[tau] /. Last[solution]],{tau,0,200},AspectRatio->Automatic]
Plot[Evaluate[p[tau] /. Last[solution]],{tau,0,200},AspectRatio->Automatic]


Oumuamua
Perihelion	0.25534  AU  ->  25868606.6
Semi-major axis −1.2798  AU  -> 129657095.4
Eccentricity	1.19951 
L=sqrt(GM a (e^2-1))         -> 7543.
AU                           ->101310435.6
interstellar speed 26.33 km/s
e=sqrt(1+(v/c)^2)            ->1.000000003856828

? 1.2798*au/(gn*msun/c^2)
* ? sqrt(res*(1.19951^2-1))
    7542.988557846278    
* ? 0.25534*au/(gn*msun/c^2)
    25868606.61966126    
* ? au/(gn*msun/c^2)
    101310435.5747680    
*  ? sqrt((26.33e3/c)^2+1)
    1.000000003856828  
* ? ly/(gn*msun/c^2)
    6406845038888.813    
* ? sqrt(1.2798*au/(gn*Msun/c^2)*(1.19951^2-1))
    7542.988557846278    

L=7543.
e=1.000000003856828
r0=101310435.5*2
dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
T=6406845038888.8

solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,.5 T},AspectRatio->Automatic]
Plot[Evaluate[r[tau] /. Last[solution]],{tau,.1 T,.3 T}]
Plot[Evaluate[r[tau] /. Last[solution]],{tau,0, .5 T}]

min distance: 2.6e6 (check)
FindMinValue[Evaluate[r[tau] /. Last[solution]],{tau,.2T}]                                                

                    7
Out[123]= 2.58677 10

reduce L for sharper curve:
L=3000.
e=1.000000003856828
r0=101310435.5*2
dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]

T=6406845038888.8
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,.5 T},AspectRatio->Automatic]



for circular orbit: r=L^2(1+-Sqrt[1-12/L^2])/2
e^2=(1 +L^2/r0^2) (1-2/r0) 


L=3
r0=9/2
e=Sqrt[(1 +L^2/r0^2) (1-2/r0) ]
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,3,10},PlotRange->All]

L=Sqrt[12]
r0=6
e=L (1-2/r0)/Sqrt[r0]
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,3,10},PlotRange->All]

L=4
rA=L^2(1-Sqrt[1-12/L^2])/2
rB=L^2(1+Sqrt[1-12/L^2])/2
e=L (1-2/rB)/Sqrt[rB]
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,.9 rA,1.2 rB},PlotRange->All]

r0=rB
e=L (1-2/r0)/Sqrt[r0]+.001
dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]

T=200
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,10*T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,T},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,10 T},AspectRatio->Automatic]

r0=rA
e=L (1-2/r0)/Sqrt[r0]+.000001
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,10*T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,T},AspectRatio->Automatic]

dr0=-Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
T=32
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,31.22},AspectRatio->Automatic]


L=8
rA=L^2(1-Sqrt[1-12/L^2])/2
rB=L^2(1+Sqrt[1-12/L^2])/2
e=L (1-2/rB)/Sqrt[rB]
Plot[(L^2/r^2 +1) (1-2/r) -e^2,{r,.75 rA,1.2 rB},PlotRange->All]
r0=rB
e=L (1-2/r0)/Sqrt[r0]+.0001
T=2 Pi r0^2/L
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,10 T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,T},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,10 T},AspectRatio->Automatic]


L=16
r0=L^2(1+Sqrt[1-12/L^2])/2
e=L (1-2/r0)/Sqrt[r0]+.0001
T=2 Pi r0^2/L
dr0=Sqrt[-(L^2/r0^2 +1) (1-2/r0)+e^2]
solution =NDSolve[{D[eqr3,tau]==0,r'[0]==dr0,eqp3==0,eqct3==0,r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{tau,0,100 T}]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,T},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,0,10 T},AspectRatio->Automatic]
ParametricPlot[Evaluate[{r[tau] Cos[p[tau]], r[tau] Sin[p[tau]]} /. Last[solution]],{tau,99 T,100 T},AspectRatio->Automatic]
Plot[Evaluate[r[tau] /. Last[solution]],{tau,0,9.9 T}]
Plot[Evaluate[p[tau] /. Last[solution]],{tau,0,9.9 T}]

find perihelion advance
pp[tau_]=Evaluate[p[tau] /. Last[solution]]
rr[tau_]=Evaluate[r[tau] /. Last[solution]]
FindRoot[rr[x]==rr[0],{x, 9.9 T}]
pp[x]/(2 Pi) /. First[%]
Out[258]= 9.1087

Plot[{(rr[x]-rr[0])/10,Mod[pp[x], 2 Pi]},{x,0,247472.}]
FindMaximum[rr[x],{x, 10000}]

ps=ConstantArray[0,{3,9}]
Do[ps[[1,i]]=x /. Last[FindMaximum[rr[x],{x, -18668.8+ 247472./9 i},Method->PrincipalAxis]];
ps[[2,i]]=First[FindMaximum[rr[x],{x, -18668.8+ 247472./9 i},Method->PrincipalAxis]];
ps[[3,i]]=Mod[pp[ps[[1,i]]], 2 Pi], {i,9}]


--classic tests: perihelion advance Mercury

sigma[r_]=1-2/r
e2=(rb^2-ra^2)/(rb^2/sigma[rb]-ra^2/sigma[ra])
L2=(sigma[rb]-sigma[ra])/(sigma[ra]/ra^2-sigma[rb]/rb^2)
Integrate[1/(r^2 Sqrt[e2-(L2/r^2+1)(1-2/r)]),{r,ra,rb}]
* ? gn*Msun/c^2
    1476.628442581272    
* set rb=69816900e3 ra=46001200e3
* ? rb/(gn*Msun/c^2)
    47281291.61453379    
* ? ra/(gn*Msun/c^2)
    31152860.57986665    
Sqrt[L2] %[[1]] -Pi /. {ra->31152860.6, rb->47281291.6}                              

                  -7
Out[9]= 2.50934 10

In[10]:= % 2 180/Pi 3600                                                        

Out[10]= 0.103518
? 415.2030758243857 * 0.103518
    42.98099200318875    arcsec per century

--classic test: light bending near Sun

NoJoy:
Integrate[1/Sqrt[(1-2/ro) x^2-(1-2/(ro x))]/x,{x,1,Infinity}]
Integrate[1/( Sqrt[1-2/x] Sqrt[(1-2/ro)/(1-2/x) (x/ro)^2-1] x),{x,ro,Infinity}]
Integrate[1/(x*Sqrt[-1 + x^2]) + (-1 + x^3)/(ro*x^2*(-1 + x^2)^(3/2)) ,{x,1,Infinity}]
Integrate[1/Sqrt[(1-2/ro)/ro^2-u^2(1-2 u)],{u,0,1/ro}]

try series expand
Series[1/Sqrt[(1-2/ro) x^2-(1-2/(ro x))]/x,{ro,Infinity,1}]
Integrate[SeriesCoefficient[%,0],{x,1,Infinity}]
Integrate[SeriesCoefficient[%%,1],{x,1,Infinity}]

* ? 4/(Rsun/1476.628)*180/pi*3600
    1.751193838258085    arcsec


