Am I seeing a particle orbiting a Morris-Thorne Wormhole?

I tried to calculate the time-like geodesics of the Morris-Thorne Wormhole $ [1]$ , $ [2]$ , $ [3]$ for redshift function $ \Phi(r) = 0$ and $ b(r) = \sqrt{r_{0} r}$ . But I don’t know for sure if all the plots that I made are correct.

So, the metric is:

$ $ ds^2 = -e^{2\Phi(r)}dt^2+\Biggr\{ 1-\frac{b(r)}{r}\Biggr\}^{-1}dr^2+r^2[d\theta^2+sin^2(\theta)d\phi^2]\tag{1}$ $

And the plots are:

enter image description here

enter image description here

The various spheres are the 2-spheres due to spherical symmetry of the metric, and the central sphere is the sphere of the radius $ r_{0}$ , therefore is the throat.

My doubt is:

Am I seeing a particle orbiting the throat of a Morris-Thorne Wormhole? In other words is it correct to say that the particle approches from infinity, orbits the throat and returns to the same universe?

$ $ * * * $ $

$ [1]$ LOBO.F.S.N; Exotic solutions in General Relativity,arxiv:0710.4474v1, 2007

$ [2]$ LOBO.F.S.N; Wormholes, Warp Drives and Energy Conditions,Springer, Vol.189, 2017

$ [3]$ VISSER.M; Lorentzian Wormholes, Springer, AIP Press, 1995

Appendix: Mathematica code:

