Not able to get plot graph from differential equation

Question

I'm trying to recreate the plots in Figure 2 in https://arxiv.org/pdf/1302.1705.pdf While I'm able to get the embedding diagram, I'm not able to plot the trajectory for the particle orbiting the wormhole. The equations I've used from the paper are equation 3 for the differential equation dr/dΦ, equation 12 for z(r), and x = rcos(phi), y = rsin(phi). Here's the code that I've used to solve the differential equation numerically and then trying to plot a 3D graph for the trajectory.

rulesDS = {L -> Sqrt[4.5], Energy -> Sqrt[0.949], b -> 1.01};
z[r_] := 2 Sqrt[r - 1] - (2*Sqrt[b - 1] /. rulesDS);
RevolutionPlot3D[{{r, z[r]}, {r, -z[r]}}, {r, 1.01, 
  2 }] 
(*this gives the embedding diagram*)

eq1 = D[r[phi], 
     phi]^2 == (r[
        phi]^4 (Energy^2 - (1 + (L/r[phi])^2) (1 - 1/r[phi])))/L^2;
eq1n = eq1 /. rulesDS;
solDS1N = 
  NDSolve[{eq1n, r[0] == 1.01} /. rulesDS, r, {phi, 0, 2 Pi}, 
      Method -> {"EquationSimplification" -> "Residual", 
        "DAEInitialization" -> {"Collocation", 
          "DefaultStartingValue" -> #}}][[1]] & /@ {1, -1};
solDS1N[[1]] 
(*this solves our differential equation and takes the \
real part of the equation*)

Show[RevolutionPlot3D[{{r, z[r]}, {r, -z[r]}}, {r, 1.01, 2}, 
  PlotStyle -> 
   Directive[Magenta, Opacity[0.4], Specularity[White, 10]], 
  MeshStyle -> Blue, PlotPoints -> 30, Boxed -> False, 
  PlotRange -> All, Axes -> False], 
 ParametricPlot3D[
  Evaluate[{Cos[phi] r[phi], Sin[phi] r[phi], 
     4 Sqrt[-1 + h[phi]/2]} /.     solDS1N[[1]]], {phi, 0, 2 Pi}, 
  PlotRange -> All, PlotStyle -> Red], 
 ParametricPlot3D[
  Evaluate[{Cos[phi] r[phi], 
     Sin[phi] r[phi], -4 Sqrt[-1 + h[phi]/2]} /.     
    solDS1N[[1]]], {phi, 0, 2 Pi}, PlotRange -> All, 
  PlotStyle -> Red]] 
(*only get the embedding diagram*)
```

Answer 1

The problem seems to be with the z-coordinate in the ParametricPlot3D

Show[RevolutionPlot3D[{{r, z[r]}, {r, -z[r]}}, {r, 1.01, 2}, 
  PlotStyle -> Directive[Magenta, Opacity[0.4], Specularity[White, 10]], 
  MeshStyle -> Blue, PlotPoints -> 30, Boxed -> False, 
  PlotRange -> All, Axes -> False], 
 ParametricPlot3D[{Cos[phi] r[phi], Sin[phi] r[phi], z[r[phi]]} /. 
   solDS1N[[1]], {phi, 0, 2 Pi}, PlotStyle -> Red], 
 ParametricPlot3D[{Cos[phi] r[phi], Sin[phi] r[phi], -z[r[phi]]} /. 
   solDS1N[[1]], {phi, 0, 2 Pi}, PlotStyle -> Red]]