Solve the equations for the velocities

In[96]:=

  
  eqrules = Solve[eqlist,{rdot[r],thdot[r],tdot[r]}]

Out[96]=

                            2
                      2    J
  {{rdot[r] -> -Sqrt[E  - -----], tdot[r] -> E, 
                           2  2
                          c  r
   
                                                 2
                   J                       2    J
     thdot[r] -> -----}, {rdot[r] -> Sqrt[E  - -----], 
                  2  2                          2  2
                 c  r                          c  r
   
                                 J
     tdot[r] -> E, thdot[r] -> -----}}
                                2  2
                               c  r

Now we face a sign choice for our path. The derivative rdot[r] can be positive (travelling away from the central mass as time increases) or negative (travelling towards from the central mass as time increases).
So we'll make two functions, one for incoming and one for outgoing geodesics:

In[97]:=

  
  rdotin[r_] = rdot[r] /. eqrules[[1]][[1]]

Out[97]=

               2
         2    J
  -Sqrt[E  - -----]
              2  2
             c  r

In[98]:=

  
  rdotout[r_] = rdot[r] /. eqrules[[2]][[1]]

Out[98]=

              2
        2    J
  Sqrt[E  - -----]
             2  2
            c  r

In[99]:=

  
  thdot[r_] = thdot[r] /. eqrules[[1]][[3]]

Out[99]=

    J
  -----
   2  2
  c  r

In[100]:=

  
  tdot[r_] = tdot[r] /. eqrules[[1]][[2]]

Out[100]=

  E

The energy and angular momentum of a photon's trajectory around the central point mass can be absorbed into the "impact parameter", which measure the distance of closest approach to the point mass. The impact parameter is the ratio b = J/E, or angular momentum per unit energy. Notice that a photon trajectory with zero angular momentum relative to the central point mass has zero impact parameter. This is correct, because it's just a direct path to the point mass! The distance of closest approach is zero.

We institute this reparametrization as follows:

In[101]:=

  
  paramrules = {E -> 1, J -> b}

Out[101]=

  {E -> 1, J -> b}

In[102]:=

  
  rdotout[r_] = rdotout[r] /. paramrules

Out[102]=

             2
            b
  Sqrt[1 - -----]
            2  2
           c  r

In[103]:=

  
  rdotin[r_] = rdotin[r] /. paramrules

Out[103]=

              2
             b
  -Sqrt[1 - -----]
             2  2
            c  r

In[104]:=

  
  thdot[r_] = thdot[r] /. paramrules

Out[104]=

    b
  -----
   2  2
  c  r

In[105]:=

  
  
  tdot[r_] = tdot[r] /. paramrules

Out[105]=

  1

Up to