n = Length[coords]; a = 0; r0; \[CapitalPhi][r] = 0; B[r] = (Sqrt[r*r0]);  tt = -(Exp[     2*\[CapitalPhi][       r]   ]); rr = (1)/(1 - ((B[r])/(r))); \[Theta]\[Theta] =   r^2; \[CurlyPhi]\[CurlyPhi] =   r^2*(Sin[\[Theta]]*Sin[\[Theta]]); t\[CurlyPhi] = 0; metric = {{tt,     0, 0, 0}, {0, rr, 0, 0}, {0, 0, \[Theta]\[Theta], 0}, {0, 0,     0, \[CurlyPhi]\[CurlyPhi]}}; metric // MatrixForm inversemetric = Simplify[Inverse[metric]]; inversemetric // MatrixForm christoffel :=   christoffel =    Simplify[Table[     1/2*Sum[inversemetric[[i,          s]]*(D[metric[[s, j]], coords[[k]]] +           D[metric[[s, k]], coords[[j]]] \[Minus]            D[metric[[j, k]], coords[[s]]]), {s, 1, n}], {i, 1, n}, {j,       1, n}, {k, 1, n}]] listchristoffel :=   Table[If[UnsameQ[christoffel[[i, j, k]],       0], {ToString[\[CapitalGamma][i, j, k]],       christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1,      j}] TableForm[    Partition[DeleteCases[Flatten[listchristoffel], Null], 2],     TableSpacing -> {2, 2}] geodesic :=   geodesic =    Simplify[Table[\[Minus]Sum[       christoffel[[i, j, k]] coords[[j]]' coords[[k]]', {j, 1, n}, {k,         1, n}], {i, 1, n}]] listgeodesic :=   Table[{"d/d\[Tau]" ToString[coords[[i]]'], " =", geodesic[[i]]}, {i,     1, n}] TableForm[listgeodesic, TableSpacing -> {2}]  max\[Tau] = 750; ivs = {0, 0, 0.088}; ics = {0, 6.5, \[Pi]/2,    0}; r0 = 1; computeSoln[max\[Tau]i_, ivsi_, icsi_] :=   Block[{ivs, ics, i, \[Chi], tmp, soln}, ics = icsi;   ivs = Join[{\[Chi]}, ivsi];   op1 = Table[coords[[i]] -> ics[[i]], {i, 0, n}];   tm = metric /. op1;   tmp = ivs.(tm.ivs); \[Chi]slv = Solve[tmp == uinvar, \[Chi]];   ivs[[1]] = Last[\[Chi] /. \[Chi]slv];   op = {Derivative[1][t] -> Derivative[1][t][\[Tau]],      Derivative[1][r] -> Derivative[1][r][\[Tau]],      Derivative[1][\[Theta]] -> Derivative[1][\[Theta]][\[Tau]],      Derivative[1][\[CurlyPhi]] -> Derivative[1][\[CurlyPhi]][\[Tau]],      t -> t[\[Tau]],      r -> r[\[Tau]], \[Theta] -> \[Theta][\[Tau]], \[CurlyPhi] -> \ \[CurlyPhi][\[Tau]]};   deq = Table[     coords[[i]]''[\[Tau]] == Simplify[geodesic[[i]] /. op], {i, 1,       n}]; deq =     Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}],      Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];   soln = NDSolve[deq, coords, {\[Tau], 0, max\[Tau]i}]; soln] uinvar = \[Minus]1; sphslnToCartsln[soln_] :=   Block[{xs, ys, zs},    xs = r[\[Tau]] Sin[\[Theta][\[Tau]]] Cos[\[CurlyPhi][\[Tau]]] /.      soln; ys =     r[\[Tau]] Sin[\[Theta][\[Tau]]] Sin[\[CurlyPhi][\[Tau]]] /. soln;   zs = r[\[Tau]] Cos[\[Theta][\[Tau]]] /. soln; {xs, ys, zs}] udotu[solni_, \[Tau]val_] :=   Block[{x\[Alpha], u\[Alpha]},    x\[Alpha] =     Table[coords[[i]][\[Tau]] /. solni, {i, 1, n}] // Flatten;   u\[Alpha] = D[x\[Alpha], \[Tau]];   x\[Alpha] = x\[Alpha] /. \[Tau] -> \[Tau]val;   u\[Alpha] = u\[Alpha] /. \[Tau] -> \[Tau]val;   u\[Alpha].((metric /.         Table[coords[[i]] -> x\[Alpha][[i]], {i, 1, n}]).u\[Alpha])] coordlist[\[Tau]in_] :=   Table[ToString[coords[[i]]] <>     " = " <> {ToString[coords[[i]][\[Tau]in] /. soln // First]}, {i, 1,     n}]  soln = computeSoln[max\[Tau], ivs, ics];  xyzsoln = sphslnToCartsln[soln]; Join[{"Final Coordinates:"}, coordlist[max\[Tau]]] // TableForm Join[{{"", "", "", "u.u values"}},    Table[{"\[Tau]=", ToString[i], "->", udotu[soln, i]}, {i, 0,      max\[Tau], max\[Tau]/5}]] // TableForm {Plot[Evaluate[    Table[coords[[i]][\[Tau]] /. soln, {i, 1, n}]], {\[Tau], 0,     max\[Tau]}, AxesLabel -> {"\[Tau]", "Coordinate"},    PlotLegends -> {"t", "r", "\[Theta]", "\[CurlyPhi]"},    PlotRange -> {0, 30}],   Show[ParametricPlot[    Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], {\[Tau], 0,      max\[Tau]}, AspectRatio -> 1, PlotStyle -> Gray],    Graphics[{Red, Circle[{0, 0}, 2]}]]}  Shape FUNCTION;  (1)/(((r)/((Sqrt[r*r0]))) - 1)  Integrate[(1)/(((r)/((Sqrt[r*r0]))) - 1), r]  max\[Tau] = 2000; ivs = {\[Minus]0.08, .035, .0359}; ics = {0, 2, Pi/4, 0.2}; xyzsoln = sphslnToCartsln[computeSoln[max\[Tau], ivs, ics]]; angle = ParametricPlot3D[    Evaluate[Re[xyzsoln] // Flatten], {\[Tau], 0, max\[Tau]},     AxesLabel -> {x, y, z}, DisplayFunction -> Identity]; sphhoriz =    SphericalPlot3D[{1, 2, 3, 4}, {\[Theta], 0, Pi}, {\[Phi], 0, 2*Pi},     PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}},     BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2,     Mesh -> None, PlotStyle -> Opacity[0.2]]; sphhoriz2 =    SphericalPlot3D[{1, 2, 3, 4}, {\[Theta], 0, Pi}, {\[Phi],      0, (3/2)*Pi}, PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}},     BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2,     Mesh -> None, PlotStyle -> Opacity[0.3]]; Show[angle, sphhoriz, DisplayFunction -> $  DisplayFunction,   PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}] Show[angle, sphhoriz2 , DisplayFunction -> $  DisplayFunction,   PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}]