This is some matterial that could be a helpfull companion to the paper, 
"An all-purpose metric for the exterior of any kind of rotating neutron star", by G. Pappas and T. A. Apostolatos.

metric functions for the spacetimes used

Two - Soliton functions

<< Algebra`ReIm`
Off[General::"spell1"]
Off[Plot::"plnr"]

Clear[z, M, ρ, a, b, k]
z /: Im[z] = 0;
M /: Im[M] = 0;
a /: Im[a] = 0;
ρ /: Im[ρ] = 0;
b /: Im[b] = 0;
k /: Im[k] = 0;

\(\(A = \(-16\)\ d\ k\ M\^2\ \((\(r\_-\)\ \(r\_+\) + \(R\_-\)\ \(R\_+\))\) - \((\(-4\)\ I\ \((a\^2\ b - a\ \((b\^2 + d + k)\) + 
                      b\ \((d + k - 
                            M\^2)\))\)\ \(κ\_-\)\ \((\(R\_-\)\ \(r\_+\) - \(r\_-\)\ \(R\_+\))\) + \((4\ \((a - b)\)\^2\ \((a\ b - d - k)\) - 
                      4\ \((\((a - 2\ b)\)\ b - 
                            d)\)\ M\^2)\)\ \((\(R\_-\)\ \(r\_+\) + \(r\_-\)\ \(R\_+\))\))\)\ κ\_+\%2 + κ\_-\%2\ \((\((4\ \((a - b)\)\^2\ \((a\ b + d - k)\) - 
                      4\ \((\((a - 2\ b)\)\ b + 
                            d)\)\ M\^2)\)\ \((\(r\_-\)\ \(R\_-\) + \(r\_+\)\ \(R\_+\))\) - 4\ I\ \((\((a - b)\)\ \((a\ b + d - k)\) - 
                      b\ M\^2)\)\ \((\(r\_-\)\ \(R\_-\) - \(r\_+\)\ \(R\_+\))\)\ \(κ\_+\))\);\)\)
\(\(B = \(κ\_-\)\ \(κ\_+\)\ \((\(-2\)\ I\ M\ \((\(R\_-\) - \(R\_+\))\)\ \((\(-\(κ\_-\)\) + \(κ\_+\))\)\ \((\(-a\^2\)\ b + 2\ a\ k + b\ \((b\^2 - 2\ k + M\^2)\) - 
                    b\ \(κ\_-\)\ \(κ\_+\))\) + 
              2\ I\ M\ \((\(r\_-\) - \(r\_+\))\)\ \((\(κ\_-\) + \(κ\_+\))\)\ \((\(-a\^2\)\ b + 2\ a\ k + 
                    b\ \((b\^2 - 2\ k + M\^2)\) + 
                    b\ \(κ\_-\)\ \(κ\_+\))\) + 
              1\/2\ \((8\ d\ M\ \((\(-a\^2\) + b\^2 + 
                          M\^2)\)\ \((\(-\(r\_-\)\) + \(R\_-\) - \(r\_+\) + \(R\_+\))\) - 
                    8\ d\ M\ \(κ\_-\)\ \((\(r\_-\) + \(R\_-\) + \(r\_+\) + \(R\_+\))\)\ \(κ\_+\))\))\);\)\)

\(\(L = 
        4\ I\ d\ \((a + I\ M)\)\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\) - 2\ \((a\^2 - b\^2 - 2\ k)\)\ M\^2 + 
                M\^4)\)\ \((\(r\_-\) + \(R\_-\) + \(r\_+\) + \(R\_+\))\) + 
          4\ I\ \((\(-I\)\ b\^4\ M\^2 + I\ b\^2\ M\^2\ \((d + 2\ k - 
                      M\^2)\) + I\ M\^2\ \((2\ k\ \((d + k)\) + 
                      d\ M\^2)\) + a\^4\ b\ \((M - z)\) + 
                a\^3\ \((b\^2 - d - k)\)\ \((M - z)\) + 
                b\^3\ \((d + k + M\^2)\)\ \((M - z)\) - 
                b\ \((4\ k\ \((d + k)\) - \((d + k)\)\ M\^2 - M\^4)\)\ \((M - 
                      z)\) + a\ \((\(-4\)\ I\ b\ k\ M\^2 + 
                      b\^2\ \((d + 5\ k - M\^2)\)\ \((M - z)\) + \((d + 
                            k)\)\ \((4\ k + M\^2)\)\ \((M - z)\) + 
                      b\^4\ \((\(-M\) + z)\))\) + 
                a\^2\ \((I\ b\^2\ M\^2 - I\ d\ M\^2 + 
                      b\^3\ \((\(-M\) + z)\) + 
                      b\ \((d + 5\ k + 2\ M\^2)\)\ \((\(-M\) + 
                            z)\))\))\)\ \(κ\_-\)\ \((\(-\(R\_-\)\)\ \(r\_+\) + \(r\_-\)\ \(R\_+\))\) + 
          4\ \((I\ b\^5\ M\^2 - I\ b\^3\ \((3\ k\ M\^2 - 2\ M\^4)\) + I\ b\ M\^2\ \((\(-2\)\ k\ \((d + k)\) + 
                      3\ k\ M\^2 + M\^4)\) + 
                a\^3\ \((\(-I\)\ k\ M\^2 - 2\ b\^3\ \((M - z)\) - 
                      2\ b\ \((2\ k + M\^2)\)\ \((M - z)\))\) + 
                a\^2\ \((\(-2\)\ I\ b\^3\ M\^2 - I\ b\ \((3\ k\ M\^2 + 2\ M\^4)\) + 
                      2\ b\^2\ \((d + 5\ k)\)\ \((M - 
                            z)\) + \((4\ k\ \((d + k)\) + \((2\ d + 
                                  k)\)\ M\^2)\)\ \((M - z)\))\) + 
                a\ \((7\ I\ b\^2\ k\ M\^2 + I\ k\ M\^2\ \((2\ \((d + k)\) + M\^2)\) + b\^5\ \((M - z)\) - 
                      2\ b\^3\ \((2\ k - M\^2)\)\ \((M - z)\) - 
                      b\ \((8\ k\ \((d + k)\) - 4\ k\ M\^2 - M\^4)\)\ \((M - 
                            z)\))\) + a\^5\ b\ \((M - z)\) - 
                b\^4\ \((d + k)\)\ \((M - z)\) - 
                M\^2\ \((2\ k\ \((d + k)\) + d\ M\^2)\)\ \((M - z)\) + 
                b\^2\ \((4\ k\ \((d + k)\) - \((2\ d + 
                            3\ k)\)\ M\^2)\)\ \((M - z)\) + 
                a\^4\ \((\(-M\)\ \((d + k - I\ b\ M)\) + \((d + 
                            k)\)\ z)\))\)\ \((\(R\_-\)\ \(r\_+\) + \(r\_-\)\ \(R\_+\))\) - 
          16\ d\ k\ M\^2\ \((I\ a - I\ b + M - 
                z)\)\ \((\(r\_-\)\ \(r\_+\) + \(R\_-\)\ \(R\_+\))\) - 
          4\ \((I\ b\^5\ M\^2 - I\ b\^3\ \((3\ k\ M\^2 - 2\ M\^4)\) + I\ b\ M\^2\ \((2\ \((d - k)\)\ k + 
                      3\ k\ M\^2 + M\^4)\) + 
                a\^3\ \((\(-I\)\ k\ M\^2 - 2\ b\^3\ \((M - z)\) - 
                      2\ b\ \((2\ k + M\^2)\)\ \((M - z)\))\) + 
                a\ \((7\ I\ b\^2\ k\ M\^2 + I\ k\ M\^2\ \((\(-2\)\ d + 2\ k + M\^2)\) + b\^5\ \((M - z)\) - 
                      2\ b\^3\ \((2\ k - M\^2)\)\ \((M - z)\) + 
                      b\ \((8\ \((d - k)\)\ k + 4\ k\ M\^2 + M\^4)\)\ \((M - 
                            z)\))\) + 
                a\^2\ \((\(-2\)\ I\ b\^3\ M\^2 - I\ b\ \((3\ k\ M\^2 + 2\ M\^4)\) - 
                      2\ b\^2\ \((d - 5\ k)\)\ \((M - 
                            z)\) - \((2\ d\ \((2\ k + M\^2)\) - 
                            k\ \((4\ k + M\^2)\))\)\ \((M - z)\))\) + 
                a\^5\ b\ \((M - z)\) + b\^4\ \((d - k)\)\ \((M - z)\) + 
                M\^2\ \((2\ \((d - k)\)\ k + d\ M\^2)\)\ \((M - z)\) - 
                b\^2\ \((4\ d\ k - 4\ k\^2 - 2\ d\ M\^2 + 
                      3\ k\ M\^2)\)\ \((M - z)\) + 
                a\^4\ \((I\ b\ M\^2 + d\ \((M - z)\) + 
                      k\ \((\(-M\) + 
                            z)\))\))\)\ \((\(r\_-\)\ \(R\_-\) + \(r\_+\)\ \(R\_+\))\) - 
          4\ d\ M\ \((I\ a\^3 + 2\ I\ b\ k - 
                a\^2\ M + b\^2\ M + 
                M\^3 - I\ a\ \((b\^2 + 2\ k + 
                      M\^2)\))\)\ \(κ\_-\)\ \((\(r\_-\) - \(R\_-\) + \(r\_+\) - \(R\_+\))\)\ \(κ\_+\) + 
          4\ I\ \((\(-I\)\ b\^4\ M\^2 - I\ b\^2\ M\^2\ \((d - 2\ k + 
                      M\^2)\) - I\ M\^2\ \((2\ \((d - k)\)\ k + 
                      d\ M\^2)\) + a\^4\ b\ \((M - z)\) + 
                a\^3\ \((b\^2 + d - k)\)\ \((M - z)\) - 
                b\^3\ \((d - k - M\^2)\)\ \((M - z)\) + 
                b\ \((4\ \((d - k)\)\ k - \((d - k)\)\ M\^2 + M\^4)\)\ \((M - 
                      z)\) + a\^2\ \((I\ b\^2\ M\^2 + I\ d\ M\^2 + b\ \((d - 5\ k - 2\ M\^2)\)\ \((M - z)\) + 
                      b\^3\ \((\(-M\) + z)\))\) + 
                a\ \((\(-4\)\ I\ b\ k\ M\^2 - \((d - 
                            k)\)\ \((4\ k + M\^2)\)\ \((M - z)\) + 
                      b\^4\ \((\(-M\) + z)\) + 
                      b\^2\ \((d - 5\ k + M\^2)\)\ \((\(-M\) + 
                            z)\))\))\)\ \((\(r\_-\)\ \(R\_-\) - \(r\_+\)\ \(R\_+\))\)\ \(κ\_+\) + 
          4\ M\ \((\(R\_-\) - \(R\_+\))\)\ \((\(-\((\((a - b)\)\ \((a\ b - 
                              k)\)\ \((\((a + b)\)\ \((a\^2 - 
                                    d)\) + \((\(-3\)\ a + 
                                    b)\)\ k)\) + I\ \((a - 
                              b)\)\ \((b\ \((a + b)\)\ \((a\^2 - 
                                    d)\) - \((a\^2 + 3\ a\ b - 2\ d)\)\ k + 
                              2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                              k\ \((2\ b\^2 - d + k)\) - 
                              a\ b\ \((b\^2 + d + 
                                    3\ k)\))\)\ M\^2 - I\ \((2\ a\^2\ b - a\ k - b\ \((b\^2 + d + 2\ k)\))\)\ M\^3 + 
                        a\ b\ M\^4 + I\ b\ M\^5)\)\)\ \(κ\_-\) + \((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 + 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) + I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 + 
                                  d)\) - \((a\^2 + 3\ a\ b + 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 + d + k)\) - 
                            a\ b\ \((b\^2 - d + 
                                  3\ k)\))\)\ M\^2 + I\ \((\(-2\)\ a\^2\ b + b\^3 - b\ d + a\ k + 2\ b\ k)\)\ M\^3 + 
                      a\ b\ M\^4 + I\ b\ M\^5)\)\ \(κ\_+\))\) + 
          4\ M\ \((\(r\_-\) - \(r\_+\))\)\ \((\((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 - 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) + I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 - 
                                  d)\) - \((a\^2 + 3\ a\ b - 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 - d + k)\) - 
                            a\ b\ \((b\^2 + d + 
                                  3\ k)\))\)\ M\^2 - I\ \((2\ a\^2\ b - a\ k - b\ \((b\^2 + d + 2\ k)\))\)\ M\^3 + 
                      a\ b\ M\^4 + I\ b\ M\^5)\)\ \(κ\_-\) + \((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 + 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) + I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 + 
                                  d)\) - \((a\^2 + 3\ a\ b + 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 + d + k)\) - 
                            a\ b\ \((b\^2 - d + 
                                  3\ k)\))\)\ M\^2 + I\ \((\(-2\)\ a\^2\ b + b\^3 - b\ d + a\ k + 2\ b\ k)\)\ M\^3 + 
                      a\ b\ M\^4 + I\ b\ M\^5)\)\ \(κ\_+\))\);\)\)

\(\(Ε = 
        4\ d\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\) - 
                2\ \((a\^2 - b\^2 - 2\ k)\)\ M\^2 + 
                M\^4)\)\ \((I\ a + M - 
                z)\)\ \((\(r\_-\) + \(R\_-\) + \(r\_+\) + \(R\_+\))\) - 
          4\ I\ M\ \((\((a - b)\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b - d - 
                      k)\) - I\ \((\((a - b)\)\ \((a + 
                            b)\)\ \((b\^2 - d)\) + 
                      2\ \((\(-2\)\ a\ b + b\^2 + d)\)\ k + 
                      2\ k\^2)\)\ M - \((a + b)\)\ \((2\ a\ b - b\^2 - d - 
                      k)\)\ M\^2 + I\ \((b\^2 - d)\)\ M\^3 + 
                b\ M\^4)\)\ \(κ\_-\)\ \((\(-\(R\_-\)\)\ \(r\_+\) + \(r\_-\)\ \(R\_+\))\) + 
          4\ M\ \((\(-\((a - b)\)\^2\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b - d - k)\) + I\ \((a - 
                      b)\)\ \((\((a - b)\)\ b\ \((a + b)\)\^2 - \((a\^2 + 
                            4\ a\ b - 3\ b\^2 - 2\ d)\)\ k + 
                      2\ k\^2)\)\ M + \((2\ a\^3\ b + 2\ k\ \((d + k)\) - 
                      a\^2\ \((2\ d + k)\) + b\^2\ \((2\ d + 3\ k)\) - 
                      2\ a\ \((b\^3 + 
                            2\ b\ k)\))\)\ M\^2 - I\ \((2\ a\^2\ b - 2\ b\^3 - a\ k - 3\ b\ k)\)\ M\^3 - \((a\ b - 
                      d)\)\ M\^4 + I\ b\ M\^5)\)\ \((\(R\_-\)\ \(r\_+\) + \(r\_-\)\ \(R\_+\))\) + 
          16\ d\ k\ M\^2\ \((\(-I\)\ a + I\ b + 
                M)\)\ \((\(r\_-\)\ \(r\_+\) + \(R\_-\)\ \(R\_+\))\) + 
          4\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\)\ \((a\ b + d - 
                      k)\) - I\ \((a - 
                      b)\)\ \((\((a - b)\)\ b\ \((a + b)\)\^2 - \((a\^2 + 
                            4\ a\ b - 3\ b\^2 + 2\ d)\)\ k + 
                      2\ k\^2)\)\ M - \((2\ \((a - b)\)\ \((a + 
                            b)\)\ \((a\ b + d)\) - \((a\^2 + 4\ a\ b - 
                            3\ b\^2 + 2\ d)\)\ k + 
                      2\ k\^2)\)\ M\^2 + I\ \((2\ a\^2\ b - 
                      2\ b\^3 - a\ k - 3\ b\ k)\)\ M\^3 + \((a\ b + 
                      d)\)\ M\^4 - I\ b\ M\^5)\)\ \((\(r\_-\)\ \(R\_-\) + \(r\_+\)\ \(R\_+\))\) + 
          4\ d\ M\ \((\(-I\)\ a\^3 - 2\ I\ b\ k + 
                b\^2\ M + 
                M\^3 + I\ a\ \((b\^2 + 2\ k + M\^2)\) - \((b\^2 + 
                      M\^2)\)\ z + 
                a\^2\ \((\(-M\) + 
                      z)\))\)\ \(κ\_-\)\ \((\(r\_-\) - \(R\_-\) + \(r\_+\) - \(R\_+\))\)\ \(κ\_+\) - 
          4\ I\ M\ \((\((a - b)\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b + d - 
                      k)\) - I\ \((\((a - b)\)\ \((a + 
                            b)\)\ \((b\^2 + d)\) - 
                      2\ \((2\ a\ b - b\^2 + d)\)\ k + 2\ k\^2)\)\ M - \((a + 
                      b)\)\ \((2\ a\ b - b\^2 + d - 
                      k)\)\ M\^2 + I\ \((b\^2 + d)\)\ M\^3 + 
                b\ M\^4)\)\ \((\(r\_-\)\ \(R\_-\) - \(r\_+\)\ \(R\_+\))\)\ \(κ\_+\) - 
          4\ M\ \((\(R\_-\) - \(R\_+\))\)\ \((\((a\^5\ b + \((d - 
                            k)\)\ k\ M\^2 + 
                      b\^2\ k\ \((\(-d\) + k - 2\ M\^2)\) - 
                      a\^4\ \((k + I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((d + 3\ k + 
                                  2\ M\^2)\) - I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d - k + M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 + \((d + 3\ k)\)\ M\^2 + 
                                  M\^4)\) - 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) - I\ k\ \((2\ \((d + 
                                        k)\) + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((d + 3\ k + 
                                  M\^2)\) + I\ b\^3\ \((M - 
                                  z)\) + I\ b\ \((d + 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) - I\ b\^3\ \((d + 
                            M\^2)\)\ \((M - 
                            z)\) + I\ b\ \((2\ k\ \((d + 
                                  k)\) - \((d + 2\ k)\)\ M\^2 - 
                            M\^4)\)\ \((M - 
                            z)\))\)\ \(κ\_-\) + \((\(-a\^5\)\ b + 
                      k\ \((d + k)\)\ M\^2 - b\^2\ k\ \((d + k - 2\ M\^2)\) + 
                      a\^4\ \((k + I\ b\ \((M - z)\))\) + 
                      a\^3\ \((b\^3 + 
                            b\ \((\(-d\) + 3\ k + 
                                  2\ M\^2)\) - I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d + k - M\^2)\) + 
                            b\ \((d\ M\^2 + \((k - M\^2)\)\ \((4\ k + 
                                        M\^2)\))\) + 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) + I\ k\ \((\(-2\)\ d + 
                                  2\ k + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((\(-5\)\ b\^2\ k + 
                            k\ \((d - 3\ k - 
                                  M\^2)\) - I\ b\^3\ \((M - 
                                  z)\) + I\ b\ \((d - 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) - I\ b\^3\ \((d - 
                            M\^2)\)\ \((M - 
                            z)\) + I\ b\ \((2\ \((d - 
                                  k)\)\ k - \((d - 2\ k)\)\ M\^2 + 
                            M\^4)\)\ \((M - z)\))\)\ \(κ\_+\))\) + 
          4\ M\ \((\(r\_-\) - \(r\_+\))\)\ \((\((a\^5\ b + \((d - 
                            k)\)\ k\ M\^2 + 
                      b\^2\ k\ \((\(-d\) + k - 2\ M\^2)\) - 
                      a\^4\ \((k + I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((d + 3\ k + 
                                  2\ M\^2)\) - I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d - k + M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 + \((d + 3\ k)\)\ M\^2 + 
                                  M\^4)\) - 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) - I\ k\ \((2\ \((d + 
                                        k)\) + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((d + 3\ k + 
                                  M\^2)\) + I\ b\^3\ \((M - 
                                  z)\) + I\ b\ \((d + 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) - I\ b\^3\ \((d + 
                            M\^2)\)\ \((M - 
                            z)\) + I\ b\ \((2\ k\ \((d + 
                                  k)\) - \((d + 2\ k)\)\ M\^2 - 
                            M\^4)\)\ \((M - 
                            z)\))\)\ \(κ\_-\) + \((a\^5\ b - 
                      k\ \((d + k)\)\ M\^2 + b\^2\ k\ \((d + k - 2\ M\^2)\) - 
                      a\^4\ \((k + I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((\(-d\) + 3\ k + 
                                  2\ M\^2)\) - I\ k\ \((M - 
                                  z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((\(-d\) + 3\ k + 
                                  M\^2)\) + I\ b\^3\ \((M - 
                                  z)\) - I\ b\ \((d - 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) + I\ b\^3\ \((d - 
                            M\^2)\)\ \((M - 
                            z)\) - I\ b\ \((2\ \((d - 
                                  k)\)\ k - \((d - 2\ k)\)\ M\^2 + 
                            M\^4)\)\ \((M - z)\) + 
                      a\ \((\(-b\^3\)\ \((d + k - M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 - \((d - 3\ k)\)\ M\^2 + 
                                  M\^4)\) - 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) + I\ k\ \((\(-2\)\ d + 
                                  2\ k + M\^2)\)\ \((\(-M\) + 
                                  z)\))\))\)\ \(κ\_+\))\);\)\)

\(\(Ko = 16\ d\ κ\_-\%2\ κ\_+\%2;\)\)

Kerr solution

wKerr = ((2 a M ((y^2) - 1) (M kk + (kk^2) x))/(kk ((M^2) ((x^2) - 1) + (a^2) ((y^2) - (x^2))))) //. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            kk -> (Sqrt[(M^2) - (a^2)])};
fKerr = 1 - ((2 M (M + x kk))/((a^2) (y^2) + (((M + x kk))^2))) //. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            kk -> (Sqrt[(M^2) - (a^2)])};
gammaKerr = 
        (1/2) (Log[(((M^2) ((x^2) - 1) + (a^2) ((y^2) - (x^2)))/((x^2) - (y^2)))] - 
                Log[(kk^2)]) //. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            kk -> (Sqrt[(M^2) - (a^2)])};
gttKerr[M_, a_, ρ_, 
          z_] = -fKerr;
gtφKerr[M_, 
          a_, ρ_, z_] = 
        fKerr wKerr;
gφφKerr[
          M_, a_, ρ_, 
          z_] = -fKerr (wKerr^2) + (ρ^2) (fKerr^(-1));
ΩKerr[M_, 
            a_, ρ_] = (1/D)[
                  gφφKerr[M, a, ρ, 
                    0], ρ] (-D[
                    gtφKerr[M, a, ρ, 
                      0], ρ] + √(D[gtφKerr[M, a, ρ, 0], ρ]^2 - 
                      D[gttKerr[M, a, ρ, 0], ρ] D[
                          gφφKerr[M, a, ρ, 
                            0], ρ]));
EtildaKerr[M_, 
          a_, ρ_] = (-gttKerr[M, a, ρ, 0] - 
              gtφKerr[M, a, ρ, 0] ΩKerr[M, 
                  a, ρ])/(√(-gttKerr[M, a, ρ, 0] - 
                  2  gtφKerr[M, a, ρ, 
                      0]  ΩKerr[M, a, ρ] - 
                  gφφKerr[M, a, ρ, 
                      0] ((( ΩKerr[M, a, ρ]))^2)));
DEtildaKerr[M_, 
          a_, ρ_] = -ΩKerr[M, 
              a, ρ] D[EtildaKerr[M, a, ρ], ρ]/D[ΩKerr[M, a, ρ], ρ];

Hartle - Thorne solution 
(on equatorial plane + 
        functions for the transformation to the z - axis)

Am1 = ((15 r (r - 2 Mm) (1 - 3 (u^2)))/(16 (Mm^2))) Log[(r/(r - 2 Mm))];
Am2 = ((15 ((r^2) - 2 (Mm^2)) (3 (u^2) - 1))/(16 (Mm^2))) Log[(r/(r - 2 Mm))];
L1 = 
        80 (Mm^6) + 8 (Mm^4) (r^2) + 10 (Mm^3) (r^3) + 20 (Mm^2) (r^4) - 
          45 Mm (r^5) + 15 (r^6);
p1 = (((8 Mm (r^4) (r - 2 Mm)))^(-1));
W1 = (r - Mm) (16 (Mm^5) + 8 (Mm^4) r - 10 (Mm^2) (r^3) - 
                30 Mm (r^4) + 15 (r^5)) + 
          (u^2) (48 (Mm^6) - 8 (Mm^5) r - 24 (Mm^4) (r^2) - 
                30 (Mm^3) (r^3) - 60 (Mm^2) (r^4) + 135 Mm (r^5) - 
                45 (r^6));
Fm1[Mm_, r_, u_] = -p1 W1 + Am1;
Fm2[Mm_, r_, u_] = 
        5 (r^3) p1 (3 (u^2) - 1) (r - Mm) (2 (Mm^2) + 
                6 Mm r - 3 (r^2)) - Am1;
Gm1[Mm_, r_, u_] = 
        p1 ((L1 - 72 (Mm^5) r) - 
                3 (u^2) (L1 - 56 (Mm^5) r)) - 
          Am1;
Hm1[Mm_, r_, u_] = 
        Am2 + (((8 Mm (r^4)))^(-1)) (1 - 3 (u^2)) (16 (Mm^5) + 
                8 (Mm^4) r - 10 (Mm^2) (r^3) + 15 Mm (r^4) + 
                15 (r^5));
Hm2[Mm_, r_, 
          u_] = -Am2 + (((8 Mm r))^(-1)) 5 (1 - 
                3 (u^2)) (2 (Mm^2) - 3 Mm r - 
                3 (r^2)); 
gttHT[Mm_, jj_, q_, r_, 
          u_] = -(1 - 2 (Mm/r)) (1 + (jj^2) Fm1[Mm, r, u] - 
              q Fm2[Mm, r, u]);
grrHT[Mm_, jj_, q_, r_, 
          u_] = (((1 - 2 (Mm/r)))^(-1)) (1 + (jj^2) Gm1[Mm, r, u] + 
              q Fm2[Mm, r, u]);
gθθHT[Mm_, jj_, q_, r_, u_] = 
        (r^2) (1 + (jj^2) Hm1[Mm, r, u] - 
              q Hm2[Mm, r, u]);
gφφHT[Mm_, jj_, q_, r_, u_] = 
        (r^2) (1 + (jj^2) Hm1[Mm, r, u] - q Hm2[Mm, r, u]) (1 - 
              (u^2));
gtφHT[Mm_, jj_, r_, 
          u_] = -2 ((Mm^2)/r) jj (1 - 
              (u^2));
RcircHTNew[Mm_, jj_, q_, 
          r_] = (Sqrt[g])φφHT[Mm, jj, q, r, 0];
ΩHTNew[Mm_, jj_, q_, 
          r_] = (-D[gtφHT[Mm, jj, r, 0], 
                  r]  + √((((D[gtφHT[Mm, jj, r, 0], r]))^2) - D[gttHT[Mm, jj, q, r, 0], r] D[
                        gφφHT[Mm, jj, q, r, 0], r]))/
          D[gφφHT[Mm, jj, q, r, 0], 
            r];
EtildaHTNew[Mm_, jj_, q_, 
          r_] = (-gttHT[Mm, jj, q, r, 0] - 
              gtφHT[Mm, jj, r, 0] ΩHTNew[Mm, jj, q, 
                  r])/(√(-gttHT[Mm, jj, q, r, 0] - 
                  2 gtφHT[Mm, jj, r, 0] ΩHTNew[Mm, 
                      jj, q, r] - 
                  gφφHT[Mm, jj, q, r, 
                      0] (((ΩHTNew[Mm, jj, q, r]))^2)));
DEtildaHTNew[Mm_, jj_, q_, 
          r_] = -ΩHTNew[Mm, jj, q, 
              r] D[EtildaHTNew[Mm, jj, q, r], r]/D[ΩHTNew[Mm, jj, q, r], r];
ρHT[Mm_, jj_, q_, r_, 
          u_] = (((gtφHT[Mm, jj, r, u]^2 - gttHT[Mm, jj, q, r, u] gφφHT[Mm, jj, q, r, u]))^(1/2));
\(\(gzzHT[Mm_, jj_, q_, r_, 
          u_] = \((D[ρHT[Mm, jj, q, r, u], r]\^2\/grrHT[Mm, jj, q, r, u] + \(\((1 - u\^2)\)\ D[ρHT[Mm, jj, q, r, u], u]\^2\)\/gθθHT[Mm, jj, q, r, u])\)\^\(-1\);\)\)

dzdμHT[Mm_, jj_, q_, r_, 
          u_] = ((Sqrt[1 - (u^2)]) √(((1 - ((2 Mm)/r)) (r^2) ((1/(r^4)) (2 Mm ((jj^2) (16 (Mm^5) + 8 (Mm^4) r - 
                                        10 (Mm^2) (r^3) + 15 Mm (r^4) + 
                                        15 (r^5)) (1 - 3 (u^2)) + 
                                      (r^3) (15 q (r^2) (1 - 
                                        3 (u^2)) + 
                                        10 (Mm^2) q (-1 + 
                                        3 (u^2)) + 
                                        Mm r (8 + 15 q - 
                                        45 q (u^2))))) - 
                            15 ((jj^2) + q) (2 (Mm^2) - 
                                  (r^2)) (-1 + 3 (u^2)) Log[
                                (r/(-2 Mm + r))]))/((Mm^2) ((1/(8 Mm (2 Mm - r) (r^4))) ((jj^2) (80 (Mm^6) (-1 + 
                                        3 (u^2)) + 
                                        8 (Mm^4) (r^2) (-1 + 
                                        3 (u^2)) + 
                                        10 (Mm^3) (r^3) (-1 + 
                                        3 (u^2)) + 
                                        20 (Mm^2) (r^4) (-1 + 
                                        3 (u^2)) - 
                                        45 Mm (r^5) (-1 + 
                                        3 (u^2)) + 
                                        15 (r^6) (-1 + 3 (u^2)) - 
                                        24 (Mm^5) r (-3 + 
                                        7 (u^2))) + 
                                  (r^3) (10 (Mm^3) q (-1 + 
                                        3 (u^2)) + 
                                        15 q (r^3) (-1 + 3 (u^2)) + 
                                        4 (Mm^2) r (4 - 5 q + 
                                        15 q (u^2)) + 
                                        Mm (r^2) (-8 + 
                                        45 q (1 - 
                                        3 (u^2))))) - ((15 ((jj^2) + q) (2 Mm - r) r (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))])/(16 (Mm^2)))))) ((1/((2 Mm - 
                              r) (r^8))) (4 (Mm^2) (-(jj^4) (-1 + 3 (u^2)) (675 (r^12) (1 - 3 (u^2)) + 
                                  1792 (Mm^12) (-1 + 3 (u^2)) - 
                                  300 (Mm^6) (r^6) (-1 + 3 (u^2)) + 
                                  220 (Mm^5) (r^7) (-1 + 3 (u^2)) + 
                                  1380 (Mm^4) (r^8) (-1 + 3 (u^2)) - 
                                  3000 (Mm^3) (r^9) (-1 + 3 (u^2)) - 
                                  675 (Mm^2) (r^10) (-1 + 3 (u^2)) + 
                                  2250 Mm (r^11) (-1 + 3 (u^2)) + 
                                  96 (Mm^8) (r^4) (-3 + 4 (u^2)) - 
                                  128 (Mm^11) r (-7 + 9 (u^2)) + 
                                  160 (Mm^7) (r^5) (-8 + 27 (u^2)) - 
                                  64 (Mm^10) (r^2) (-15 + 
                                        47 (u^2)) - 
                                  32 (Mm^9) (r^3) (-33 + 
                                        109 (u^2))) + 
                            (r^6) (-100 (Mm^6) (q^2) (((1 - 3 (u^2)))^2) - 100 (Mm^5) (q^2) r (((1 - 3 (u^2)))^2) - 
                                  2250 Mm (q^2) (r^5) (((1 - 3 (u^2)))^2) + 675 (q^2) (r^6) (((1 - 3 (u^2)))^2) + 
                                  (Mm^2) (r^4) (-64 + 
                                        45 q (-1 + 3 (u^2)) (8 - 
                                        15 q + 45 q (u^2))) + 
                                  4 (Mm^4) (r^2) (-32 - 
                                        15 q (-1 + 
                                        3 (u^2)) (-12 - 25 q + 
                                        75 q (u^2))) + 
                                  24 (Mm^3) (r^3) (8 + 
                                        5 q (-1 + 
                                        3 (u^2)) (-9 - 25 q + 
                                        75 q (u^2)))) + 
                            2 (jj^2) (r^3) (640 (Mm^9) q (((1 - 3 (u^2)))^2) - 160 (Mm^5) q (r^4) (((1 - 3 (u^2)))^2) - 
                                  2250 Mm q (r^8) (((1 - 3 (u^2)))^2) + 
                                  675 q (r^9) (((1 - 3 (u^2)))^2) + 
                                  16 (Mm^8) (-12 + 
                                        5 q) r (-1 + 3 (u^2)) - 
                                  360 (Mm^4) (r^5) (-1 + 
                                        3 (u^2)) (-1 - 4 q + 
                                        12 q (u^2)) + 
                                  45 (Mm^2) (r^7) (-1 + 
                                        3 (u^2)) (4 - 15 q + 
                                        45 q (u^2)) + 
                                  4 (Mm^6) (r^3) (25 q (((1 - 3 (u^2)))^2) - 24 (-1 + (u^2))) + 
                                  60 (Mm^3) (r^6) (-1 + 
                                        3 (u^2)) (-9 + 
                                        50 q (-1 + 3 (u^2))) + 
                                  16 (Mm^7) (r^2) (6 (-3 + 
                                        5 (u^2)) - 
                                        5 q (8 - 51 (u^2) + 
                                        81 (u^4)))))) + (1/(r^4)) (15 ((jj^2) + q) (-1 + 
                            3 (u^2)) (4 Mm ((jj^2) ((Mm^8) (24 - 
                                        72 (u^2)) + 
                                        8 (Mm^7) r (3 - 11 (u^2)) + 
                                        16 (Mm^4) (r^4) (1 - 3 (u^2)) + 
                                        30 (Mm^2) (r^6) (1 - 3 (u^2)) + 
                                        16 (Mm^6) (r^2) (-1 + 
                                        3 (u^2)) + 
                                        130 (Mm^3) (r^5) (-1 + 
                                        3 (u^2)) - 
                                        105 Mm (r^7) (-1 + 
                                        3 (u^2)) + 
                                        45 (r^8) (-1 + 3 (u^2))) + 
                                  (r^4) (-105 Mm q (r^3) (-1 + 
                                        3 (u^2)) + 
                                        45 q (r^4) (-1 + 3 (u^2)) + 
                                        (Mm^4) (8 + 20 q - 
                                        60 q (u^2)) + 
                                        6 (Mm^2) (r^2) (2 + 5 q - 
                                        15 q (u^2)) + 
                                        2 (Mm^3) r (-12 + 
                                        65 q (-1 + 
                                        3 (u^2))))) + 
                            15 ((jj^2) + q) (2 Mm - 
                                  r) (r^5) (4 (Mm^3) - 4 (Mm^2) r - 
                                  4 Mm (r^2) + 3 (r^3)) (-1 + 
                                  3 (u^2)) (Log[-(r/Mm)] - 
                                  Log[2 - 
                                      (r/Mm)])) (Log[-(r/Mm)] - 
                            Log[2 - 
                                (r/Mm)]))))/(1024 (Mm^4) √((1/(r^2)) ((1 - (u^2)) (-4 (jj^2) (Mm^4) (-1 + 
                                  (u^2)) - (1/(16 (Mm^2))) ((-1 + ((2 Mm)/r)) (r^4) ((1/(8 Mm (2 Mm - 
                                        r) (r^4))) ((jj^2) (-8 (Mm^5) r (-1 + u) (1 + u) + 8 (Mm^4) (r^2) (1 - 3 (u^2)) + 
                                        10 (Mm^3) (r^3) (1 - 3 (u^2)) + 
                                        20 (Mm^2) (r^4) (1 - 3 (u^2)) + 
                                        15 (r^6) (1 - 3 (u^2)) + 
                                        16 (Mm^6) (-1 + 3 (u^2)) + 
                                        45 Mm (r^5) (-1 + 
                                        3 (u^2))) + 
                                        (r^3) (10 (Mm^3) q (1 - 
                                        3 (u^2)) + 
                                        15 q (r^3) (1 - 3 (u^2)) + 
                                        4 (Mm^2) r (4 + 5 q - 
                                        15 q (u^2)) + 
                                        Mm (r^2) (-8 + 
                                        45 q (-1 + 
                                        3 (u^2))))) + ((15 ((jj^2) + q) (2 Mm - r) r (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))])/(16 (Mm^2)))) ((1/(r^4)) (2 Mm ((jj^2) (16 (Mm^5) + 
                                        8 (Mm^4) r - 10 (Mm^2) (r^3) + 
                                        15 Mm (r^4) + 15 (r^5)) (1 - 
                                        3 (u^2)) + 
                                        (r^3) (15 q (r^2) (1 - 
                                        3 (u^2)) + 
                                        10 (Mm^2) q (-1 + 
                                        3 (u^2)) + Mm r (8 + 15 q - 
                                        45 q (u^2))))) - 
                                      15 ((jj^2) + q) (2 (Mm^2) - 
                                        (r^2)) (-1 + 3 (u^2)) Log[
                                        (r/(-2 Mm + r))]))))));
dzdrHT[Mm_, jj_, q_, r_, 
          u_] = -(u (Sqrt[1 - (u^2)]) (2048 (jj^2) (Mm^8) (r^5) (-1 + (u^2)) - 3 (-1 + 
                          (u^2)) (2 Mm (16 (jj^2) (Mm^5) + 
                                8 (jj^2) (Mm^4) r - 
                                10 (Mm^2) ((jj^2) + q) (r^3) + 
                                15 Mm ((jj^2) + q) (r^4) + 
                                15 ((jj^2) + q) (r^5)) - 
                          15 ((jj^2) + q) (r^4) (-2 (Mm^2) + 
                                (r^2)) Log[
                              (r/(-2 Mm + r))]) (2 Mm ((jj^2) (8 (Mm^4) (r^2) (1 - 3 (u^2)) + 10 (Mm^3) (r^3) (1 - 3 (u^2)) + 
                                      20 (Mm^2) (r^4) (1 - 3 (u^2)) + 
                                      15 (r^6) (1 - 3 (u^2)) - 
                                      8 (Mm^5) r (-1 + (u^2)) + 
                                      16 (Mm^6) (-1 + 3 (u^2)) + 
                                      45 Mm (r^5) (-1 + 
                                        3 (u^2))) + 
                                (r^3) (10 (Mm^3) q (1 - 3 (u^2)) + 
                                      15 q (r^3) (1 - 3 (u^2)) + 
                                      4 (Mm^2) r (4 + 
                                        q (5 - 15 (u^2))) + 
                                      Mm (r^2) (-8 + 
                                        45 q (-1 + 
                                        3 (u^2))))) + 
                          15 ((jj^2) + 
                                q) (r^5) (((-2 Mm + r))^2) (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))]) + 
                    16 (Mm^2) (-1 + ((2 Mm)/r)) (r^9) ((1/(8 Mm (2 Mm - r) (r^4))) ((jj^2) (-8 (Mm^5) r (-1 + 
                                        u) (1 + u) + 
                                      8 (Mm^4) (r^2) (1 - 3 (u^2)) + 
                                      10 (Mm^3) (r^3) (1 - 3 (u^2)) + 
                                      20 (Mm^2) (r^4) (1 - 3 (u^2)) + 
                                      15 (r^6) (1 - 3 (u^2)) + 
                                      16 (Mm^6) (-1 + 3 (u^2)) + 
                                      45 Mm (r^5) (-1 + 
                                        3 (u^2))) + 
                                (r^3) (10 (Mm^3) q (1 - 3 (u^2)) + 
                                      15 q (r^3) (1 - 3 (u^2)) + 
                                      4 (Mm^2) r (4 + 5 q - 
                                        15 q (u^2)) + 
                                      Mm (r^2) (-8 + 
                                        45 q (-1 + 
                                        3 (u^2))))) + ((15 ((jj^2) + q) (2 Mm - r) r (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))])/(16 (Mm^2)))) ((1/(r^4)) (2 Mm ((jj^2) (16 (Mm^5) + 
                                        8 (Mm^4) r - 10 (Mm^2) (r^3) + 
                                        15 Mm (r^4) + 15 (r^5)) (1 - 
                                        3 (u^2)) + 
                                    (r^3) (15 q (r^2) (1 - 3 (u^2)) + 
                                        10 (Mm^2) q (-1 + 
                                        3 (u^2)) + 
                                        Mm r (8 + 15 q - 
                                        45 q (u^2))))) - 
                          15 ((jj^2) + q) (2 (Mm^2) - 
                                (r^2)) (-1 + 3 (u^2)) Log[
                              (r/(-2 Mm + r))]) - (-1 + 
                          (u^2)) (2 Mm (48 (jj^2) (Mm^6) - 
                                8 (jj^2) (Mm^5) r - 24 (jj^2) (Mm^4) (r^2) - 
                                30 (Mm^3) ((jj^2) + q) (r^3) - 
                                60 (Mm^2) ((jj^2) + q) (r^4) + 
                                135 Mm ((jj^2) + q) (r^5) - 
                                45 ((jj^2) + q) (r^6)) + 
                          45 ((jj^2) + 
                                q) (r^5) (((-2 Mm + r))^2) Log[
                              (r/(-2 Mm + r))]) (2 Mm ((jj^2) (16 (Mm^5) + 8 (Mm^4) r - 10 (Mm^2) (r^3) + 15 Mm (r^4) + 
                                      15 (r^5)) (-1 + 3 (u^2)) + 
                                (r^3) (10 (Mm^2) q (1 - 3 (u^2)) + 
                                      15 q (r^2) (-1 + 3 (u^2)) + 
                                      Mm r (-8 - 15 q + 
                                        45 q (u^2)))) - 
                          15 ((jj^2) + q) (r^4) (-2 (Mm^2) + 
                                (r^2)) (-1 + 3 (u^2)) Log[
                              (r/(-2 Mm + r))])))/(64 (Mm^4) (r^7) √(((1 - ((2 Mm)/r)) (r^2) ((1/(r^4)) (2 Mm ((jj^2) (16 (Mm^5) + 8 (Mm^4) r - 10 (Mm^2) (r^3) + 15 Mm (r^4) + 
                                        15 (r^5)) (1 - 3 (u^2)) + 
                                      (r^3) (15 q (r^2) (1 - 
                                        3 (u^2)) + 
                                        10 (Mm^2) q (-1 + 
                                        3 (u^2)) + 
                                        Mm r (8 + 15 q - 
                                        45 q (u^2))))) - 
                            15 ((jj^2) + q) (2 (Mm^2) - 
                                  (r^2)) (-1 + 3 (u^2)) Log[
                                (r/(-2 Mm + r))]))/((Mm^2) ((1/(8 Mm (2 Mm - r) (r^4))) ((jj^2) (80 (Mm^6) (-1 + 
                                        3 (u^2)) + 
                                        8 (Mm^4) (r^2) (-1 + 
                                        3 (u^2)) + 
                                        10 (Mm^3) (r^3) (-1 + 
                                        3 (u^2)) + 
                                        20 (Mm^2) (r^4) (-1 + 
                                        3 (u^2)) - 
                                        45 Mm (r^5) (-1 + 
                                        3 (u^2)) + 
                                        15 (r^6) (-1 + 3 (u^2)) - 
                                        24 (Mm^5) r (-3 + 
                                        7 (u^2))) + 
                                  (r^3) (10 (Mm^3) q (-1 + 
                                        3 (u^2)) + 
                                        15 q (r^3) (-1 + 3 (u^2)) + 
                                        4 (Mm^2) r (4 - 5 q + 
                                        15 q (u^2)) + 
                                        Mm (r^2) (-8 + 
                                        45 q (1 - 
                                        3 (u^2))))) - ((15 ((jj^2) + q) (2 Mm - r) r (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))])/(16 (Mm^2)))))) √((1/(r^2)) ((1 - 
                            (u^2)) (-4 (jj^2) (Mm^4) (-1 + 
                                  (u^2)) - (1/(16 (Mm^2))) ((-1 + ((2 Mm)/r)) (r^4) ((1/(8 Mm (2 Mm - 
                                        r) (r^4))) ((jj^2) (-8 (Mm^5) r (-1 + u) (1 + u) + 8 (Mm^4) (r^2) (1 - 3 (u^2)) + 
                                        10 (Mm^3) (r^3) (1 - 3 (u^2)) + 
                                        20 (Mm^2) (r^4) (1 - 3 (u^2)) + 
                                        15 (r^6) (1 - 3 (u^2)) + 
                                        16 (Mm^6) (-1 + 3 (u^2)) + 
                                        45 Mm (r^5) (-1 + 
                                        3 (u^2))) + 
                                        (r^3) (10 (Mm^3) q (1 - 
                                        3 (u^2)) + 
                                        15 q (r^3) (1 - 3 (u^2)) + 
                                        4 (Mm^2) r (4 + 5 q - 
                                        15 q (u^2)) + 
                                        Mm (r^2) (-8 + 
                                        45 q (-1 + 
                                        3 (u^2))))) + ((15 ((jj^2) + q) (2 Mm - r) r (-1 + 3 (u^2)) Log[(r/(-2 Mm + r))])/(16 (Mm^2)))) ((1/(r^4)) (2 Mm ((jj^2) (16 (Mm^5) + 
                                        8 (Mm^4) r - 10 (Mm^2) (r^3) + 
                                        15 Mm (r^4) + 15 (r^5)) (1 - 
                                        3 (u^2)) + 
                                        (r^3) (15 q (r^2) (1 - 
                                        3 (u^2)) + 
                                        10 (Mm^2) q (-1 + 
                                        3 (u^2)) + Mm r (8 + 15 q - 
                                        45 q (u^2))))) - 
                                      15 ((jj^2) + q) (2 (Mm^2) - 
                                        (r^2)) (-1 + 3 (u^2)) Log[
                                        (r/(-2 Mm + r))]))))));

AΩ[Mm_, 
          r_] = ((15 ((r^3) - 2 (Mm^3)))/(32 (Mm^3))) Log[
            (r/(r - 2 Mm))];
FΩ1[Mm_, 
          r_] = (48 (Mm^7) - 80 (Mm^6) r + 4 (Mm^5) (r^2) - 
                18 (Mm^4) (r^3) + 40 (Mm^3) (r^4) + 10 (Mm^2) (r^5) + 
                15 Mm (r^6) - 
                15 (r^7)) (((16 (Mm^2) (r - 2 Mm) (r^4)))^(-1)) + 
          AΩ[Mm, r];
FΩ2[Mm_, r_] = 
        5 (6 (Mm^4) - 8 (Mm^3) r - 2 (Mm^2) (r^2) - 3 Mm (r^3) + 
                3 (r^4)) (((16 (Mm^2) (r - 2 Mm) r))^(-1)) - 
          AΩ[Mm, r];
ΩHT[Mm_, jj_, q_, 
          r_] = ((Mm^(1/2))/(r^(3/2))) (1 - 
              jj ((Mm^(3/2))/(r^(3/2))) + (jj^2) FΩ1[Mm, r] + 
              q FΩ2[Mm, r]);
Ar[Mm_, 
          r_] = (((16 (Mm^3) (r - 6 Mm)))^(-1)) (15 r (r - 
                  2 Mm) (2 (Mm^2) + 13 Mm r - 4 (r^2))) Log[
            (r/(r - 2 Mm))];
Fr1[Mm_, 
          r_] = ((6 (Mm^(3/2)) (r + 2 Mm))/((r^(3/2)) (r - 6 Mm)));
Fr2[Mm_, 
        r_] := (((8 (Mm^2) (r^4) (r - 2 Mm) (r - 6 Mm)))^(-1)) (384 (Mm^8) - 720 (Mm^7) r - 112 (Mm^6) (r^2) - 76 (Mm^5) (r^3) - 
              138 (Mm^4) (r^4) - 130 (Mm^3) (r^5) + 635 (Mm^2) (r^6) - 
              375 Mm (r^7) + 60 (r^8)) + 
        Ar[Mm, r]
Fr3[Mm_, 
        r_] := (((8 (Mm^2) r (r - 2 Mm) (r - 6 Mm)))^(-1)) (5 (48 (Mm^5) + 30 (Mm^4) r + 26 (Mm^3) (r^2) - 127 (Mm^2) (r^3) + 
                  75 Mm (r^4) - 12 (r^5))) - 
        Ar[Mm, r]
vHTr[Mm_, jj_, q_, 
        r_] := (((Mm (r - 6 Mm) (r^(-4)) (1 + jj Fr1[Mm, r] - (jj^2) Fr2[Mm, r] - q Fr3[Mm, r])))^(1/2))
rISCOHT[Mm_, jj_, q_] := 
      6 Mm (1 - jj ((((2/3)))^(3/2)) + 
            (jj^2) ((251647/2592) - 240 Log[3/2]) + 
            q (((-9325)/96) + 
                  240 Log[3/2]))
Bz[Mm_, 
        r_] := ((15 (2  r - Mm) (((r - 2 Mm))^2))/(16 (Mm^3))) Log[(r/(r - 2 Mm))]
G1[Mm_, 
        r_] := ((6 (Mm^(3/2)))/(r^(3/2)))
G2[Mm_, 
        r_] := (((8 (Mm^2) (r^4) (r - 2 Mm)))^(-1)) (48 (Mm^7) - 
              244 (Mm^6) r + 28 (Mm^5) (r^2) + 6 (Mm^4) (r^3) - 
              170 (Mm^3) (r^4) + 295 (Mm^2) (r^5) - 165 Mm (r^6) + 
              30 (r^7)) - Bz[Mm, r]
G3[Mm_, 
        r_] := (((8 (Mm^2) r (r - 2 Mm)))^(-1)) (5 (6 (Mm^4) + 34 (Mm^3) r - 59 (Mm^2) (r^2) + 33 Mm (r^3) - 6 (r^4))) + 
        Bz[Mm, r]
vHTz[Mm_, jj_, q_, 
        r_] := (((Mm (r^(-3)) (1 - jj G1[Mm, r] + (jj^2) G2[Mm, r] + q G3[Mm, r])))^(1/2))
A2[Mm_, 
        r_] := -15 ((((r^2) - 2 (Mm^2)))/(16 (Mm^2))) Log[
          (r/(r - 2 Mm))]
H2[Mm_, 
        r_] := (((8 Mm r))^(-1)) 5 (2 (Mm^2) - 3 Mm r - 
              3 (r^2)) - A2[Mm, r]
H1[Mm_, 
        r_] := (((8 Mm (r^4)))^(-1)) (16 (Mm^5) + 8 (Mm^4) r - 
              10 (Mm^2) (r^3) + 15 Mm (r^4) + 15 (r^5)) + 
        A2[Mm, r]
RcircHT[Mm_, jj_, q_, 
          r_] = (Sqrt[(r^2) (1 + (jj^2) H1[Mm, r] + q H2[Mm, r])]);
Be[Mm_, 
        r_] := ((15 r (8 (Mm^2) - 7 Mm r + 2 (r^2)))/(32 (Mm^2) (r - 3 Mm))) Log[(r/(r - 2 Mm))]
Fe1[Mm_, r_] := 
      ((Mm^(5/2))/((Sqrt[r]) (r - 2 Mm) (r - 3 Mm)))
Fe2[Mm_, 
        r_] := (((16 Mm (r^4) (r - 2 Mm) (((r - 3 Mm))^2)))^(-1)) (144 (Mm^8) - 144 (Mm^7) r - 28 (Mm^6) (r^2) - 
              58 (Mm^5) (r^3) - 176 (Mm^4) (r^4) + 685 (Mm^3) (r^5) - 
              610 (Mm^2) (r^6) + 225 Mm (r^7) - 30 (r^8)) + 
        Be[Mm, r]
Fe3[Mm_, 
        r_] := ((5 (r - Mm) (6 (Mm^3) - 20 (Mm^2) r - 21 Mm (r^2) + 6 (r^3)))/(16 Mm r (r - 2 Mm) (r - 3 Mm))) - 
        Be[Mm, r]
EtildaHT[Mm_, jj_, q_, 
        r_] := ((r - 2 Mm)/(Sqrt[r (r - 3 Mm)])) (1 - 
            jj Fe1[Mm, r] + (jj^2) Fe2[Mm, r] + 
            q Fe3[Mm, r])
DEtildaHT[Mm_, jj_, q_, 
          r_] = -ΩHT[Mm, jj, q, 
              r] D[EtildaHT[Mm, jj, q, r], r]/D[ΩHT[Mm, jj, q, r], r];

reading the numerical data and calculating everything in Weyl-Papapetrou coordinates

pathread = "D:\mathematica notebook\NS-model";
SetDirectory[pathread];
\(\(eosFilesTab = 
        Sort[FileNames[], \((StringCount[#1, "M"]\ StringCount[#1, "l"] + 
                  StringCount[#1, "L"]\ StringCount[#1, "p"])\) < \((StringCount[#2, "M"]\ StringCount[#2, "l"] + 
                  StringCount[#2, "L"]\ StringCount[#2, "p"])\) &];\)\)
Dimensions[eosFilesTab]
Table[(StringCount[eosFilesTab[[i]], "M"] StringCount[
                eosFilesTab[[i]], "l"] + 
            StringCount[eosFilesTab[[i]], "L"] StringCount[
                eosFilesTab[[i]], "p"]), {i, 
          1, Dimensions[
              eosFilesTab][[1]]}];
Count[%, 0]
eosFilesTabNew = Table[eosFilesTab[[i]], {i, 1, %}]

<< NumericalMath`ListIntegrate`

overallstarttime = 
        SessionTime[];
<< 
      Statistics`LinearRegression`
Off[
      General::"spell1"]; Off[FindFit::"fmgz"]; 
    Off[FindRoot::"lstol"]; Off[FindRoot::"jsing"]; 
    Off[FindRoot::"cvmit"]; Off[NIntegrate::"slwcon"]; 
    Off[NIntegrate::"ncvb"]; Off[FindRoot::"nlnum"]; 
    Off[InterpolatingFunction::"dmval"]; 
    Clear[M]
rhoSeries = 
        Series[2 Log[((1 - (M/(2 r)))/(1 + (M/(2 r))))] - 
            Log[1 - ((M^2)/(4 (r^2)))], {r, ∞, 
            6}];
rhoSeriesTab = 
        Table[SeriesCoefficient[rhoSeries, i], {i, 0, 
            6}];
TableForm[
        Table[(Subscript[GegeP, n])[z_] = 
            FullSimplify[(Sqrt[(π/2)]) (((((-1))^n) Gamma[n + 2])/((2^(n + 1/2)) n! Gamma[n + 3/2])) (((1 - (z^2)))^(-1/2)) D[(((1 - (z^2)))^(n + 1/2)), {z, n}]], {n, 0, 
            5}]];
assympt = 
        300;
DRiscoTabP = 
        Table[{Null, Null, Null, Null, Null, Null, Null, Null}, {i, 
            1, Dimensions[
                eosFilesTabNew][[1]]}];
DRiscoTabM = 
        Table[{Null, Null, Null, Null, Null, Null, Null, Null}, {i, 
            1, Dimensions[
                eosFilesTabNew][[1]]}];
ModelTypeTab = 
        Table[{Null, Null, Null, Null, Null, Null, Null, Null, Null, Null, 
            Null, Null, Null, Null, Null, Null, Null, Null, Null, Null, Null, 
            Null, Null, Null, Null, Null, Null, Null, Null, Null, Null, 
            Null}, {i, 
            1, Dimensions[eosFilesTabNew][[1]] + 
              1}];
DEtildaPlot = 
        Table[Null, {i, 
            1, Dimensions[
                eosFilesTabNew][[1]]}];
\(ModelTypeTab[\([1]\)] = {"model", "type", \*""d=\!\(\@\(\((\(-a\)\\ b + k)\)\^2 - b\^2\\ M\^2\)\)"", \*""(1/2)(\!\(\(κ\_+\)\)-\!\(\(κ\_-\)\))"", \*""(1/2)(\!\(\(κ\_+\)\)+\!\(\(κ\_-\)\))"", "rotation j", "Manko", "b_-", "M", "J", "M_2", "J_3", "M_4", "h_+", "DRisco_+", "a", "b", "k", "ρ_c", "Ωisco", \*""(Subscript[u, ISCO])"", \*""L/\!\(M\& . \)"", "B0/M^2", \*""(Subscript[σ, tt])"", \*""(Subscript[σ, t])φ"", \*""(Subscript[σ, φ])φ"", \*""(Subscript[σ, zz])"", \*""Axis-(Subscript[σ, tt])"", \*""(Subscript[σ, Δ])E"", \*""(Subscript[σ, Ω])"", \*""(Subscript[σ, Ω])z"", \*""(Subscript[σ, Ω])r""}\)
modIn = 
        1;
starttime = 
        SessionTime[];
SetDirectory[
        pathread];
dataFileName = 
        StringDrop[
          eosFilesTabNew[[modIn]], -4];
Print[dataFileName];
ModelTypeTab[[modIn + 1,1]]=dataFileName;
fileToRead1 = 
        dataFileName <> "-p9.txt";
fileToRead2 = dataFileName <> ".txt";
tableMoments=Import[fileToRead1,"Table"];
tabM1 = 
          Flatten[tableMoments[[Dimensions[tableMoments]]]];
tabM2=Flatten[tableMoments[[4]]];
Print[
        Table[{tabM1[[i]], tabM2[[i]]}, {i, 
              1, Dimensions[tabM1][[1]]}] // MatrixForm];tableData=Import[fileToRead2,"Table"];
tableData[[Dimensions[tableMoments][[1]]+21]];
tableData[[Dimensions[tableMoments][[1]]+23]];
numtable=Table[tableData[[Dimensions[tableMoments][[1]]+22+i]],{i, 
                1, Dimensions[tableData][[1]] - Dimensions[
                      tableMoments][[1]] - 22}];
numericaldata = numtable;
Jgeo = tabM1[[33]];Mgeo = tabM1[[31]];ageo = Jgeo/Mgeo;Rs = tabM1[[4]];M2num = tabM1[[32]];If[Jgeo == N[0], M2num = 0. , M2num = tabM1[[32]]];S3num = tabM1[[34]];hISCOminus = tabM1[[11]];hISCOplus = tabM1[[10]];renum = tabM1[[16]];re = renum;ModelTypeTab[[modIn + 1,9]]=Mgeo;ModelTypeTab[[modIn + 1,10]]=Jgeo;ModelTypeTab[[modIn + 1,11]]=S3num;ModelTypeTab[[modIn + 1,12]]=M2num;ModelTypeTab[[modIn + 1,14]]=hISCOplus;ModelTypeTab[[modIn + 1,19]]=tabM1[[1]];"ρ_c";
Clear[a, M, a2, 
        c];
M = Mgeo;a = ageo;J = a M;M2 = -((a^2) + a2f + a cf) M;S3 = -((a^3) + 2 (a^2) cf + a2f cf + 
                a (2 a2f + (cf^2))) M;If[Jgeo == N[0], 
        twosolitonpar = {{a2f -> 0. , cf -> 3. }}, 
        twosolitonpar = 
          Simplify[
            Solve[{M2 == M2num, S3 == S3num}, {a2f, cf}]]];a = ageo;"km";M = Mgeo;"km";Clear[c, a2];c = (cf //. Flatten[twosolitonpar]);"km";a2 = (a2f //. Flatten[twosolitonpar]);"(km^2)";ModelTypeTab[[modIn + 1,16]]=a;ModelTypeTab[[modIn + 1,17]]=-c;ModelTypeTab[[modIn + 1,18]]=-a2 - a c;ModelTypeTab[[modIn + 1,13]]=((a^4)-(3 (a^2) - 2 a b + (b^2)) k+(k^2)-(1/7) k (M^2)) M//.{k -> -a2 - a c, b -> -c};irad = 601;jpol = 301;kN = irad*jpol;re = renum;xlist = Table[0, {i, kN}];rlist = Table[0, {i, kN}];coslist = Table[0, {i, kN}];rholist = Table[0, {i, kN}];gammalist = Table[0, {i, kN}];alphalist = Table[0, {i, kN}];omegalist = Table[0, {i, kN}];pressurelist = Table[0, {i, kN}];
For[i = 1, i < kN + 1, 
        xlist[[i]] = numericaldata[[i, 1]]; 
        rlist[[i]] = 
          numericaldata[[i, 1]] re/(1 - numericaldata[[i, 1]]); 
        coslist[[i]] = numericaldata[[i, 2]]; 
        rholist[[i]] = numericaldata[[i, 3]]; 
        gammalist[[i]] = numericaldata[[i, 4]]; 
        alphalist[[i]] = numericaldata[[i, 5]]; 
        omegalist[[i]] = numericaldata[[i, 6]]; 
        pressurelist[[i]] = numericaldata[[i, 7]]; 
        i = i + 1];
Clear[
        numericaldata];
rhoOFm = 
      Table[{coslist[[i]], rholist[[i]]}, {i, kN - 1 jpol + 1, 
          kN - 0 jpol}]; 
    omegaOFm = 
      Table[{coslist[[i]], omegalist[[i]]/29.979}, {i, 
          kN - 1 jpol + 1, kN - 0 jpol}];
Nωm = Interpolation[omegaOFm]; 
      Nrhom = Interpolation[rhoOFm];
 (*==========================================================================*) 
Print["=============== Correcting the moments =================="];
Clear[b];
BTab = 
        Table[{1/rlist[[(i - 1)*jpol + j]], 
            coslist[[(i - 1)*jpol + j]], 
            Exp[gammalist[[(i - 1)*jpol + j]]]}, {i, 601 - assympt, 
            601}, {j, 1, 301}];
BTabIntG0 = 
        Table[{BTab[[k, 1, 1]], 
            (4/π) Sum[(1/2) ((BTab[[k, i, 3]] (Subscript[GegeP, 0])[
                              BTab[[k, i, 
                                  2]]] (((1 - BTab[[k, i, 2]]^2))^(1/2)) + BTab[[k, i + 1, 3]] (Subscript[GegeP, 0])[
                              BTab[[k, i + 1, 
                                  2]]] (((1 - BTab[[k, i + 1, 2]]^2))^(1/2))) (BTab[[k, i + 1, 2]] - 
                          BTab[[k, i, 2]])), {i, 1, 300}]}, {k, 
            1, Dimensions[
                BTab][[1]]}];
kmax = 2; 
    vFit = (∑_(k = 0)^kmax) (Subscript[b, k]) (x^k); 
    BCoeffG0 = 
      FindFit[BTabIntG0, vFit, Table[(Subscript[b, k]), {k, 0, kmax}], x, 
        MaxIterations -> 
          100000];
\(Print[\*""Coefficients of (Subscript[G, 0])(μ) for the function B: "", 
        BCoeffG0];\)
\(Print[\*""-((M^2)/4) = "", \(-\(M\^2\/4\)\) /. {M -> 
              Mgeo}, \*"" -> (Subscript[b, 2])""];\)
\(Print["b = ", \(b\_2/M\^2 /. {M -> Mgeo}\) /. 
          BCoeffG0, " , b based on Laarakkers & Poisson should be -1/4"];\)
ModelTypeTab[[modIn + 1,23]]=(Subscript[b, 2])/(M^2)/.{M -> Mgeo}/.BCoeffG0;ModelTypeTab[[modIn + 1,12]]=S3num - (12/5) (((M^2)/4) + (Subscript[b, 2])) Jgeo /. {M -> 
                  Mgeo} /. BCoeffG0;If[Jgeo == N[0], S3New = 0, 
          S3New = S3num - (12/5) (((M^2)/4) + (Subscript[b, 2])) Jgeo /. {M -> 
                    Mgeo} /. BCoeffG0];\(Print[\*""\!\(S\_3\) corrected= "", S3New]\);ModelTypeTab[[modIn + 1,11]]=M2num - (4/3) (((M^2)/4) + (Subscript[b, 2])) M /. {M -> Mgeo} /. 
            BCoeffG0;If[Jgeo == N[0], M2New = 0, 
          M2New = M2num - (4/3) (((M^2)/4) + (Subscript[b, 2])) M /. {M -> 
                    Mgeo} /. BCoeffG0];\(Print[\*""\!\(M\_2\) corrected= "", M2New]\);  
If[Jgeo == N[0], 
          twosolitonpar2 = {{a2f -> 0. , cf -> 3. }}, 
          twosolitonpar2 = 
            Simplify[
              Solve[{M2 == 
                        M2num - (4/3) (((M^2)/4) + (Subscript[b, 2])) M /. {M -> 
                          Mgeo} /. 
                    BCoeffG0, S3 == 
                        S3num - (12/5) (((M^2)/4) + (Subscript[b, 2])) Jgeo /. {M -> 
                          Mgeo} /. BCoeffG0}, {a2f, cf}]]];a = ageo;"km";M = Mgeo;"km";cNew = (cf //. Flatten[twosolitonpar2]);"km";a2New = (a2f //. Flatten[twosolitonpar2]);"(km^2)";ModelTypeTab[[modIn + 1,16]]=a;ModelTypeTab[[modIn + 1,17]]=-cNew;ModelTypeTab[[modIn + 1,18]]=-a2New - a cNew;ModelTypeTab[[modIn + 1,13]]=((a^4)-(3 (a^2) - 2 a b + (b^2)) k+(k^2)-(1/7) k (M^2)) M//.{k -> -a2New - a cNew, b -> -cNew};Print["a= ", a];Print["M= ", M];Print["(Subscript[M, 4]) two-soliton= ",ModelTypeTab[[modIn + 1,13]]];
  
Print["============= Introducing the numerical \
solution on the equatorial plane ===================="];
Nflist = Table[{0, 0}, {i, irad}]; 
    Ngtt = Table[{0, 0}, {i, irad}]; 
    NDflist = Table[{0, 0}, {i, irad - 303}]; 
    NDDflist = Table[{0, 0}, {i, irad - 303}]; 
    Nωlist = Table[{0, 0}, {i, irad}]; 
    Ngtφ = Table[{0, 0}, {i, irad}]; 
    NDωlist = Table[{0, 0}, {i, irad - 303}]; 
    NDDωlist = Table[{0, 0}, {i, irad - 303}]; 
    pressurelistρ = Table[{0, 0}, {i, irad}]; 
    glistρ = Table[{0, 0}, {i, irad}]; 
    rhoOFr = Table[{0, 0}, {i, irad - 471}]; 
    omegaOFr = 
      Table[{0, 0}, {i, irad - 471}];
j = 
        1;
\(For[i = 1, i < kN, 
        If[j > 471, omegaOFr[\([j - 471, 1]\)] = 1/rlist[\([i]\)]; 
          omegaOFr[\([j - 471, 2]\)] = omegalist[\([i]\)]/29.979; 
          rhoOFr[\([j - 471, 1]\)] = 1/rlist[\([i]\)]; 
          rhoOFr[\([j - 471, 2]\)] = rholist[\([i]\)]]; 
        Nflist[\([j, 
              2]\)] = \((Exp[
                rholist[\([i]\)] + 
                  gammalist[\([i]\)]] - \((omegalist[\([i]\)]/29.979)\)\^2\ rlist[\([i]\)]\^2\ \((1 - coslist[\([i]\)]\^2)\)\ Exp[
                  gammalist[\([i]\)] - rholist[\([i]\)]])\); 
        Nflist[\([j, 1]\)] = 
          rlist[\([i]\)]\ Sin[ArcCos[coslist[\([i]\)]]]\ Exp[
              gammalist[\([i]\)]]; 
        If[j > 302 && j < irad, 
          NDflist[\([j - 302, 
                2]\)] = \((re/\((rlist[\([i]\)] + re)\)\^2)\)\ \(\((Exp[
                        gammalist[\([i]\)]] + 
                      rlist[\([i]\)]\ \((re/\((rlist[\([i]\)] + re)\)\^2)\)\ \((Exp[gammalist[\([i + jpol]\)]] - 
                              Exp[gammalist[\([i - 
                                      jpol]\)]])\)/\((xlist[\([i + jpol]\)] - 
                              xlist[\([i - jpol]\)])\))\)\^\(-1\)\) \((\((Exp[
                          rholist[\([i + jpol]\)] + 
                            gammalist[\([i + 
                                  jpol]\)]] - \((omegalist[\([i + jpol]\)]/29.979)\)\^2\ rlist[\([i + jpol]\)]\^2\ \((1 - 
                              coslist[\([i + jpol]\)]\^2)\)\ Exp[
                            gammalist[\([i + jpol]\)] - 
                              rholist[\([i + jpol]\)]])\) - 
                    Nflist[\([j - 1, 2]\)])\)/\((xlist[\([i + jpol]\)] - 
                    xlist[\([i - jpol]\)])\); 
          NDflist[\([j - 302, 1]\)] = 
            Nflist[\([j, 
                1]\)]; \*""ddrddρ=-(2 (re/(((rlist[[i]] + re))^2)) (Exp[gammalist[[i+jpol]]]-Exp[gammalist[[i-jpol]]])/(xlist[[i+jpol]]-xlist[[i-jpol]])+rlist[[i]] ((((re/(((rlist[[i]] + re))^2)))^2) (Exp[gammalist[[i+jpol]]]+Exp[gammalist[[i-jpol]]]-2\ Exp[gammalist[[i]]])/((((xlist[[i + jpol]] - xlist[[i - jpol]])/2))^2)+(Exp[gammalist[[i+jpol]]]-Exp[gammalist[[i-jpol]]])/(xlist[[i+jpol]]-xlist[[i-jpol]]) (-2 re/(((rlist[[i]] + re))^3))))/\!\(\((Exp[gammalist[\([i]\)]] + rlist[\([i]\)]\\\\ \((re/\((rlist[\([i]\)] + re)\)\^2)\)\\\\ \((Exp[gammalist[\([i + jpol]\)]] - Exp[gammalist[\([i - jpol]\)]])\)/\((xlist[\([i + jpol]\)] - xlist[\([i - jpol]\)])\))\)\^3\)""; "NDDflist[[j-302,2]]="; 
          NDDflist[\([j - 302, 1]\)] = Nflist[\([j, 1]\)], Null]; 
        Ngtt[\([j, 2]\)] = \(-Nflist[\([j, 2]\)]\); 
        Ngtt[\([j, 1]\)] = Nflist[\([j, 1]\)]; 
        Nωlist[\([j, 
              2]\)] = \(-\((\((omegalist[\([i]\)]/
                      29.979)\)\ rlist[\([i]\)]\^2\ \((1 - 
                      coslist[\([i]\)]\^2)\)\ Exp[
                    gammalist[\([i]\)] - rholist[\([i]\)]])\)\)/
            Nflist[\([j, 2]\)]; 
        Nωlist[\([j, 1]\)] = Nflist[\([j, 1]\)]; 
        If[j > 302 && j < irad, 
          NDωlist[\([j - 302, 
                2]\)] = \((re/\((rlist[\([i]\)] + re)\)\^2)\)\ \((Exp[gammalist[\([i]\)]] + rlist[\([i]\)]\ \((re/\((rlist[\([i]\)] + re)\)\^2)\)\ \((Exp[gammalist[\([i + jpol]\)]] - Exp[gammalist[\([i - jpol]\)]])\)/\((xlist[\([i + jpol]\)] - xlist[\([i - jpol]\)])\))\)\^\(-1\)\ \((\(-\((\((omegalist[\([i + jpol]\)]/
                                29.979)\)\ rlist[\([i + jpol]\)]\^2\ \((1 - 
                                coslist[\([i + jpol]\)]\^2)\)\ Exp[
                              gammalist[\([i + jpol]\)] - 
                                rholist[\([i + jpol]\)]])\)\)/\((Exp[
                            rholist[\([i + jpol]\)] + 
                              gammalist[\([i + 
                                    jpol]\)]] - \((omegalist[\([i + jpol]\)]/29.979)\)\^2\ rlist[\([i + jpol]\)]\^2\ \((1 - 
                                coslist[\([i + jpol]\)]\^2)\)\ Exp[
                              gammalist[\([i + jpol]\)] - 
                                rholist[\([i + jpol]\)]])\) - 
                    Nωlist[\([j - 1, 
                        2]\)])\)/\((xlist[\([i + jpol]\)] - 
                    xlist[\([i - jpol]\)])\); 
          NDωlist[\([j - 302, 1]\)] = Nωlist[\([j, 1]\)], 
          Null]; Ngtφ[\([j, 2]\)] = 
          Nflist[\([j, 2]\)]\ Nωlist[\([j, 2]\)]; 
        Ngtφ[\([j, 1]\)] = Nflist[\([j, 1]\)]; 
        pressurelistρ[\([j, 2]\)] = pressurelist[\([i]\)]; 
        pressurelistρ[\([j, 1]\)] = Nflist[\([j, 1]\)]; 
        glistρ[\([j, 2]\)] = gammalist[\([i]\)]; 
        glistρ[\([j, 1]\)] = Nflist[\([j, 1]\)]; j = j + 1; 
        i = i + jpol]; Clear[j]; re = renum; SDIV = 601; MDIV = 301; 
      Ngzz = Table[{0, 0}, {i, SDIV - 303}]; j = 1; 
      For[i = \((302\ MDIV + 1)\), i < \((SDIV - 1)\)\ MDIV + 1, 
        Ngzz[\([j, 1]\)] = 
          rlist[\([i]\)]\ Sin[ArcCos[coslist[\([i]\)]]]\ Exp[
              gammalist[\([i]\)]]; 
        Ngzz[\([j, 
              2]\)] = \((Exp[gammalist[\([i]\)]] + rlist[\([i]\)]\ \((1 - coslist[\([i]\)]\^2)\)\ \((re\/\((re + rlist[\([i]\)])\)\^2\ \((Exp[gammalist[\([i + MDIV]\)]] - Exp[gammalist[\([i - MDIV]\)]])\)/\((xlist[\([i + MDIV]\)] - xlist[\([i - MDIV]\)])\))\))\)\^\(-2\)\ Exp[2\ alphalist[\([i]\)]]; j = j + 1; 
        i = i + jpol]; InNgzz = Interpolation[Ngzz];\)
 (*==========================================================================*) 
Print["================ introducing the numerical solution on axis z \
================================"];
re = renum; SDIV = 601; MDIV = 301; Le = 29.979; 
      surfaceTabR = 
        Table[{((re*numtable[[i, 1]])/(1 - numtable[[i, 1]])), 
            numtable[[i, 2]], Exp[numtable[[i, 4]]], 
            numtable[[i, 7]]}, {i, 1, SDIV*MDIV}]; 
      surfaceTabX = 
        Table[{numtable[[i, 1]], numtable[[i, 2]], 
            Exp[numtable[[i, 4]]], numtable[[i, 7]]}, {i, 1, 
            SDIV*MDIV}]; rOfTheta = Table[Null, {i, 1, MDIV}]; 
      rOfTheta[[1]] = 301; 
      For[i = 2, i < 302, j = 0; 
        While[surfaceTabR[[(((rOfTheta[[i - 1]] - j) MDIV) + 
                    i), 4]] == 0, j = j + 1]; 
        rOfTheta[[i]] = rOfTheta[[i - 1]] - j + 1; i = i + 1]; 
      rOfTheta; surfaceTabR[[(rOfTheta[[1]] - 1)*MDIV + 1, 4]]; 
      surfaceTabR[[rOfTheta[[1]]*MDIV + 1, 4]];
 (*==========================================================================*) 
\(listdBdr = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 
                2]\)], \((\((surfaceTabX[\([\((\((rOfTheta[\([1]\)] + 2)\)*
                                MDIV)\) + \((i)\), 3]\)] - 
                      surfaceTabX[\([\((\((\((rOfTheta[\([1]\)])\)*
                                  MDIV)\) + \((i)\))\), 
                          3]\)])\)/\((surfaceTabX[\([\((\((\((rOfTheta[\([1]\)] + 2)\)*MDIV)\) + \((i)\))\), 1]\)] - 
                      surfaceTabX[\([\((\((\((\((rOfTheta[\([1]\)])\)*
                                    MDIV)\))\) + \((i)\))\), 
                          1]\)])\))\)\ \((re/\((re + surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 1]\)])\)\^2)\)}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdBdμ = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            If[i == 1, 0, 
              If[i == MDIV, 
                0, \((surfaceTabR[\([\((\((rOfTheta[\([1]\)] + 1)\)*
                                MDIV)\) + \((i + 1)\), 3]\)] - 
                      surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                  MDIV)\) + \((i - 1)\))\), 
                          3]\)])\)/\((surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\) + \((i + 1)\))\), 2]\)] - 
                      surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                    MDIV)\))\) + \((i - 1)\))\), 
                          2]\)])\)]]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdBdω = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 
                2]\)], √\((1 - 
                      surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 2]\)]\^2)\)\ listdBdr[\([i, 3]\)] - 
              surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\))\) + \((i)\))\), 
                    2]\)]\ \((√\((1 - 
                          surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 2]\)]\^2)\))\)/
                  surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                MDIV)\))\) + \((i)\))\), 
                      1]\)]\ listdBdμ[\([i, 3]\)]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdBdζ = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\))\) + \((i)\))\), 2]\)]\ listdBdr[\([i, 
                    3]\)] + \((1 - 
                      surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 2]\)]\^2)\)/
                  surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                MDIV)\))\) + \((i)\))\), 
                      1]\)]\ listdBdμ[\([i, 3]\)]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdρdω = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                            MDIV)\))\) + \((i)\))\), 3]\)]\  + 
              surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\))\) + \((i)\))\), 
                    1]\)]\ √\((1 - 
                      
                      surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 2]\)]\^2)\)\ listdBdω[\([i, 3]\)]}, {i, 
            1, MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; listdzdω = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 
                2]\)], \(-surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\))\) + \((i)\))\), 1]\)]\)\ √\((1 - 
                    surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 2]\)]\^2)\)\ listdBdζ[\([i, 3]\)]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdrdθ = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            If[i == 1, 0, 
              If[i == MDIV, 
                0, \((surfaceTabR[\([\((\((rOfTheta[\([1]\)] + 1)\)*
                                MDIV)\) + \((i + 1)\), 1]\)] - 
                      surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                  MDIV)\) + \((i - 1)\))\), 
                          1]\)])\)/\((ArcCos[
                        surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                    MDIV)\) + \((i + 1)\))\), 2]\)]] - 
                      ArcCos[surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 
                                        1)\)*MDIV)\))\) + \((i - 1)\))\), 
                            2]\)]])\)]]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdρdθ = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            If[i == 1, 
              0, \((surfaceTabR[\([\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\) + \((i)\), 
                        1]\)]\ surfaceTabR[\([\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\) + \((i)\), 
                        2]\)]\ surfaceTabR[\([\((\((rOfTheta[\([1]\)] + 1)\)*
                              MDIV)\) + \((i)\), 3]\)] - 
                  surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*
                                MDIV)\) + \((i)\))\), 1]\)]\ \((1 - 
                        
                        surfaceTabR[\([\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\) + \((i)\))\), 2]\)]\^2)\)\ listdBdμ[\([i, 3]\)])\)]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[i]]+1)*MDIV))+(i))"; 
      listdzdθ = 
        Table[{surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 1]\)], 
            surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*
                          MDIV)\))\) + \((i)\))\), 2]\)], 
            If[i == MDIV, 
              0, \((listdρdω[\([i, 3]\)]\^2 + 
                      listdzdω[\([i, 3]\)]\^2)\)\ \((listdrdθ[\([i, 3]\)]\^2 + 
                      surfaceTabR[\([\((\((\((\((rOfTheta[\([1]\)] + 1)\)*MDIV)\))\) + \((i)\))\), 1]\)]\^2)\) - 
                listdρdθ[\([i, 3]\)]\^2]}, {i, 1, 
            MDIV}]; "orismos sto shmeio ((((rOfTheta[[1]]+1)*MDIV))+(i))";\)
 (*==========================================================================*) 
zSurfaceList = 
        Table[If[i == 1, 
            0, -(1/2) (ArcCos[
                    surfaceTabR[[((((rOfTheta[[1]] + 1)*
                                  MDIV)) + (i)), 2]]] - 
                  ArcCos[surfaceTabR[[((((rOfTheta[[1]] + 1)*
                                  MDIV)) + (i - 1)), 
                        2]]]) (listdzdθ[[i, 3]]^(1/2) + 
                  listdzdθ[[i - 1, 3]]^(1/2))], {i, 1, 
            MDIV}]; zSurface = 
        Table[{surfaceTabR[[((((rOfTheta[[1]] + 1)*
                          MDIV)) + (i)), 1]], 
            surfaceTabR[[((((rOfTheta[[1]] + 1)*
                          MDIV)) + (i)), 2]], 
            surfaceTabR[[((((rOfTheta[[1]] + 1)*
                            MDIV)) + (i)), 
                  1]] √(1 - 
                    surfaceTabR[[((((rOfTheta[[1]] + 1)*MDIV)) + (i)), 2]]^2) surfaceTabR[[((((rOfTheta[[1]] + 1)*MDIV)) + (i)), 3]], 
            If[i == 1, 0, Sum[zSurfaceList[[k]], {k, 1, i}]]}, {i, 
            1, MDIV}]; 
      zOnAxis = Table[Null, {i, rOfTheta[[MDIV]] + 2, SDIV}]; 
      k = rOfTheta[[1]] + 2 - (rOfTheta[[MDIV]] + 1); 
      For[i = rOfTheta[[1]] + 2, i < SDIV + 1, 
        zOnAxis[[k]] = 
          If[k == 
              rOfTheta[[1]] + 2 - (rOfTheta[[MDIV]] + 1), 
            zSurface[[MDIV, 4]], 
            zOnAxis[[k - 1]] + 
              (1/2) (surfaceTabR[[i*MDIV, 1]] - 
                    surfaceTabR[[(i - 1)*MDIV, 
                        1]]) (surfaceTabR[[i*MDIV, 3]] + 
                    surfaceTabR[[(i - 1)*MDIV, 3]])]; i = i + 1; 
        k = k + 1]; Clear[k]; 
      k = rOfTheta[[1]] + 2 - (rOfTheta[[MDIV]] + 1); 
      For[i = rOfTheta[[1]] + 2, i > rOfTheta[[MDIV]] + 1, 
        zOnAxis[[k]] = 
          If[i == rOfTheta[[1]] + 2, zSurface[[MDIV, 4]], 
            zOnAxis[[k + 1]] - 
              (1/2) (surfaceTabR[[(i + 1)*MDIV, 1]] - 
                    surfaceTabR[[i*MDIV, 
                        1]]) (surfaceTabR[[(i + 1)*MDIV, 3]] + 
                    surfaceTabR[[i*MDIV, 3]])]; i = i - 1; k = k - 1]; 
      Dimensions[zOnAxis]; Clear[k];
 (*==========================================================================*) 
Nflistz = 
      Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}]; 
    Ngttz = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}]; 
    Nωlistz = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}]; 
    Ngtφz = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}]; 
    j = 1; For[i = (rOfTheta[[MDIV]] + 2)*MDIV, i < kN + 1, 
      Nflistz[[j, 
            2]] = (Exp[
              rholist[[i]] + 
                gammalist[[i]]] - (((omegalist[[i]]/29.979))^2) rlist[[i]]^2 (1 - coslist[[i]]^2) Exp[
                gammalist[[i]] - rholist[[i]]]); 
      Nflistz[[j, 1]] = zOnAxis[[j]]; 
      Ngttz[[j, 2]] = -Nflistz[[j, 2]]; 
      Ngttz[[j, 1]] = Nflistz[[j, 1]]; 
      Nωlistz[[j, 
            2]] = -((omegalist[[i]]/
                    29.979) rlist[[i]]^2 (1 - 
                    coslist[[i]]^2) Exp[
                  gammalist[[i]] - rholist[[i]]])/
          Nflistz[[j, 2]]; 
      Nωlistz[[j, 1]] = Nflistz[[j, 1]]; 
      Ngtφz[[j, 2]] = 
        Nflistz[[j, 2]] Nωlistz[[j, 2]]; 
      Ngtφz[[j, 1]] = Nflistz[[j, 1]]; j = j + 1; 
      i = i + jpol]; Clear[j]; Rsp = zOnAxis[[1]]; 
    NRs = Nflist[[303, 1]];
Clear[j];
 (*==========================================================================*) 
Nωρ = 
        Interpolation[
          Table[{Nωlist[[i, 1]], {Nωlist[[i, 2]], 
                NDωlist[[i - 302, 2]]}}, {i, 303, 
              SDIV - 1}]];
Nfρ = 
        Interpolation[
          Table[{Nflist[[i, 1]], {Nflist[[i, 2]], 
                NDflist[[i - 302, 2]]}}, {i, 303, 
              SDIV - 1}]];
Nωρ2 = 
        Interpolation[
          Table[{Nωlist[[i, 1]], Nωlist[[i, 2]]}, {i, 
              301, SDIV}]];
Nfρ2 = 
        Interpolation[
          Table[{Nflist[[i, 1]], Nflist[[i, 2]]}, {i, 301, 
              SDIV}]];
Nωρ3 = 
        Interpolation[
          Table[{Nωlist[[i, 1]], Nωlist[[i, 2]]}, {i, 
              1, SDIV}]];
Nfρ3 = 
        Interpolation[
          Table[{Nflist[[i, 1]], Nflist[[i, 2]]}, {i, 1, 
              SDIV}]];
Nωr = 
        Interpolation[omegaOFr];
Nrhor = 
        Interpolation[
          rhoOFr];
Print["======== Weyl \
coordinates of the numerical ISCOs and the eq. surface ========="];
rhoISCO2 = 
      FindRoot[(Sqrt[Nfρ2[ρ]^(-1) (ρ^2) - Nfρ2[ρ] Nωρ2[ρ]^2]) == (Rs + hISCOminus), {ρ, 
          NRs}]
rhoISCO1 = 
      FindRoot[(Sqrt[Nfρ2[ρ]^(-1) (ρ^2) - Nfρ2[ρ] Nωρ2[ρ]^2]) == (Rs + hISCOplus), {ρ, 
          NRs}]
rhoSurface = 
        FindRoot[(Sqrt[Nfρ2[ρ]^(-1) (ρ^2) - Nfρ2[ρ] Nωρ[ρ]^2]) == Rs, {ρ, NRs}]
 (*==========================================================================*) 
Print["============== computing the frequencies ================"];
ΩN[ρ_] = (-D[
                  Nfρ[ρ] Nωρ[ρ], ρ] + √(D[Nfρ[ρ] Nωρ[ρ], ρ]^2 - 
                    D[-Nfρ[ρ], ρ] D[
                        Nfρ[ρ]^(-1) (ρ^2) - 
                          Nfρ[ρ] Nωρ[ρ]^2, ρ]))/
          D[Nfρ[ρ]^(-1) (ρ^2) - 
              Nfρ[ρ] Nωρ[ρ]^2, ρ];
\(Off[Power::"infy"]; Off[∞::"indet"]; 
      NgttTab = 
        Table[\(-Exp[gammalist[\([i]\)] + rholist[\([i]\)]]\) + 
            Exp[gammalist[\([i]\)] - 
                  rholist[\([i]\)]]\ \((\(xlist[\([i]\)]\ re\)\/\(1 - xlist[\([i]\)]\))\)\^2\ \((1 - 
                  coslist[\([i]\)]\^2)\) \((omegalist[\([i]\)]/29.979)\)\^2, {i, 1, MDIV*SDIV}]; 
      NgφφTab = 
        Table[Exp[
              gammalist[\([i]\)] - 
                rholist[\([i]\)]]\ \((\(xlist[\([i]\)]\ re\)\/\(1 - xlist[\([i]\)]\))\)\^2\ \((1 - coslist[\([i]\)]\^2)\), {i, 1, MDIV*SDIV}]; 
      NgtφTab = 
        Table[\(-Exp[
                gammalist[\([i]\)] - 
                  rholist[\([i]\)]]\)\ \((\(xlist[\([i]\)]\ re\)\/\(1 - xlist[\([i]\)]\))\)\^2\ \((1 - coslist[\([i]\)]\^2)\) \((omegalist[\([i]\)]/
                29.979)\), {i, 1, MDIV*SDIV}]; 
      NgxxTab = 
        Table[\((Exp[2\ alphalist[\([i]\)]]\ \((re\/\((1 - xlist[\([i]\)])\)\^2)\)\^2\ )\)\^\(-1\), {i, 1, MDIV*SDIV}]; 
      NgμμTab = 
        Table[\((Exp[2\ alphalist[\([i]\)]]\ \((\(xlist[\([i]\)]\ re\)\/\((1 - xlist[\([i]\)])\))\)\^2)\)\^\(-1\)\ \((1 - coslist[\([i]\)]\^2)\), {i, 1, 
            MDIV*SDIV}]; 
      NΩTab = 
        Table[{rlist[\([i*MDIV + 1]\)]\ Sin[
                ArcCos[coslist[\([i*MDIV + 1]\)]]]\ Exp[
                gammalist[\([i*MDIV + 
                      1]\)]], \((\(-\((\((NgtφTab[\([\((i + 1)\)*
                                    MDIV + 1]\)] - 
                            NgtφTab[\([\((i - 1)\)*MDIV + 
                                  1]\)])\)/\((xlist[\([\((i + 1)\)*MDIV + 
                                  1]\)] - 
                            xlist[\([\((i - 1)\)*MDIV + 1]\)])\))\)\) + 
                  Sqrt[\((\((NgtφTab[\([\((i + 1)\)*MDIV + 1]\)] - NgtφTab[\([\((i - 1)\)*MDIV + 1]\)])\)/\((xlist[\([\((i + 1)\)*MDIV + 1]\)] - xlist[\([\((i - 1)\)*MDIV + 1]\)])\))\)\^2 - \((\((NgttTab[\([\((i + 1)\)*MDIV + 1]\)] - NgttTab[\([\((i - 1)\)*MDIV + 
                                      1]\)])\)/\((xlist[\([\((i + 1)\)*MDIV + 
                                      1]\)] - 
                                xlist[\([\((i - 1)\)*MDIV + 
                                      1]\)])\))\)\ \((\((NgφφTab[\([\((i + 1)\)*MDIV + 1]\)] - 
                                NgφφTab[\([\((i - 1)\)*
                                        MDIV + 
                                      1]\)])\)/\((xlist[\([\((i + 1)\)*MDIV + 
                                      1]\)] - 
                                xlist[\([\((i - 1)\)*MDIV + 
                                      1]\)])\))\)])\)/\((\((NgφφTab[\([\((i + 1)\)*MDIV + 1]\)] - 
                      NgφφTab[\([\((i - 1)\)*MDIV + 
                            1]\)])\)/\((xlist[\([\((i + 1)\)*MDIV + 1]\)] - 
                      xlist[\([\((i - 1)\)*MDIV + 1]\)])\))\)}, {i, 301, 
            SDIV - 3}];\)
(*---------------------------------calculating \Delta\tilde{E}---------------------------------*)
NΩTab2 = 
        Table[{xlist[[i*MDIV + 1]], 
            rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 
                      1]]], (-((NgtφTab[[(i + 1)*
                                    MDIV + 1]] - 
                            NgtφTab[[(i - 1)*MDIV + 
                                  1]])/(xlist[[(i + 1)*MDIV + 
                                  1]] - 
                            xlist[[(i - 1)*MDIV + 1]])) + 
                  Sqrt[((((NgtφTab[[(i + 1)*MDIV + 1]] - NgtφTab[[(i - 1)*MDIV + 1]])/(xlist[[(i + 1)*MDIV + 1]] - xlist[[(i - 1)*MDIV + 1]])))^2) - ((NgttTab[[(i + 1)*MDIV + 1]] - NgttTab[[(i - 1)*MDIV + 
                                      1]])/(xlist[[(i + 1)*MDIV + 
                                      1]] - 
                                xlist[[(i - 1)*MDIV + 
                                      1]])) ((NgφφTab[[(i + 1)*MDIV + 1]] - 
                                NgφφTab[[(i - 1)*
                                        MDIV + 
                                      1]])/(xlist[[(i + 1)*MDIV + 
                                      1]] - 
                                xlist[[(i - 1)*MDIV + 
                                      1]]))])/((NgφφTab[[(i + 1)*MDIV + 1]] - 
                      NgφφTab[[(i - 1)*MDIV + 
                            1]])/(xlist[[(i + 1)*MDIV + 1]] - 
                      xlist[[(i - 1)*MDIV + 1]]))}, {i, 301, 
            SDIV - 3}]; 
      NEtildaTab = 
        Table[{xlist[[i*MDIV + 1]], 
            rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 
                      1]]], (-NgttTab[[(i)*MDIV + 1]] - 
                  NgtφTab[[(i)*MDIV + 
                          1]] (NΩTab2[[i - 300, 
                          3]]))/(((-NgttTab[[(i)*MDIV + 1]] - 2 NgtφTab[[(i)*MDIV + 1]] (NΩTab2[[i - 300, 3]]) - NgφφTab[[(i - 1)*MDIV + 1]] (((NΩTab2[[i - 300, 3]]))^2)))^(1/2))}, {i, 301, 
            SDIV - 3}]; 
      NDEtildaTab = 
        Table[{(NEtildaTab[[i, 2]] + NEtildaTab[[i + 1, 2]])/
              2, -((NEtildaTab[[i + 1, 3]] - 
                        NEtildaTab[[i, 
                            3]])/(NΩTab2[[i + 1, 
                            3]] - 
                        NΩTab2[[i, 
                            3]])) (((NΩTab2[[i + 1, 3]] + NΩTab2[[i, 3]])/2))}, {i, 
            2, Dimensions[NEtildaTab][[1]] - 1}];
  
 (*======================= calculating the frequencies ========= = =====================*) 
NΩxTab = 
        Table[{rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 1]]], 
            NΩTab[[i - 300, 
                  2]] - (((((-NgxxTab[[i*MDIV + 1]])/2) ((((NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (NgφφTab[[(i + 1)*MDIV + 1]]/NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgφφTab[[(i)*MDIV + 1]]/NΩTab[[(i) - 300, 1]]^2 + NgφφTab[[(i - 1)*MDIV + 1]]/NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2) - 2 (NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]) (NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]) (NgtφTab[[(i + 1)*MDIV + 1]]/NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgtφTab[[(i)*MDIV + 1]]/NΩTab[[(i) - 300, 1]]^2 + NgtφTab[[(i - 1)*MDIV + 1]]/NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2) + (((NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (NgttTab[[(i + 1)*MDIV + 1]]/NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgttTab[[(i)*MDIV + 1]]/NΩTab[[(i) - 300, 1]]^2 + NgttTab[[(i - 1)*MDIV + 1]]/NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2))))^(1/2))}, {i, 302, SDIV - 4}]; 
      NKx2Tab = 
        Table[{rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 
                      1]]], (((-NgxxTab[[i*MDIV + 1]])/2) ((((NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (NgφφTab[[(i + 1)*MDIV + 1]]/
                              NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgφφTab[[(i)*MDIV + 1]]/
                                NΩTab[[(i) - 300, 1]]^2 + NgφφTab[[(i - 1)*MDIV + 1]]/
                              NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2) - 2 (NgttTab[[i*MDIV + 1]] + 
                          NgtφTab[[i*MDIV + 
                                  1]] NΩTab[[i - 300, 
                                2]]) (NgtφTab[[i*MDIV + 
                                1]] + 
                          NgφφTab[[i*MDIV + 
                                  1]] NΩTab[[i - 300, 
                                2]]) (NgtφTab[[(i + 1)*
                                      MDIV + 1]]/
                              NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgtφTab[[(i)*MDIV + 1]]/
                                NΩTab[[(i) - 300, 1]]^2 + NgtφTab[[(i - 1)*MDIV + 1]]/
                              NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2) + (((NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (NgttTab[[(i + 1)*MDIV + 1]]/
                              NΩTab[[(i + 1) - 300, 1]]^2 - 2 NgttTab[[(i)*MDIV + 1]]/
                                NΩTab[[(i) - 300, 1]]^2 + NgttTab[[(i - 1)*MDIV + 1]]/
                              NΩTab[[(i - 1) - 300, 1]]^2)/(((xlist[[(i + 1)*MDIV + 1]] - xlist[[(i)*MDIV + 1]]))^2)))}, {i, 302, SDIV - 4}]; 
      NΩμTab = 
        Table[{rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 1]]], 
            NΩTab[[i - 300, 
                  2]] - (((((-NgμμTab[[i*MDIV + 1]])/2) ((((NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (2 NgφφTab[[(i)*MDIV + 2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgφφTab[[(i)*MDIV + 1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2) - 2 (NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]) (NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]) (2 NgtφTab[[(i)*MDIV + 2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgtφTab[[(i)*MDIV + 1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2) + (((NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (2 NgttTab[[(i)*MDIV + 2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgttTab[[(i)*MDIV + 1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2))))^(1/2))}, {i, 301, SDIV - 3}]; 
      NKμ2Tab = 
        Table[{rlist[[i*MDIV + 1]] Sin[
                ArcCos[coslist[[i*MDIV + 1]]]] Exp[
                gammalist[[i*MDIV + 
                      1]]], (((-NgμμTab[[i*MDIV + 1]])/2) ((((NgttTab[[i*MDIV + 1]] + NgtφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (2 NgφφTab[[(i)*MDIV + 
                                      2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgφφTab[[(i)*MDIV + 
                                      1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2) - 2 (NgttTab[[i*MDIV + 1]] + 
                          NgtφTab[[i*MDIV + 
                                  1]] NΩTab[[i - 300, 
                                2]]) (NgtφTab[[i*MDIV + 
                                1]] + 
                          NgφφTab[[i*MDIV + 
                                  1]] NΩTab[[i - 300, 
                                2]]) (2 NgtφTab[[(i)*
                                        MDIV + 
                                      2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgtφTab[[(i)*MDIV + 
                                      1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2) + (((NgtφTab[[i*MDIV + 1]] + NgφφTab[[i*MDIV + 1]] NΩTab[[i - 300, 2]]))^2) (2 NgttTab[[(i)*MDIV + 
                                      2]]/(((rlist[[i*MDIV + 2]] Sin[ArcCos[coslist[[i*MDIV + 2]]]] Exp[gammalist[[i*MDIV + 2]]]))^2) - 2 NgttTab[[(i)*MDIV + 
                                      1]]/(((rlist[[i*MDIV + 1]] Sin[ArcCos[coslist[[i*MDIV + 1]]]] Exp[gammalist[[i*MDIV + 1]]]))^2))/(((coslist[[(i)*MDIV + 2]] - coslist[[(i)*MDIV + 1]]))^2)))}, {i, 301, SDIV - 3}]; 
      NΩInt = Interpolation[NΩTab]; 
      NΩxInt = Interpolation[NΩxTab]; 
      NΩμInt = Interpolation[NΩμTab]; 
      AvNΩxTab = 
        Table[{NΩxTab[[i, 
                1]], (NΩxTab[[i - 1, 2]] + 
                  2 NΩxTab[[i, 2]] + 
                  NΩxTab[[i + 1, 2]])/4}, {i, 
            2, Dimensions[NΩxTab][[1]] - 1}]; 
      AvNKx2Tab = 
        Table[{NKx2Tab[[i, 
                1]], (NKx2Tab[[i - 1, 2]] + 2 NKx2Tab[[i, 2]] + 
                  NKx2Tab[[i + 1, 2]])/4}, {i, 
            2, Dimensions[NKx2Tab][[1]] - 1}]; 
      AvNKμ2Tab = 
        Table[{NKμ2Tab[[i, 
                1]], (NKμ2Tab[[i - 1, 2]] + 
                  2 NKμ2Tab[[i, 2]] + NKμ2Tab[[i + 1, 2]])/
              4}, {i, 2, Dimensions[NKμ2Tab][[1]] - 1}]; 
      AvNΩxInt = Interpolation[AvNΩxTab]; 
      AvNΩμTab = 
        Table[{NΩμTab[[i, 
                1]], (NΩμTab[[i - 1, 2]] + 
                  2 NΩμTab[[i, 2]] + 
                  NΩμTab[[i + 1, 2]])/4}, {i, 
            2, Dimensions[NΩμTab][[1]] - 1}]; 
      AvNΩμInt = Interpolation[AvNΩμTab]; 
      AvNKxInt = Interpolation[AvNKx2Tab]; 
      AvNKμInt = Interpolation[AvNKμ2Tab];
 (*==========================================================================*) 
overallendtime = 
      SessionTime[]; overalllapsedtime = overallendtime - overallstarttime; 
    Print["total time = ", 
      IntegerPart[overalllapsedtime/(60 60)], " h ", 
      Mod[IntegerPart[overalllapsedtime/60], 60], " min ", 
      Mod[IntegerPart[overalllapsedtime], 60], " sec"];

calculating the specific two-soliton metric

Print["================ the two-soliton parameters \
=================="];
Clear[d, c, a2];a = ageo;"km";M = Mgeo;"km";c=(cf//.Flatten[twosolitonpar2]);"km";a2=(a2f//.Flatten[twosolitonpar2]);"(km^2)";ρmax = 100*Rs;zmax = 100*Rs;b = -c;k = -a2 - a c;Print["b= ", b, "    k= ", k];Print["d = ", (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])];ModelTypeTab[[modIn + 1,3]]=(Sqrt[(((-a b + k))^2) - (b^2) (M^2)]);ModelTypeTab[[modIn + 1,4]]=\(\(1\/2\) \((\(κ\_+\) - \(κ\_-\))\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}\);ModelTypeTab[[modIn + 1,5]]=\(\(1\/2\) \((\(κ\_+\) + \(κ\_-\))\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}\);\(Print[\*""\!\(\(κ\_+\)\) = "", \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\)}]\);\(Print[\*""\!\(\(κ\_-\)\) = "", \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\)}]\);\(Print[{\(R\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(R\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(r\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\), \(r\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\)} //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}]\);If[(((-a b + k))^2) - (b^2) (M^2) > 0,
If[(-(a^2) - (b^2) + 2 d + 2 k + 
                        (M^2) //. {d -> (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])}) > 0 && (-(a^2) - (b^2) - 2 d + 2 k + 
                        (M^2) //. {d -> (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])}) > 0,
Print[
                  StyleForm["Case Ia", FontFamily -> "Times", 
                    FontSize -> 18, FontWeight -> "Bold", 
                    FontColor -> RGBColor[0.3, 0.3, 1]]];ModelTypeTab[[modIn + 1,2]]="Ia";,
If[(-(a^2) - (b^2) + 2 d + 2 k + 
                          (M^2) //. {d -> (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])}) > 0,
Print[
                      StyleForm["Case IIa", 
                        FontFamily -> "Times", FontSize -> 18, 
                        FontWeight -> "Bold", 
                        FontColor -> RGBColor[0.3, 0.3, 1]]];ModelTypeTab[[modIn + 1,2]]="IIa";,
If[(-(a^2) - (b^2) + 2 d + 2 k + 
                              (M^2) //. {d -> (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])}) == 0,Print[
                          StyleForm["Ekfilismos", 
                            FontFamily -> "Times", 
                            FontSize -> 18, 
                            FontWeight -> "Bold", 
                            FontColor -> RGBColor[0.3, 0.3, 1]]];ModelTypeTab[[modIn + 1,2]]="IIb";,
Print[
                          StyleForm["Case IIc", 
                            FontFamily -> "Times", 
                            FontSize -> 18, 
                            FontWeight -> "Bold", 
                            FontColor -> RGBColor[0.3, 0.3, 1]]];ModelTypeTab[[modIn + 1,2]]="IIc";]]],
Print[
              StyleForm["Case III", FontFamily -> "Times", 
                FontSize -> 18, FontWeight -> "Bold", 
                FontColor -> RGBColor[0.3, 0.3, 1]]];ModelTypeTab[[modIn + 1,2]]="III";];Clear[d, c, a2];a = ageo;"km";M = Mgeo;"km";c=(cf//.Flatten[twosolitonpar2]);"km";a2=(a2f//.Flatten[twosolitonpar2]);"(km^2)";ρmax = 100*Rs;zmax = 100*Rs;b = -c;k = -a2 - a c;
Print["Case IIa"];Clear[d];d /: Im[d] = 0;\(ConjA = \(-16\)\ d\ k\ M\^2\ \((Conjugate[\(r\_-\)]\ Conjugate[\(r\_+\)] + 
                Conjugate[\(R\_-\)]\ Conjugate[\(R\_+\)])\) - \((4\ I\ \((a\^2\ b - a\ \((b\^2 + d + k)\) + 
                      b\ \((d + k - 
                            M\^2)\))\)\ Conjugate[\(κ\_-\)]\ \((Conjugate[\(R\_-\)]\ Conjugate[\(r\_+\)] - 
                      Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\) + \((4\ \((a - b)\)\^2\ \((a\ b - d - k)\) - 
                      4\ \((\((a - 2\ b)\)\ b - 
                            d)\)\ M\^2)\)\ \((Conjugate[\(R\_-\)]\ Conjugate[\(r\_+\)] + 
                      Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\))\)\ Conjugate[\(κ\_+\)]\^2 + 
          Conjugate[\(κ\_-\)]\^2\ \((\((4\ \((a - b)\)\^2\ \((a\ b + 
                            d - k)\) - 
                      4\ \((\((a - 2\ b)\)\ b + 
                            d)\)\ M\^2)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] + Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\) + 
                4\ I\ \((\((a - b)\)\ \((a\ b + d - k)\) - 
                      b\ M\^2)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] - Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)])\)\);\(ConjB = 
        Conjugate[\(κ\_-\)]\ Conjugate[\(κ\_+\)]\ \((2\ I\ M\ \((Conjugate[\(R\_-\)] - 
                    Conjugate[\(R\_+\)])\)\ \((\(-Conjugate[\(κ\_-\)]\) + 
                    Conjugate[\(κ\_+\)])\)\ \((\(-a\^2\)\ b + 
                    2\ a\ k + b\ \((b\^2 - 2\ k + M\^2)\) - 
                    b\ Conjugate[\(κ\_-\)]\ Conjugate[\(κ\_+\)])\) - 
              2\ I\ M\ \((Conjugate[\(r\_-\)] - 
                    Conjugate[\(r\_+\)])\)\ \((Conjugate[\(κ\_-\)] + 
                    Conjugate[\(κ\_+\)])\)\ \((\(-a\^2\)\ b + 
                    2\ a\ k + b\ \((b\^2 - 2\ k + M\^2)\) + 
                    b\ Conjugate[\(κ\_-\)]\ Conjugate[\(κ\_+\)])\) + 
              1\/2\ \((8\ d\ M\ \((\(-a\^2\) + b\^2 + 
                          M\^2)\)\ \((\(-Conjugate[\(r\_-\)]\) + 
                          Conjugate[\(R\_-\)] - Conjugate[\(r\_+\)] + 
                          Conjugate[\(R\_+\)])\) - 
                    8\ d\ M\ Conjugate[\(κ\_-\)]\ \((Conjugate[\(r\_-\)] + Conjugate[\(R\_-\)] + Conjugate[\(r\_+\)] + 
                          Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)])\))\)\);\(ConjA1 = 
        ConjA //. {Conjugate[\(κ\_-\)] -> \(-\(κ\_-\)\), 
            Conjugate[\(κ\_+\)] -> \(κ\_+\), 
            Conjugate[\(r\_-\)] -> \(R\_-\), 
            Conjugate[\(R\_-\)] -> \(r\_-\), 
            Conjugate[\(r\_+\)] -> \(R\_+\), 
            Conjugate[\(R\_+\)] -> \(r\_+\)}\);\(ConjB1 = 
        ConjB //. {Conjugate[\(κ\_-\)] -> \(-\(κ\_-\)\), 
            Conjugate[\(κ\_+\)] -> \(κ\_+\), 
            Conjugate[\(r\_-\)] -> \(R\_-\), 
            Conjugate[\(R\_-\)] -> \(r\_-\), 
            Conjugate[\(r\_+\)] -> \(R\_+\), 
            Conjugate[\(R\_+\)] -> \(r\_+\)}\);ConjKo = Conjugate[Ko];\(ConjKo1 = 
        ConjKo //. {Conjugate[\(κ\_-\)] -> \(-\(κ\_-\)\), 
            Conjugate[\(κ\_+\)] -> \(κ\_+\), 
            Conjugate[\(r\_-\)] -> \(R\_-\), 
            Conjugate[\(R\_-\)] -> \(r\_-\), 
            Conjugate[\(r\_+\)] -> \(R\_+\), 
            Conjugate[\(R\_+\)] -> \(r\_+\)}\);\(ConjL = \(-4\)\ I\ d\ \((a - I\ M)\)\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\) - 
                2\ \((a\^2 - b\^2 - 2\ k)\)\ M\^2 + 
                M\^4)\)\ \((Conjugate[\(r\_-\)] + Conjugate[\(R\_-\)] + 
                Conjugate[\(r\_+\)] + Conjugate[\(R\_+\)])\) - 
          4\ I\ \((I\ b\^4\ M\^2 - I\ b\^2\ M\^2\ \((d + 2\ k - 
                      M\^2)\) - I\ M\^2\ \((2\ k\ \((d + k)\) + 
                      d\ M\^2)\) + a\^4\ b\ \((M - z)\) + 
                a\^3\ \((b\^2 - d - k)\)\ \((M - z)\) + 
                b\^3\ \((d + k + M\^2)\)\ \((M - z)\) - 
                b\ \((4\ k\ \((d + k)\) - \((d + k)\)\ M\^2 - M\^4)\)\ \((M - 
                      z)\) + a\ \((4\ I\ b\ k\ M\^2 + 
                      b\^2\ \((d + 5\ k - M\^2)\)\ \((M - z)\) + \((d + 
                            k)\)\ \((4\ k + M\^2)\)\ \((M - z)\) + 
                      b\^4\ \((\(-M\) + z)\))\) + 
                a\^2\ \((\(-I\)\ b\^2\ M\^2 + I\ d\ M\^2 + b\^3\ \((\(-M\) + z)\) + 
                      b\ \((d + 5\ k + 2\ M\^2)\)\ \((\(-M\) + 
                            z)\))\))\)\ Conjugate[\(κ\_-\)]\ \((\(-Conjugate[\(R\_-\)]\)\ Conjugate[\(r\_+\)] + 
                Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\) + 
          4\ \((\(-I\)\ b\^5\ M\^2 + I\ b\^3\ \((3\ k\ M\^2 - 2\ M\^4)\) - I\ b\ M\^2\ \((\(-2\)\ k\ \((d + k)\) + 
                      3\ k\ M\^2 + M\^4)\) + 
                a\^3\ \((I\ k\ M\^2 - 2\ b\^3\ \((M - z)\) - 
                      2\ b\ \((2\ k + M\^2)\)\ \((M - z)\))\) + 
                a\^2\ \((2\ I\ b\^3\ M\^2 + I\ b\ \((3\ k\ M\^2 + 2\ M\^4)\) + 
                      2\ b\^2\ \((d + 5\ k)\)\ \((M - 
                            z)\) + \((4\ k\ \((d + k)\) + \((2\ d + 
                                  k)\)\ M\^2)\)\ \((M - z)\))\) + 
                a\ \((\(-7\)\ I\ b\^2\ k\ M\^2 - I\ k\ M\^2\ \((2\ \((d + k)\) + M\^2)\) + b\^5\ \((M - z)\) - 
                      2\ b\^3\ \((2\ k - M\^2)\)\ \((M - z)\) - 
                      b\ \((8\ k\ \((d + k)\) - 4\ k\ M\^2 - M\^4)\)\ \((M - 
                            z)\))\) + a\^5\ b\ \((M - z)\) - 
                b\^4\ \((d + k)\)\ \((M - z)\) - 
                M\^2\ \((2\ k\ \((d + k)\) + d\ M\^2)\)\ \((M - z)\) + 
                b\^2\ \((4\ k\ \((d + k)\) - \((2\ d + 
                            3\ k)\)\ M\^2)\)\ \((M - z)\) + 
                a\^4\ \((\(-M\)\ \((d + k + I\ b\ M)\) + \((d + 
                            k)\)\ z)\))\)\ \((Conjugate[\(R\_-\)]\ Conjugate[\(r\_+\)] + Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\) - 
          16\ d\ k\ M\^2\ \((\(-I\)\ a + I\ b + M - 
                z)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(r\_+\)] + 
                Conjugate[\(R\_-\)]\ Conjugate[\(R\_+\)])\) - 
          4\ \((\(-I\)\ b\^5\ M\^2 + I\ b\^3\ \((3\ k\ M\^2 - 2\ M\^4)\) - I\ b\ M\^2\ \((2\ \((d - k)\)\ k + 
                      3\ k\ M\^2 + M\^4)\) + 
                a\^3\ \((I\ k\ M\^2 - 2\ b\^3\ \((M - z)\) - 
                      2\ b\ \((2\ k + M\^2)\)\ \((M - z)\))\) + 
                a\ \((\(-7\)\ I\ b\^2\ k\ M\^2 - I\ k\ M\^2\ \((\(-2\)\ d + 2\ k + M\^2)\) + b\^5\ \((M - z)\) - 
                      2\ b\^3\ \((2\ k - M\^2)\)\ \((M - z)\) + 
                      b\ \((8\ \((d - k)\)\ k + 4\ k\ M\^2 + M\^4)\)\ \((M - 
                            z)\))\) + 
                a\^2\ \((2\ I\ b\^3\ M\^2 + I\ b\ \((3\ k\ M\^2 + 2\ M\^4)\) - 
                      2\ b\^2\ \((d - 5\ k)\)\ \((M - 
                            z)\) - \((2\ d\ \((2\ k + M\^2)\) - 
                            k\ \((4\ k + M\^2)\))\)\ \((M - z)\))\) + 
                a\^5\ b\ \((M - z)\) + b\^4\ \((d - k)\)\ \((M - z)\) + 
                M\^2\ \((2\ \((d - k)\)\ k + d\ M\^2)\)\ \((M - z)\) - 
                b\^2\ \((4\ d\ k - 4\ k\^2 - 2\ d\ M\^2 + 
                      3\ k\ M\^2)\)\ \((M - z)\) + 
                a\^4\ \((\(-I\)\ b\ M\^2 + d\ \((M - z)\) + 
                      k\ \((\(-M\) + 
                            z)\))\))\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] + Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\) - 
          4\ d\ M\ \((\(-I\)\ a\^3 - 2\ I\ b\ k - 
                a\^2\ M + b\^2\ M + 
                M\^3 + I\ a\ \((b\^2 + 2\ k + 
                      M\^2)\))\)\ Conjugate[\(κ\_-\)]\ \((Conjugate[\(r\_-\)] - Conjugate[\(R\_-\)] + Conjugate[\(r\_+\)] - 
                Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)] - 
          4\ I\ \((I\ b\^4\ M\^2 + I\ b\^2\ M\^2\ \((d - 2\ k + 
                      M\^2)\) + I\ M\^2\ \((2\ \((d - k)\)\ k + 
                      d\ M\^2)\) + a\^4\ b\ \((M - z)\) + 
                a\^3\ \((b\^2 + d - k)\)\ \((M - z)\) - 
                b\^3\ \((d - k - M\^2)\)\ \((M - z)\) + 
                b\ \((4\ \((d - k)\)\ k - \((d - k)\)\ M\^2 + M\^4)\)\ \((M - 
                      z)\) + a\^2\ \((\(-I\)\ b\^2\ M\^2 - I\ d\ M\^2 + b\ \((d - 5\ k - 2\ M\^2)\)\ \((M - z)\) + 
                      b\^3\ \((\(-M\) + z)\))\) + 
                a\ \((4\ I\ b\ k\ M\^2 - \((d - k)\)\ \((4\ k + 
                            M\^2)\)\ \((M - z)\) + b\^4\ \((\(-M\) + z)\) + 
                      b\^2\ \((d - 5\ k + M\^2)\)\ \((\(-M\) + 
                            z)\))\))\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] - Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)] + 4\ M\ \((Conjugate[\(r\_-\)] - 
                Conjugate[\(r\_+\)])\)\ \((\((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 - 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) - I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 - 
                                  d)\) - \((a\^2 + 3\ a\ b - 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 - d + k)\) - 
                            a\ b\ \((b\^2 + d + 
                                  3\ k)\))\)\ M\^2 + I\ \((2\ a\^2\ b - a\ k - b\ \((b\^2 + d + 2\ k)\))\)\ M\^3 + 
                      a\ b\ M\^4 - I\ b\ M\^5)\)\ Conjugate[\(κ\_-\)] + \((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 + 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) - I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 + 
                                  d)\) - \((a\^2 + 3\ a\ b + 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 + d + k)\) - 
                            a\ b\ \((b\^2 - d + 
                                  3\ k)\))\)\ M\^2 - I\ \((\(-2\)\ a\^2\ b + b\^3 - b\ d + a\ k + 2\ b\ k)\)\ M\^3 + 
                      a\ b\ M\^4 - I\ b\ M\^5)\)\ Conjugate[\(κ\_+\)])\) + 
          4\ M\ \((Conjugate[\(R\_-\)] - 
                Conjugate[\(R\_+\)])\)\ \((\((\(-\((a - b)\)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 - 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) + I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 - 
                                  d)\) - \((a\^2 + 3\ a\ b - 2\ d)\)\ k + 
                            2\ k\^2)\)\ M + \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 - d + k)\) - 
                            a\ b\ \((b\^2 + d + 
                                  3\ k)\))\)\ M\^2 - I\ \((2\ a\^2\ b - a\ k - b\ \((b\^2 + d + 2\ k)\))\)\ M\^3 - 
                      a\ b\ M\^4 + I\ b\ M\^5)\)\ Conjugate[\(κ\_-\)] + \((\((a - b)\)\ \((a\ b - 
                            k)\)\ \((\((a + b)\)\ \((a\^2 + 
                                  d)\) + \((\(-3\)\ a + 
                                  b)\)\ k)\) - I\ \((a - 
                            b)\)\ \((b\ \((a + b)\)\ \((a\^2 + 
                                  d)\) - \((a\^2 + 3\ a\ b + 2\ d)\)\ k + 
                            2\ k\^2)\)\ M - \((2\ a\^3\ b - a\^2\ k + 
                            k\ \((2\ b\^2 + d + k)\) - 
                            a\ b\ \((b\^2 - d + 
                                  3\ k)\))\)\ M\^2 - I\ \((\(-2\)\ a\^2\ b + b\^3 - b\ d + a\ k + 2\ b\ k)\)\ M\^3 + 
                      a\ b\ M\^4 - I\ b\ M\^5)\)\ Conjugate[\(κ\_+\)])\)\);\(ConjE = 
        4\ d\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\) - 
                2\ \((a\^2 - b\^2 - 2\ k)\)\ M\^2 + 
                M\^4)\)\ \((\(-I\)\ a + M - 
                z)\)\ \((Conjugate[\(r\_-\)] + Conjugate[\(R\_-\)] + 
                Conjugate[\(r\_+\)] + Conjugate[\(R\_+\)])\) + 
          4\ I\ M\ \((\((a - b)\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b - d - 
                      k)\) + I\ \((\((a - b)\)\ \((a + 
                            b)\)\ \((b\^2 - d)\) + 
                      2\ \((\(-2\)\ a\ b + b\^2 + d)\)\ k + 
                      2\ k\^2)\)\ M - \((a + b)\)\ \((2\ a\ b - b\^2 - d - 
                      k)\)\ M\^2 - I\ \((b\^2 - d)\)\ M\^3 + 
                b\ M\^4)\)\ Conjugate[\(κ\_-\)]\ \((\(-Conjugate[\(R\_-\)]\)\ Conjugate[\(r\_+\)] + 
                Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\) + 
          4\ M\ \((\(-\((a - b)\)\^2\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b - d - k)\) - I\ \((a - 
                      b)\)\ \((\((a - b)\)\ b\ \((a + b)\)\^2 - \((a\^2 + 
                            4\ a\ b - 3\ b\^2 - 2\ d)\)\ k + 
                      2\ k\^2)\)\ M + \((2\ a\^3\ b + 2\ k\ \((d + k)\) - 
                      a\^2\ \((2\ d + k)\) + b\^2\ \((2\ d + 3\ k)\) - 
                      2\ a\ \((b\^3 + 
                            2\ b\ k)\))\)\ M\^2 + I\ \((2\ a\^2\ b - 2\ b\^3 - a\ k - 3\ b\ k)\)\ M\^3 - \((a\ b - 
                      d)\)\ M\^4 - I\ b\ M\^5)\)\ \((Conjugate[\(R\_-\)]\ Conjugate[\(r\_+\)] + 
                Conjugate[\(r\_-\)]\ Conjugate[\(R\_+\)])\) + 
          16\ d\ k\ M\^2\ \((I\ a - I\ b + 
                M)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(r\_+\)] + 
                Conjugate[\(R\_-\)]\ Conjugate[\(R\_+\)])\) + 
          4\ M\ \((\((a - b)\)\^2\ \((\((a + b)\)\^2 - 4\ k)\)\ \((a\ b + d - 
                      k)\) + I\ \((a - 
                      b)\)\ \((\((a - b)\)\ b\ \((a + b)\)\^2 - \((a\^2 + 
                            4\ a\ b - 3\ b\^2 + 2\ d)\)\ k + 
                      2\ k\^2)\)\ M - \((2\ \((a - b)\)\ \((a + 
                            b)\)\ \((a\ b + d)\) - \((a\^2 + 4\ a\ b - 
                            3\ b\^2 + 2\ d)\)\ k + 
                      2\ k\^2)\)\ M\^2 - I\ \((2\ a\^2\ b - 
                      2\ b\^3 - a\ k - 3\ b\ k)\)\ M\^3 + \((a\ b + 
                      d)\)\ M\^4 + I\ b\ M\^5)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] + 
                Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\) + 
          4\ d\ M\ \((I\ a\^3 + 2\ I\ b\ k + 
                b\^2\ M + 
                M\^3 - I\ a\ \((b\^2 + 2\ k + M\^2)\) - \((b\^2 + 
                      M\^2)\)\ z + 
                a\^2\ \((\(-M\) + 
                      z)\))\)\ Conjugate[\(κ\_-\)]\ \((Conjugate[\(r\_-\)] - Conjugate[\(R\_-\)] + Conjugate[\(r\_+\)] - 
                Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)] + 
          4\ I\ M\ \((\((a - b)\)\ \((\((a + b)\)\^2 - 
                      4\ k)\)\ \((a\ b + d - 
                      k)\) + I\ \((\((a - b)\)\ \((a + 
                            b)\)\ \((b\^2 + d)\) - 
                      2\ \((2\ a\ b - b\^2 + d)\)\ k + 2\ k\^2)\)\ M - \((a + 
                      b)\)\ \((2\ a\ b - b\^2 + d - 
                      k)\)\ M\^2 - I\ \((b\^2 + d)\)\ M\^3 + 
                b\ M\^4)\)\ \((Conjugate[\(r\_-\)]\ Conjugate[\(R\_-\)] - 
                Conjugate[\(r\_+\)]\ Conjugate[\(R\_+\)])\)\ Conjugate[\(κ\_+\)] - 
          4\ M\ \((Conjugate[\(R\_-\)] - 
                Conjugate[\(R\_+\)])\)\ \((\((a\^5\ b + \((d - 
                            k)\)\ k\ M\^2 + 
                      b\^2\ k\ \((\(-d\) + k - 2\ M\^2)\) - 
                      a\^4\ \((k - I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((d + 3\ k + 
                                  2\ M\^2)\) + I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d - k + M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 + \((d + 3\ k)\)\ M\^2 + 
                                  M\^4)\) + 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) + I\ k\ \((2\ \((d + 
                                        k)\) + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((d + 3\ k + 
                                  M\^2)\) - I\ b\^3\ \((M - 
                                  z)\) - I\ b\ \((d + 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) + I\ b\^3\ \((d + 
                            M\^2)\)\ \((M - 
                            z)\) - I\ b\ \((2\ k\ \((d + 
                                  k)\) - \((d + 2\ k)\)\ M\^2 - 
                            M\^4)\)\ \((M - 
                            z)\))\)\ Conjugate[\(κ\_-\)] + \((\(-a\^5\)\ b + k\ \((d + k)\)\ M\^2 - b\^2\ k\ \((d + k - 2\ M\^2)\) + 
                      a\^4\ \((k - I\ b\ \((M - z)\))\) + 
                      a\^3\ \((b\^3 + 
                            b\ \((\(-d\) + 3\ k + 
                                  2\ M\^2)\) + I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d + k - M\^2)\) + 
                            b\ \((d\ M\^2 + \((k - M\^2)\)\ \((4\ k + 
                                        M\^2)\))\) - 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) - I\ k\ \((\(-2\)\ d + 
                                  2\ k + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((\(-5\)\ b\^2\ k + 
                            k\ \((d - 3\ k - 
                                  M\^2)\) + I\ b\^3\ \((M - 
                                  z)\) - I\ b\ \((d - 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) + I\ b\^3\ \((d - 
                            M\^2)\)\ \((M - 
                            z)\) - I\ b\ \((2\ \((d - 
                                  k)\)\ k - \((d - 2\ k)\)\ M\^2 + 
                            M\^4)\)\ \((M - 
                            z)\))\)\ Conjugate[\(κ\_+\)])\) + 
          4\ M\ \((Conjugate[\(r\_-\)] - 
                Conjugate[\(r\_+\)])\)\ \((\((a\^5\ b + \((d - 
                            k)\)\ k\ M\^2 + 
                      b\^2\ k\ \((\(-d\) + k - 2\ M\^2)\) - 
                      a\^4\ \((k - I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((d + 3\ k + 
                                  2\ M\^2)\) + I\ k\ \((M - 
                                  z)\))\) + 
                      a\ \((b\^3\ \((d - k + M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 + \((d + 3\ k)\)\ M\^2 + 
                                  M\^4)\) + 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) + I\ k\ \((2\ \((d + 
                                        k)\) + M\^2)\)\ \((M - z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((d + 3\ k + 
                                  M\^2)\) - I\ b\^3\ \((M - 
                                  z)\) - I\ b\ \((d + 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) + I\ b\^3\ \((d + 
                            M\^2)\)\ \((M - 
                            z)\) - I\ b\ \((2\ k\ \((d + 
                                  k)\) - \((d + 2\ k)\)\ M\^2 - 
                            M\^4)\)\ \((M - 
                            z)\))\)\ Conjugate[\(κ\_-\)] + \((a\^5\ b - k\ \((d + k)\)\ M\^2 + b\^2\ k\ \((d + k - 2\ M\^2)\) - 
                      a\^4\ \((k - I\ b\ \((M - z)\))\) - 
                      a\^3\ \((b\^3 + 
                            b\ \((\(-d\) + 3\ k + 
                                  2\ M\^2)\) + I\ k\ \((M - 
                                  z)\))\) + 
                      a\^2\ \((5\ b\^2\ k + 
                            k\ \((\(-d\) + 3\ k + 
                                  M\^2)\) - I\ b\^3\ \((M - 
                                  z)\) + I\ b\ \((d - 
                                  2\ \((k + M\^2)\))\)\ \((M - 
                                  z)\))\) - I\ b\^3\ \((d - 
                            M\^2)\)\ \((M - 
                            z)\) + I\ b\ \((2\ \((d - 
                                  k)\)\ k - \((d - 2\ k)\)\ M\^2 + 
                            M\^4)\)\ \((M - z)\) + 
                      a\ \((\(-b\^3\)\ \((d + k - M\^2)\) + 
                            b\ \((\(-4\)\ k\^2 - \((d - 3\ k)\)\ M\^2 + 
                                  M\^4)\) + 
                            3\ I\ b\^2\ k\ \((M - 
                                  z)\) - I\ k\ \((\(-2\)\ d + 
                                  2\ k + M\^2)\)\ \((\(-M\) + 
                                  z)\))\))\)\ Conjugate[\(κ\_+\)])\)\);\(ConjL1 = 
        ConjL //. {Conjugate[\(κ\_-\)] -> \(-\(κ\_-\)\), 
            Conjugate[\(κ\_+\)] -> \(κ\_+\), 
            Conjugate[\(r\_-\)] -> \(R\_-\), 
            Conjugate[\(R\_-\)] -> \(r\_-\), 
            Conjugate[\(r\_+\)] -> \(R\_+\), 
            Conjugate[\(R\_+\)] -> \(r\_+\)}\);\(ConjE1 = 
        ConjE //. {Conjugate[\(κ\_-\)] -> \(-\(κ\_-\)\), 
            Conjugate[\(κ\_+\)] -> \(κ\_+\), 
            Conjugate[\(r\_-\)] -> \(R\_-\), 
            Conjugate[\(R\_-\)] -> \(r\_-\), 
            Conjugate[\(r\_+\)] -> \(R\_+\), 
            Conjugate[\(R\_+\)] -> \(r\_+\)}\);\(f = \(\(A\ ConjA1 - B\ ConjB1\)\/\(\((A + B)\)\ \ \((ConjA1 + ConjB1)\)\) //. {\(R\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(R\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(r\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\), \(r\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\)}\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\)}\);\(Exp2γ = \(\(A\ ConjA1 - B\ ConjB1\)\/\(Ko\ ConjKo1\ \(R\_-\)\ \(r\_+\)\ \(r\_-\)\ \(R\_+\)\) //. {\(R\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(R\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(r\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\), \(r\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\)}\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\)}\);ω=-(-I (((A + B) ConjL1  - (ConjA1 + 
                                    ConjB1) Ε) - ((ConjA1 + ConjB1) L  - (A + B) ConjE1)))/(2 (A ConjA1  - B ConjB1))//.\({\(R\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(R\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) + \(κ\_-\))\))\)\^2\), \(r\_+\) -> \@\(ρ\^2 + \((z + \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\), \(r\_-\) -> \@\(ρ\^2 + \((z - \(1\/2\) \((\(κ\_+\) - \(κ\_-\))\))\)\^2\)}\)//.{d -> (Sqrt[(((-a b + k))^2) - (b^2) (M^2)])};
\(gtt[ρ_, 
          z_] = \(-f\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}; 
      gtφ[ρ_, z_] = 
        f\ ω //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}; 
      gφφ[ρ_, 
          z_] = \(-f\)\ ω\^2 + ρ\^2\ f\^\(-1\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)}; 
      gzz[ρ_, z_] = 
        Exp2γ\ f\^\(-1\) //. {d -> \@\(\((\(-a\)\ b + k)\)\^2 - b\^2\ M\^2\), \(κ\_+\)\  -> \@\((\(-a\^2\) - b\^2 + 2\ d + 2\ k + M\^2)\), \(κ\_-\)\  -> \@\((\(-a\^2\) - b\^2 - 2\ d + 2\ k + M\^2)\)};\)
  
Print["=================== CTCs and SLSs \
====================="];
Timing[(Subscript[plot2SLS, modIn]) = 
            ContourPlot[
              Chop[gtt[ρ, z]], {ρ, 0, 3 Mgeo}, {z, -3 Mgeo, 
                3 Mgeo}, ContourShading -> False, Contours -> {0}, 
              FrameLabel -> {"ρ", "z"}, 
              ContourStyle -> Red, PlotPoints -> 80, 
              PlotLabel -> "j= " <> ToString[a/M]];]
  
Timing[(Subscript[plot2CTC, modIn]) = 
          ContourPlot[
            Chop[gφφ[ρ, z]], {ρ, 0, 
              3 Mgeo}, {z, -3 Mgeo, 3 Mgeo}, 
            ContourShading -> False, Contours -> {0}, 
            PlotPoints -> 80, 
            FrameLabel -> {"ρ", "z"}, 
            PlotLabel -> "j= " <> 
                ToString[a/M]];]
Show[
        (Subscript[plot2SLS, modIn]), (Subscript[plot2CTC, modIn])];

calculating the corresponding Manko et al. spacetimes

Print["Manko et al."];
\(Clear[
      a, M, a2, c, b]; a = ageo; "km"; M = Mgeo; "km"; 
    c = \((cf //. Flatten[twosolitonpar2])\); "km"; 
    a2 = \((a2f //. Flatten[twosolitonpar2])\); \*""\!\(km\^2\)""; 
    J = a\ M; 
    M2 = 1\/4\ M\ \((\(-3\)\ a\^2 + 2\ a\ b - M\^2 + 
            b\^2\ \((1 + \(4\ M\^2\)\/\(\((a - b)\)\^2 - M\^2\))\))\);\)
Qplot=\(Plot[{M2New, M2}, {b, \(-10\), 10}, 
          PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, 
          Axes -> False, Frame -> True, 
          FrameLabel -> {"b", \*""\!\(M\_2\)""}]\);
Clear[QuadSol];QuadSol = Simplify[Solve[M2 == M2New, b]];Off[$RecursionLimit::"reclim"];Clear[k, z, ρ, x, y];Clear[a, M, a2, c, b, Val];a = ageo;"km";M = Mgeo;"km";J = a M;j = 3;For[i = 1, i < 5, 
        If[Im[b //. QuadSol[[i]]] == 0, 
          (Subscript[Val, j]) = b //. QuadSol[[i]]; j = j + 1, (Subscript[Val, j]) = Null]; 
        i = i + 1];If[NumberQ[(Subscript[Val, 3])], Null, (Subscript[Val, 4]) = Null];ModelTypeTab[[modIn + 1,7]]=NumberQ[(Subscript[Val, 3])];If[NumberQ[(Subscript[Val, 3])],ModelTypeTab[[modIn + 1,8]]=(Subscript[Val, 3]),ModelTypeTab[[modIn + 1,8]]=Nill];Clear[j];b = (Subscript[Val, 3]);δ = ((-(M^2) (b^2))/((M^2) - (((a - b))^2)));d = (1/4) ((M^2) - (((a - b))^2));kk = (Sqrt[d + δ]);{x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
        y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};EM = (((4 ((((kk^2) ((x^2) - 1) + δ (1 - (y^2))))^2) + (a - b) ((a - b) (d - δ) - (M^2) b) (((1 - (y^2)))^2)))^2) - 
            16 (kk^2) ((x^2) - 1) (1 - 
                  (y^2)) ((((a - b) ((kk^2) ((x^2) - (y^2)) + 2 δ (y^2)) + (M^2) b (y^2)))^2) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};DM = (((4 ((((kk^2) (x^2) - δ (y^2)))^2) + 2 kk M x (2 (kk^2) ((x^2) - 1) + (2 δ + a b - (b^2)) (1 - (y^2))) + (a - b) ((a - b) (d - δ) - (M^2) b) ((y^4) - 1) - 4 (d^2)))^2) + 
            4 (y^2) (((2 (kk^2) ((x^2) - 1) (kk x (a - b) - M b) - 2 M b δ (1 - (y^2)) + ((a - b) (d - δ) - (M^2) b) (2 kk x + M) (1 - (y^2))))^2) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};FM = 
        8 (kk^2) ((x^2) - 
                  1) ((a - b) ((kk^2) ((x^2) - (y^2)) + 
                        2 δ (y^2)) + 
                  (y^2) (M^2) b) (kk M x ((((2 kk x + M))^2) - 
                        2 (y^2) (2  δ + a b - (b^2)) - (a^2) + 
                        (b^2)) - 
                  2 (y^2) (4 δ d - 
                        (M^2) (b^2))) + (4 ((((kk^2) ((x^2) - 1) + δ (1 - (y^2))))^2) + (a - 
                        b) ((a - b) (d - δ) - 
                        (M^2) b) (((1 - (y^2)))^2)) (4 (2 kk M b x + 2 (M^2) b) ((kk^2) ((x^2) - 1) + δ (1 - 
                              (y^2))) + (1 - 
                        (y^2)) ((a - b) ((M^2) (b^2) - 
                              4 δ d) - (4 kk M x + 
                              2 (M^2)) ((a - b) (d - δ) - 
                              (M^2) b))) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};fM3 = EM/DM;ωM3 = -(1 - (y^2)) FM/
              EM /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};Exp2γM3 = 
        EM/(16 (kk^8) ((((x^2) - (y^2)))^4)) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
            y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))};Mgtt3[ρ_, z_] = -fM3;Mgtφ3[ρ_, z_] = fM3 ωM3;Mgφφ3[ρ_, 
          z_] = -fM3 ω(M3^2) + (ρ^2) (fM3^(-1));Mgzz3[ρ_, z_] = (fM3^(-1)) Exp2γM3;
Clear[k, z, ρ, x, y]; 
    b = (Subscript[Val, 4]); δ = ((-(M^2) (b^2))/((M^2) - (((a - b))^2))); 
    d = (1/4) ((M^2) - (((a - b))^2)); 
    kk = (Sqrt[d + δ]); {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
      y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; 
    EM = (((4 ((((kk^2) ((x^2) - 1) + δ (1 - (y^2))))^2) + (a - b) ((a - b) (d - δ) - (M^2) b) (((1 - (y^2)))^2)))^2) - 
          16 (kk^2) ((x^2) - 1) (1 - 
                (y^2)) ((((a - b) ((kk^2) ((x^2) - (y^2)) + 2 δ (y^2)) + (M^2) b (y^2)))^2) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
          y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; 
    DM = (((4 ((((kk^2) (x^2) - δ (y^2)))^2) + 2 kk M x (2 (kk^2) ((x^2) - 1) + (2 δ + a b - (b^2)) (1 - (y^2))) + (a - b) ((a - b) (d - δ) - (M^2) b) ((y^4) - 1) - 4 (d^2)))^2) + 
          4 (y^2) (((2 (kk^2) ((x^2) - 1) (kk x (a - b) - M b) - 2 M b δ (1 - (y^2)) + ((a - b) (d - δ) - (M^2) b) (2 kk x + M) (1 - (y^2))))^2) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; 
    FM = 8 (kk^2) ((x^2) - 
                1) ((a - b) ((kk^2) ((x^2) - (y^2)) + 
                      2 δ (y^2)) + 
                (y^2) (M^2) b) (kk M x ((((2 kk x + M))^2) - 
                      2 (y^2) (2  δ + a b - (b^2)) - (a^2) + 
                      (b^2)) - 
                2 (y^2) (4 δ d - 
                      (M^2) (b^2))) + (4 ((((kk^2) ((x^2) - 1) + δ (1 - (y^2))))^2) + (a - 
                      b) ((a - b) (d - δ) - 
                      (M^2) b) (((1 - (y^2)))^2)) (4 (2 kk M b x + 2 (M^2) b) ((kk^2) ((x^2) - 1) + δ (1 - 
                            (y^2))) + (1 - 
                      (y^2)) ((a - b) ((M^2) (b^2) - 
                            4 δ d) - (4 kk M x + 
                            2 (M^2)) ((a - b) (d - δ) - 
                            
                            (M^2) b))) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
          y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; 
    fM4 = EM/DM; ωM4 = -(1 - (y^2)) FM/
            EM /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
          y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; 
    Exp2γM4 = 
      EM/(16 (kk^8) ((((x^2) - (y^2)))^4)) /. {x -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) + (Sqrt[(((kk - z))^2) + (ρ^2)])), 
          y -> (1/(2 kk)) ((Sqrt[(((kk + z))^2) + (ρ^2)]) - (Sqrt[(((kk - z))^2) + (ρ^2)]))}; Mgtt4[ρ_, z_] = -fM4; 
    Mgtφ4[ρ_, z_] = fM4 ωM4; 
    Mgφφ4[ρ_, 
        z_] = -fM4 ω(M4^2) + (ρ^2) (fM4^(-1)); 
    Mgzz4[ρ_, z_] = (fM4^(-1)) Exp2γM4; 
    Clear[k, z, ρ, x, 
      y];
\(Print[\*""\!\(\(b\_-\)\) = "", 
        Val\_3, \*"" ,   \!\(\(b\_+\)\) = "", 
        Val\_4];\)
Off[
        Graphics::"gptn"];

comparing the metric functions between the numerical spacetime and the two-soliton, the Hartle-Thorne and the two corresponding Manko et al. solutions

Print["============== comparing the metric function gtt on \
the equatorial plane ============="];
Timing[StestVal = 0;test = Table[{0, 0}, {i, 1, 300}];test2 = Table[{0, 0, 0}, {i, 1, 300}];testM3 = Table[{0, 0}, {i, 1, 300}];testM4 = Table[{0, 0}, {i, 1, 300}];testHT = Table[{0, 0, 0}, {i, 1, 300}];For[i = 1, i < 301, 
          test[[i, 1]] = Log[10, Nflist[[301 + i, 1]]]; 
          test2[[i, 1]] = Nflist[[301 + i, 1]]; 
          testHT[[i, 1]] = Log[10, Nflist[[301 + i, 1]]]; 
          rHT = FindRoot[
              Nflist[[301 + i, 
                    1]] == (((gtφHT[Mgeo, Jgeo/(Mgeo^2), r, 0]^2 - gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] gφφHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0]))^(1/2)), {r, Rs}]; 
          If[NumberQ[
              (Subscript[Val, 3])], testM3[[i, 1]] = 
                Log[10, Nflist[[301 + i, 1]]];, Null]; 
          If[NumberQ[
              (Subscript[Val, 4])], testM4[[i, 1]] = 
                Log[10, Nflist[[301 + i, 1]]];, Null]; 
          test[[i, 2]] = 
            Log[Abs[(Nflist[[301 + i, 2]] + 
                        Chop[gtt[Nflist[[301 + i, 1]], 0]])/
                    Nflist[[301 + i, 2]]]]/Log[10]; 
          test2[[i, 2]] = Chop[gtt[Nflist[[301 + i, 1]], 0]]; 
          test2[[i, 3]] = -Nflist[[301 + i, 2]]; 
          If[i < 300, 
            StestVal = 
              StestVal + ((((Nflist[[301 + i, 2]] + Chop[gtt[Nflist[[301 + i, 1]], 0]]))^2) + (((Nflist[[301 + i + 1, 2]] + Chop[gtt[Nflist[[301 + i + 1, 1]], 0]]))^2)) (Nflist[[301 + i + 1, 1]] - Nflist[[301 + i, 1]])/2, Null]; 
          testHT[[i, 2]] = 
            Log[Abs[(Nflist[[301 + i, 2]] + 
                          gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] //. 
                        rHT)/Nflist[[301 + i, 2]]]]/Log[10]; 
          testHT[[i, 3]] = Log[Abs[r //. rHT]]/Log[10]; 
          If[NumberQ[
              (Subscript[Val, 3])], testM3[[i, 2]] = 
                Log[10, 
                  Abs[(Nflist[[301 + i, 2]] + 
                          Chop[Mgtt3[Nflist[[301 + i, 1]], 0]])/
                      Nflist[[301 + i, 2]]]];, Null]; 
          If[NumberQ[
              (Subscript[Val, 4])], testM4[[i, 2]] = 
                Log[10, 
                  Abs[(Nflist[[301 + i, 2]] + 
                          Chop[Mgtt4[Nflist[[301 + i, 1]], 0]])/
                      Nflist[[301 + i, 2]]]];, Null]; i = i + 1];If[NumberQ[(Subscript[Val, 3])], Null, Null];If[NumberQ[(Subscript[Val, 4])], Null, Null];Off[Graphics::"gptn"];ModelTypeTab[[modIn + 1,24]]=(Sqrt[StestVal]);Print["σ_tt"];Print[(Sqrt[StestVal])];If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])],plot1=Show[ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 1, 
                        300}], DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ListPlot[testM3, DisplayFunction -> Identity, 
                    PlotStyle -> Magenta],
ListPlot[test, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[testM4, DisplayFunction -> Identity, 
                    PlotStyle -> Red],ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["# H-T",{3.6, -3}],\(Text["# Two Soliton", {3.7, \(-3.5\)}]\),Text["# Manko3 b=",{3.6, -4}],Text[(Subscript[Val, 3]), {4.3, -4}],Text["# Manko4 b=",{3.6, -4.5}],Text[(Subscript[Val, 4]), {4.3, -4.5}],Text["Model with: M=", {2, -2.2}],Text[Mgeo, {2.8, -2.2}],Text["and J=", {2.21, -2.6}],Text[Jgeo, {2.8, -2.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, tt]))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}],plot1=Show[ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 1, 
                        300}], DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],
ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["# H-T",{3.6, -3}],\(Text["# Two Soliton", {3.7, \(-3.5\)}]\),Text["# testM3 b=",{3.6, -4}],Text[(Subscript[Val, 3]), {4.3, -4}],Text["# testM4 b=",{3.6, -4.5}],Text[(Subscript[Val, 4]), {4.3, -4.5}],Text["Model with: M=", {2, -2.2}],Text[Mgeo, {2.8, -2.2}],Text["and J=", {2.21, -2.6}],Text[Jgeo, {2.8, -2.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, tt]))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}]];]
Clear[i];
Print["============= comparing the metric function gt\
φ on the equatorial plane ================"];
Timing[If[Jgeo == N[0],Null,StestVal = 0;test = Table[{0, 0}, {i, 1, 300}];test2 = Table[{0, 0, 0}, {i, 1, 300}];testM3 = Table[{0, 0}, {i, 1, 300}];testM4 = Table[{0, 0}, {i, 1, 300}];testHT = Table[{0, 0, 0}, {i, 1, 300}];For[i = 1, i < 301, 
                test[[i, 1]] = Log[10, Ngtφ[[301 + i, 1]]]; 
                test2[[i, 1]] = Ngtφ[[301 + i, 1]]; 
                testHT[[i, 1]] = 
                  Log[10, Ngtφ[[301 + i, 1]]]; 
                rHT = FindRoot[
                    Ngtφ[[301 + i, 
                          1]] == (((gtφHT[Mgeo, Jgeo/(Mgeo^2), r, 0]^2 - gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] gφφHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0]))^(1/2)), {r, Rs}]; 
                If[NumberQ[
                    (Subscript[Val, 3])], testM3[[i, 1]] = 
                      Log[10, Ngtφ[[301 + i, 1]]];, Null]; 
                If[NumberQ[
                    (Subscript[Val, 4])], testM4[[i, 1]] = 
                      Log[10, Ngtφ[[301 + i, 1]]];, Null]; 
                test[[i, 2]] = 
                  Log[Abs[(Ngtφ[[301 + i, 2]] - 
                              Chop[gtφ[
                                  Ngtφ[[301 + i, 1]], 0]])/
                          Ngtφ[[301 + i, 2]]]]/Log[10]; 
                test2[[i, 2]] = 
                  Chop[gtφ[Ngtφ[[301 + i, 1]], 0]]; 
                test2[[i, 3]] = Ngtφ[[301 + i, 2]]; 
                If[i < 300, 
                  StestVal = 
                    StestVal + (((Ngtφ[[301 + i, 2]] - Chop[gtφ[Ngtφ[[301 + i, 1]], 0]]))^2) (Ngtφ[[301 + i + 1, 1]] - 
                            Ngtφ[[301 + i, 1]]), Null]; 
                testHT[[i, 2]] = 
                  Log[Abs[(Ngtφ[[301 + i, 2]] - 
                                gtφHT[Mgeo, Jgeo/(Mgeo^2), r, 0] //. 
                              rHT)/Ngtφ[[301 + i, 2]]]]/
                    Log[10]; 
                testHT[[i, 3]] = Log[Abs[r //. rHT]]/Log[10]; 
                If[NumberQ[
                    (Subscript[Val, 3])], testM3[[i, 2]] = 
                      Log[10, 
                        Abs[(Ngtφ[[301 + i, 2]] - 
                                Mgtφ3[
                                  Ngtφ[[301 + i, 1]], 0])/
                            Ngtφ[[301 + i, 2]]]];, Null]; 
                If[NumberQ[
                    (Subscript[Val, 4])], testM4[[i, 2]] = 
                      Log[10, 
                        Abs[(Ngtφ[[301 + i, 2]] - 
                                Mgtφ4[
                                  Ngtφ[[301 + i, 1]], 0])/
                            Ngtφ[[301 + i, 2]]]];, Null]; 
                i = i + 1];Off[Graphics::"gptn"];ModelTypeTab[[modIn + 1,25]]=(Sqrt[StestVal]);Print["σ_tφ"];Print[(Sqrt[StestVal])];If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])],plot2=Show[ListPlot[testM3, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> Magenta],
ListPlot[test, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> RGBColor[0, 1, 0]],ListPlot[testM4, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> Red],ListPlot[
                          Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 
                              1, 300}], DisplayFunction -> Identity, 
                          PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                            Log[t]/Log[10]}, {t, (10^(-6)), (10^(-1))}, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["# Manko3",{3.6, -1.5}],Text[(Subscript[Val, 3]), {4.3, -1.5}],Text["# Manko4 b=",{3.6, -2}],\(Text["# Two Soliton", {3.7, \(-2.5\)}]\),Text[(Subscript[Val, 4]), {4.3, -2}],Text["# HT",{3.6, -3}],\(Text["Model with: M=", {2, \(-1.2\)}]\),Text[Mgeo, {2.8, -1.2}],Text["and J=", {2.21, -1.6}],Text[Jgeo, {2.8, -1.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, t])φ)""}\),TextStyle -> {FontFamily -> "Times\
", FontSize -> 14}],plot2=Show[ListPlot[test, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                          Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 
                              1, 300}], DisplayFunction -> Identity, 
                          PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                            Log[t]/Log[10]}, {t, (10^(-6)), (10^(-1))}, 
                          DisplayFunction -> Identity, 
                          PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["# Manko3 b=",{3.6, -1.5}],Text[(Subscript[Val, 3]), {4.3, -1.5}],Text["# Manko4 b=",{3.6, -2}],\(Text["# Two Soliton", {3.7, \(-2.5\)}]\),Text[(Subscript[Val, 4]), {4.3, -2}],Text["# HT",{3.6, -3}],\(Text["Model with: M=", {2, \(-1.2\)}]\),Text[Mgeo, {2.8, -1.2}],Text["and J=", {2.21, -1.6}],Text[Jgeo, {2.8, -1.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, t])φ)""}\),TextStyle -> {FontFamily -> "Times\
", FontSize -> 14}]];];]
Clear[i];
\(Print[\*""============= comparing the metric function Rcirc = (Sqrt[g])φφ on the equatorial plane ============""]\);
Timing[StestVal = 0;test = Table[{0, 0}, {i, 1, 299}];test2 = Table[{0, 0, 0}, {i, 1, 299}];testM3 = Table[{0, 0}, {i, 1, 299}];testM4 = Table[{0, 0}, {i, 1, 299}];testHT = Table[{Null, Null, Null}, {i, 1, 299}];For[i = 1, i < 300, 
          test[[i, 1]] = Log[10, Ngtφ[[301 + i, 1]]]; 
          test2[[i, 1]] = Ngtφ[[301 + i, 1]]; 
          testHT[[i, 1]] = Log[10, Ngtφ[[301 + i, 1]]]; 
          rHT = FindRoot[
              Ngtφ[[301 + i, 
                    1]] == (((gtφHT[Mgeo, Jgeo/(Mgeo^2), r, 0]^2 - gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] gφφHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0]))^(1/2)), {r, Rs}]; 
          test[[i, 2]] = 
            Log[Abs[(√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                              Nfρ[
                                  Ngtφ[[301 + i, 
                                      1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2) - (Sqrt[Chop])[gφφ[Ngtφ[[301 + i, 1]], 0]])/√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - Nfρ[
                              Ngtφ[[301 + i, 
                                  1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2)]]/Log[10]; 
          test2[[i, 
                2]] = (Sqrt[Chop])[gφφ[Ngtφ[[301 + i, 1]], 0]]; 
          test2[[i, 
                3]] = √(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                  Nfρ[
                      Ngtφ[[301 + i, 
                          1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2); 
          If[i < 299, 
            StestVal = 
              StestVal + ((((√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - Nfρ[Ngtφ[[301 + i, 1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2) - (Sqrt[Chop])[gφφ[Ngtφ[[301 + i, 1]], 0]])))^2) (Ngtφ[[301 + i + 1, 1]] - 
                      Ngtφ[[301 + i, 1]]), Null]; 
          testHT[[i, 2]] = 
            Log[Abs[(√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                                Nfρ[
                                    Ngtφ[[301 + i, 
                                        1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2) - 
                          RcircHTNew[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r] //. 
                        rHT)/√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                          Nfρ[
                              Ngtφ[[301 + i, 
                                  1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2)]]/Log[10]; 
          testHT[[i, 3]] = Log[Abs[r //. rHT]]/Log[10]; 
          If[NumberQ[(Subscript[Val, 3])], 
            testM3[[i, 1]] = Log[10, Ngtφ[[301 + i, 1]]]; 
            testM3[[i, 2]] = 
              Log[10, Abs[(√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                              Nfρ[
                                  Ngtφ[[301 + i, 
                                      1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2) - (Sqrt[Chop])[Mgφφ3[Ngtφ[[301 + i, 1]], 0]])/√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - Nfρ[
                              Ngtφ[[301 + i, 
                                  1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2)]];, Null]; 
          If[NumberQ[(Subscript[Val, 4])], 
            testM4[[i, 1]] = Log[10, Ngtφ[[301 + i, 1]]]; 
            testM4[[i, 2]] = 
              Log[10, Abs[(√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - 
                              Nfρ[
                                  Ngtφ[[301 + i, 
                                      1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2) - (Sqrt[Chop])[Mgφφ4[Ngtφ[[301 + i, 1]], 0]])/√(Nfρ[Ngtφ[[301 + i, 1]]]^(-1) Ngtφ[[301 + i, 1]]^2 - Nfρ[
                              Ngtφ[[301 + i, 
                                  1]]] Nωρ[Ngtφ[[301 + i, 1]]]^2)]];, Null]; 
          i = i + 1];Off[Graphics::"gptn"];ModelTypeTab[[modIn + 1,26]]=(Sqrt[StestVal]);Print["σ_φφ"];Print[(Sqrt[StestVal])];If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])],ΔRcircPlot=Show[ListPlot[testM3, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 1]],
ListPlot[testM4, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 0.5]],ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 1, 
                        299}], DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-6)), (10^(-1))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{\(Text["# Manko b=", {3.6, \(-1.5\)}]\),Text[(Subscript[Val, 3]), {4.3, -1.5}],\(Text["# Manko b=", {3.6, \(-2\)}]\),\(Text["# Two Soliton", {3.7, \(-2.5\)}]\),Text["# H-T",{3.7, -3}],Text[(Subscript[Val, 4]), {4.3, -2}],Text["Model with: M=", {2, -1.2}],Text[Mgeo, {2.8, -1.2}],Text["and J=", {2.21, -1.6}],Text[Jgeo, {2.8, -1.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,FrameLabel -> {"log(ρ) (km)", "log(\
𝒟Rcirc)"},TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}],ΔRcircPlot=Show[ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 1, 
                        299}], DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-6)), (10^(-1))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{\(Text["# Manko b=", {3.6, \(-1.5\)}]\),Text["# H-T",{3.7, -3}],Text[(Subscript[Val, 3]), {4.3, -1.5}],\(Text["# Manko b=", {3.6, \(-2\)}]\),\(Text["# Two Soliton", {3.7, \(-2.5\)}]\),Text[(Subscript[Val, 4]), {4.3, -2}],Text["Model with: M=", {2, -1.2}],Text[Mgeo, {2.8, -1.2}],Text["and J=", {2.21, -1.6}],Text[Jgeo, {2.8, -1.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,FrameLabel -> {"log(ρ) (km)", "log(\
𝒟Rcirc)"},TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}]];]
Clear[i];
Print["============ comparing the metric functions gzz-g\
ρρ on the equatorial plane ==============="];
Timing[StestVal = 0;test = Table[{0, 0}, {i, 1, Dimensions[Ngzz][[1]]}];test2 = 
          Table[{0, 0, 0}, {i, 1, Dimensions[Ngzz][[1]]}];testM3 = 
          Table[{0, 0}, {i, 1, Dimensions[Ngzz][[1]]}];testM4 = 
          Table[{0, 0}, {i, 1, Dimensions[Ngzz][[1]]}];testHT = 
          Table[{0, 0, 0}, {i, 1, Dimensions[Ngzz][[1]]}];For[i = 1, i < Dimensions[Ngzz][[1]] + 1, 
          test[[i, 1]] = Log[10, Ngzz[[i, 1]]]; 
          test2[[i, 1]] = Ngzz[[i, 1]]; 
          testHT[[i, 1]] = Log[10, Ngzz[[i, 1]]]; 
          rHT = FindRoot[
              Ngzz[[i, 
                    1]] == (((gtφHT[Mgeo, Jgeo/(Mgeo^2), r, 0]^2 - gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] gφφHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0]))^(1/2)), {r, Rs}]; 
          If[NumberQ[
              (Subscript[Val, 3])], testM3[[i, 1]] = Log[10, Ngzz[[i, 1]]];, 
            Null]; If[
            NumberQ[(Subscript[Val, 4])], testM4[[i, 1]] = 
                Log[10, Ngzz[[i, 1]]];, Null]; 
          test[[i, 2]] = 
            Log[Abs[(Ngzz[[i, 2]] - Chop[gzz[Ngzz[[i, 1]], 0]])/
                    Ngzz[[i, 2]]]]/Log[10]; 
          test2[[i, 2]] = Chop[gzz[Ngzz[[i, 1]], 0]]; 
          test2[[i, 3]] = Ngzz[[i, 2]]; 
          If[i < Dimensions[Ngzz][[1]], 
            StestVal = 
              StestVal + (((Ngzz[[i, 2]] - Chop[gzz[Ngzz[[i, 1]], 0]]))^2) (Ngzz[[i + 1, 1]] - Ngzz[[i, 1]]), Null]; 
          testHT[[i, 2]] = 
            Log[Abs[(Ngzz[[i, 2]] - 
                          gzzHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 0] //. 
                        rHT)/Ngzz[[i, 2]]]]/Log[10]; 
          testHT[[i, 3]] = Log[Abs[r //. rHT]]/Log[10]; 
          If[NumberQ[
              (Subscript[Val, 3])], testM3[[i, 2]] = 
                Log[10, 
                  Abs[(Ngzz[[i, 2]] - Mgzz3[Ngzz[[i, 1]], 0])/
                      Ngzz[[i, 2]]]];, Null]; 
          If[NumberQ[
              (Subscript[Val, 4])], testM4[[i, 2]] = 
                Log[10, 
                  Abs[(Ngzz[[i, 2]] - Mgzz4[Ngzz[[i, 1]], 0])/
                      Ngzz[[i, 2]]]];, Null]; i = i + 1];Off[Graphics::"gptn"];ModelTypeTab[[modIn + 1,27]]=(Sqrt[StestVal]);Print["σ_zz"];Print[(Sqrt[StestVal])];If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])],ΔgzzEqPlot=Show[ListPlot[testM3, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 1]],
ListPlot[testM4, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 0.5]],ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 
                        1, Dimensions[Ngzz][[1]]}], 
                    DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{\(Text["# Manko b=", {3.6, \(-2.5\)}]\),Text[(Subscript[Val, 3]), {4.3, -2.5}],\(Text["# Manko b=", {3.6, \(-3\)}]\),\(Text["# Two Soliton", {3.7, \(-3.5\)}]\),Text["# H-T",{3.7, -4}],Text[(Subscript[Val, 4]), {4.3, -3}],Text["Model with: M=", {2, -2.2}],Text[Mgeo, {2.8, -2.2}],Text["and J=", {2.21, -2.6}],Text[Jgeo, {2.8, -2.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, zz]))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}],ΔgzzEqPlot=Show[ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ListPlot[
                    Table[{testHT[[i, 1]], testHT[[i, 2]]}, {i, 
                        1, Dimensions[Ngzz][[1]]}], 
                    DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ParametricPlot[{Log[NRs]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{\(Text["# Manko b=", {3.6, \(-2.5\)}]\),Text["# H-T",{3.7, -4}],Text[(Subscript[Val, 3]), {4.3, -2.5}],\(Text["# Manko b=", {3.6, \(-3\)}]\),\(Text["# Two Soliton", {3.7, \(-3.5\)}]\),Text[(Subscript[Val, 4]), {4.3, -3}],Text["Model with: M=", {2, -2.2}],Text[Mgeo, {2.8, -2.2}],Text["and J=", {2.21, -2.6}],Text[Jgeo, {2.8, -2.6}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(ρ) (km)", \*""log(𝒟(Subscript[g, zz]))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}]];]
Clear[i];
Print["============== comparing the metric function gtt on \
the axis ================"];
Timing[StestVal = 0;zofrHTtab = 
          Table[{re (( i + 1)/(601 - i)), 
              Chop[NIntegrate[
                  dzdμHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), 
                    re (( i + 1)/(601 - i)), u], {u, 0, 1}]]}, {i, 300, 
              600}];zofrHT = Interpolation[zofrHTtab];test = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}];test2 = 
          Table[{0, 0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}];testM3 = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}];testM4 = Table[{0, 0}, {i, rOfTheta[[MDIV]] + 2, SDIV}];testHT = 
          Table[{Log[10, Nflistz[[i, 1]]], 
              Log[10, Abs[((Nflistz[[i, 2]] + gttHT[Mgeo, Jgeo/(Mgeo^2), M2New/(Mgeo^3), r, 1])/Nflistz)[[i, 2]]]] //. 
                FindRoot[
                  zofrHT[r] == Nflistz[[i, 1]], {r, 10 Rsp}]}, {i, 
              1, Dimensions[Nflistz][[1]]}];For[i = 1, i < Dimensions[zOnAxis][[1]] + 1, 
          rzHT = FindRoot[zofrHT[r] == Nflistz[[i, 1]], {r, Rsp}]; 
          test[[i, 1]] = Log[10, Nflistz[[i, 1]]]; 
          test2[[i, 1]] = Nflistz[[i, 1]]; 
          test[[i, 2]] = 
            Log[10, Abs[(Nflistz[[i, 2]] + 
                      Chop[gtt[0, Nflistz[[i, 1]]]])/
                  Nflistz[[i, 2]]]]; 
          test2[[i, 2]] = Chop[gtt[0, Nflistz[[i, 1]]]]; 
          test2[[i, 3]] = -Nflistz[[i, 2]]; 
          If[i < Dimensions[zOnAxis][[1]], 
            StestVal = 
              StestVal + (((Nflistz[[i, 2]] + Chop[gtt[0, Nflistz[[i, 1]]]]))^2) (Nflistz[[i + 1, 1]] - Nflistz[[i, 1]]), Null]; 
          If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])], 
            testM3[[i, 1]] = Log[10, Nflistz[[i, 1]]]; 
            testM4[[i, 1]] = Log[10, Nflistz[[i, 1]]]; 
            testM3[[i, 2]] = 
              Log[10, Abs[(Nflistz[[i, 2]] + 
                        Mgtt3[0, Nflistz[[i, 1]]])/
                    Nflistz[[i, 2]]]]; 
testM4[[i, 
                  2]] = 
              Log[10, Abs[(Nflistz[[i, 2]] + 
                        Mgtt4[0, Nflistz[[i, 1]]])/
                    Nflistz[[i, 2]]]];, Null]; i = i + 1];ModelTypeTab[[modIn + 1,28]]=(Sqrt[StestVal]);Print["σ_tt Axis"];Print[(Sqrt[StestVal])];If[NumberQ[(Subscript[Val, 3])] && NumberQ[(Subscript[Val, 4])],ΔgttAxPlot=Show[ListPlot[testHT, DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ListPlot[testM3, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 1]],ListPlot[testM4, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 0, 0.5]],ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ParametricPlot[{Log[Rsp]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["\!\(\*
StyleBox[\"#\",\nFontColor->RGBColor[0, 0, 1]]\) Manko b=", {4, -2}],Text[(Subscript[Val, 3]), {4.7, -2}],\(Text["# Manko b=", {4, \(-2.5\)}]\),\(Text["# Two Soliton", {4, \(-3\)}]\),Text[(Subscript[Val, 4]), {4.7, -2.5}],Text["# H-T",{3.7, -4}],Text["Model with: M=", {2, -2}],Text[Mgeo, {2.8, -2}],Text["and   J=", {2.19, -2.4}],Text[Jgeo, {2.8, -2.4}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(z) (km)", \*""log(\!\(𝒟g\_tt\))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}],ΔgttAxPlot=Show[ListPlot[testHT, DisplayFunction -> Identity, 
                    PlotStyle -> Cyan],ListPlot[test, DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[0, 1, 0]],ParametricPlot[{Log[Rsp]/Log[10], 
                      Log[t]/Log[10]}, {t, (10^(-7)), (10^(-2))}, 
                    DisplayFunction -> Identity, 
                    PlotStyle -> RGBColor[1, 0, 0]],Graphics[{Text["\!\(\*
StyleBox[\"#\",\nFontColor->RGBColor[0, 0, 1]]\) Manko b=", {4, -2}],Text[(Subscript[Val, 3]), {4.7, -2}],\(Text["# Manko b=", {4, \(-2.5\)}]\),\(Text["# Two Soliton", {4, \(-3\)}]\),Text["# H-T",{3.7, -4}],Text[(Subscript[Val, 4]), {4.7, -2.5}],Text["Model with: M=", {2, -2}],Text[Mgeo, {2.8, -2}],Text["and   J=", {2.19, -2.4}],Text[Jgeo, {2.8, -2.4}]}],DisplayFunction -> $DisplayFunction,PlotRange -> All,Axes -> False,Frame -> True,\(FrameLabel -> {"log(z) (km)", \*""log(\!\(𝒟g\_tt\))""}\),TextStyle -> {FontFamily -> "Times", 
                      FontSize -> 14}]];]
Clear[i];

Show[
        GraphicsArray[{{plot1, 
              plot2}, {ΔRcircPlot, ΔgzzEqPlot}, {ΔgttAxPlot}}]];

comparing the ISCO for the two-soliton

Print["=============== calculating the ISCO \
=================="];
Clear[
      Dp, Ep, Lp, ρ]; 
    Veff[ρ_] = 
      1 - (1/(ρ^2)) ((Ep^2) gφφ[ρ, 0] + 
              2 Ep Lp gtφ[ρ, 0] + (Lp^2) gtt[ρ, 0]); 
    Clear[Ep, Lp, ρ]; 
    NVeff[ρ_] = 
      1 - (1/(ρ^2)) ((Ep^2) (Nfρ2[ρ]^(-1) (ρ^2) - Nfρ2[ρ] Nωρ2[ρ]^2) + 
              2 Ep Lp (Nfρ2[ρ] Nωρ2[ρ]) + 
              (Lp^2) (-Nfρ2[ρ]));
Print["\!\(\*
StyleBox[\"Co\",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\"-\",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\"rotating\",\nFontColor->RGBColor[1, 0, 0]]\)"];
Off[FindRoot::"lstol"]; 
    sol1 = Chop[
        FindRoot[{Veff[ρ] == 0, 
            D[Veff[ρ], ρ] == 0, 
            D[D[Veff[ρ], ρ], ρ] == 0}, {ρ, 
            NRs}, {Ep, 0.9}, {Lp, 1}]]; Print[sol1]; 
    Nsol1 = FindRoot[{NVeff[ρ] == 0, 
          D[NVeff[ρ], ρ] == 0, 
          D[D[NVeff[ρ], ρ], ρ] == 0}, {ρ, 
          Nflist[[360, 1]]}, {Ep, 0.9}, {Lp, 
          1}];
Print["\!\(\*
StyleBox[\"Counter\",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\"-\",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\"rotating\",\nFontColor->RGBColor[1, 0, 0]]\)"];
sol2 = 
      Chop[FindRoot[{Veff[ρ] == 0, 
            D[Veff[ρ], ρ] == 0, 
            D[D[Veff[ρ], ρ], ρ] == 0}, {ρ, 
            NRs}, {Ep, 1}, {Lp, -3}]]; Print[sol2]; 
    Nsol2 = FindRoot[{NVeff[ρ] == 0, 
          D[NVeff[ρ], ρ] == 0, 
          D[D[NVeff[ρ], ρ], ρ] == 0}, {ρ, 
          Nflist[[360, 1]]}, {Ep, 
          0.9}, {Lp, -1}];
\(Print[
      StyleForm[\*""(Δ((Subscript[R, ISCO])/(Subscript[R, ISCO]))) = "", 
        FontFamily -> "Times", FontSize -> 18, 
        FontWeight -> "Bold", 
        FontColor -> RGBColor[0.3, 0.3, 1]], 
      Chop[\((\((Rs + 
                    hISCOminus - \@Chop[gφφ[ρ, 0]])\)/\((Rs + hISCOminus)\) /. 
              sol2)\)]\ 100\ , "% for the counter-rotating"]; 
    DRiscoTabM[\([modIn, 4]\)] = 
      Chop[\((Rs + 
                hISCOminus - \@Chop[gφφ[ρ, 0]])\)/\((Rs + hISCOminus)\) /. sol2]; DRiscoTabM[\([modIn, 1]\)] = \(-a\)/M; 
    DRiscoTabM[\([modIn, 2]\)] = M2New/Mgeo\^3; 
    DRiscoTabM[\([modIn, 3]\)] = \(-S3New\)/
        Mgeo\^4;\)
\(Off[ReplaceAll::"reps"]; 
    Print[StyleForm[\*""(Δ((Subscript[R, ISCO])/(Subscript[R, ISCO]))) = "", 
        FontFamily -> "Times", FontSize -> 18, 
        FontWeight -> "Bold", 
        FontColor -> RGBColor[0.3, 0.3, 1]], 
      Chop[\((\((Rs + 
                    hISCOplus - \@Chop[gφφ[ρ, 0]])\)/\((Rs + hISCOplus)\) /. sol1)\)]\ 100\ , "% for the co-rotating"]; 
    DRiscoTabP[\([modIn, 4]\)] = 
      Chop[\((Rs + 
                hISCOplus - \@Chop[gφφ[ρ, 0]])\)/\((Rs + hISCOplus)\) /. sol1]; DRiscoTabP[\([modIn, 1]\)] = a/M; 
    DRiscoTabP[\([modIn, 2]\)] = M2New/Mgeo\^3; 
    DRiscoTabP[\([modIn, 3]\)] = S3New/Mgeo\^4;\)
ModelTypeTab[[modIn + 1,15]]=100 Chop[(Rs + 
                    
                    hISCOplus - (Sqrt[Chop])[gφφ[ρ, 0]])/(Rs + hISCOplus) /. sol1];Print[
        StyleForm["j = ", FontFamily -> "Times", 
          FontSize -> 18, FontWeight -> "Bold", 
          FontColor -> RGBColor[1, 0, 0]], a/M];ModelTypeTab[[modIn + 1,6]]=a/M;

\(comparing\ the\ frequencies\ and\ Δ \( Ε\&~\)\)

\(Print[\*""=============== comparing the frequencies and Δ\!\(\(E\&~\)\) ==================""]\);
Timing[Ωisco[ρ_, 
              z_] = (1/D)[
                    gφφ[ρ, z], ρ] (-D[
                      gtφ[ρ, 
                        z], ρ] + √(D[gtφ[ρ, z], ρ]^2 - D[gtt[ρ, z], ρ] D[
                            gφφ[ρ, 
                              z], ρ]));ΩiscoN[ρ_] = (-D[
                      Nfρ[ρ] Nωρ[ρ], ρ] + √(D[Nfρ[ρ] Nωρ[ρ], ρ]^2 - 
                        D[-Nfρ[ρ], ρ] D[
                            Nfρ[ρ]^(-1) (ρ^2) - 
                              Nfρ[ρ] Nωρ[ρ]^2, ρ]))/
              D[Nfρ[ρ]^(-1) (ρ^2) - 
                  Nfρ[ρ] Nωρ[ρ]^2, ρ];Etilda[ρ_] = (-gtt[ρ, 0] - 
                gtφ[ρ, 0] Ωisco[ρ, 
                    0])/(√(-gtt[ρ, 0] - 
                    2 gtφ[ρ, 0] Ωisco[ρ, 
                        0] - gφφ[ρ, 
                        0] (((Ωisco[ρ, 0]))^2)));NEtilda[ρ_] = (Nfρ[ρ] - 
                Nfρ[ρ] Nωρ[ρ] ΩiscoN[ρ])/(√(Nfρ[ρ] - 
                    2 Nfρ[ρ] Nωρ[ρ] ΩiscoN[ρ] - (Nfρ[ρ]^(-1) (ρ^2) - 
                          Nfρ[ρ] Nωρ[ρ]^2) (((ΩiscoN[ρ]))^2)));DEtilda[ρ_] = -Ωisco[ρ, 
                0] D[Etilda[ρ], ρ]/D[Ωisco[ρ, 0], ρ];NDEtilda[ρ_] = -ΩiscoN[ρ] D[NEtilda[ρ], ρ]/D[ΩiscoN[ρ], ρ];AvNDEtildaTab = 
          Table[{NDEtildaTab[[i, 
                  1]], ((NDEtilda[NDEtildaTab[[i - 1, 1]]] + 
                          NDEtilda[
                            NDEtildaTab[[i, 1]]]) (NDEtildaTab[[i, 
                              1]] - 
                          NDEtildaTab[[i - 1, 1]]) + (NDEtilda[
                            NDEtildaTab[[i + 1, 1]]] + 
                          NDEtilda[
                            NDEtildaTab[[i, 
                                1]]]) (NDEtildaTab[[i + 1, 1]] - 
                          NDEtildaTab[[i, 
                              1]]))/(2 (NDEtildaTab[[i + 1, 
                            1]] - NDEtildaTab[[i - 1, 1]]))}, {i, 
              2, Dimensions[NDEtildaTab][[1]] - 1}];AvNDEtilda = Interpolation[AvNDEtildaTab];]
\(Print[
        StyleForm[\*""\!\(Ω\_ISCO\) = "", 
          FontFamily -> "Times", FontSize -> 18, 
          FontWeight -> "Bold", 
          FontColor -> RGBColor[0.3, 0.3, 1]], 
        Chop[Ωisco[ρ, 0] /. 
            sol2], \*"" \!\(km\^\(-1\)\) for the counter-rotating""]\);\(Print[
        StyleForm[\*""\!\(Ω\_ISCO\) = "", 
          FontFamily -> "Times", FontSize -> 18, 
          FontWeight -> "Bold", 
          FontColor -> RGBColor[0.3, 0.3, 1]], 
        Chop[Ωisco[ρ, 0] /. 
            sol1], \*"" \!\(km\^\(-1\)\) for the co-rotating""]\);Print[\(StyleForm[\*""\!\(u\_ISCO\) = "", 
            FontFamily -> "Times", FontSize -> 18, 
            FontWeight -> "Bold", 
            FontColor -> RGBColor[0.3, 0.3, 1]]\),(((Mgeo Chop[Ωisco[ρ, 0] /. sol1]))^(1/3))," for the co-rotating"];Print[
        Re[Ωisco[ρ, 0] /. 
              sol2] 29.979 10,  " kHz"];Print[
        Re[Ωisco[ρ, 0] /. 
              sol1] 29.979 10,  " kHz"];Print[((2 π )/(Re[Ωisco[ρ, 0] /. sol2] 29.979 10000)) Second];Print[((2 π )/(Re[Ωisco[ρ, 0] /. sol1] 29.979 10000)) Second];ModelTypeTab[[modIn + 1,20]]=Re[Ωisco[ρ, 0] /. sol1];ModelTypeTab[[modIn + 1,21]]=(((Mgeo Chop[Ωisco[ρ, 0] /. sol1]))^(1/3));
Print["\!\(\*
StyleBox[\"Ω\",\nFontColor->RGBColor[1, 0, 0]]\)"];
\(Timing[\(Plot[{ΩiscoN[
                10\^x]\ 29.979\ 10, Ωisco[10\^x, 
                0]\ 29.979\ 10}, {x, Log[10, NDEtildaTab[\([1, 1]\)]], 3}, 
          Axes -> False, Frame -> True, 
          PlotStyle -> {Red, Green}, PlotLabel -> "Ω", 
          FrameLabel -> {"log(r) (km)", "Ω(kHz)"}];\)]\)
\(Timing[\(Plot[
          Log[Abs[Chop[\(ΩiscoN[10\^x] - Ωisco[10\^x, 0]\)\/ΩiscoN[10\^x]]]]/Log[10], {x, 
            Log[10, NDEtildaTab[\([1, 1]\)]], 3}, Axes -> False, 
          Frame -> True, PlotStyle -> Red, 
          FrameLabel -> {"log(r) (km)", "ΔΩ"}];\)]\)
\(Timing[\(Plot[{NEtilda[
              10\^x], Etilda[10\^x]}, {x, Log[10, NDEtildaTab[\([1, 1]\)]], 
            3}, Axes -> False, Frame -> True, 
          PlotStyle -> {Red, Green}, 
          PlotLabel -> \*""\!\(\(E\&~\)\)"", 
          FrameLabel -> {"log(r) (km)", \*""\!\(\(E\&~\)\)""}];\)]\)
Timing[Plot[{AvNDEtilda[(10^x)], DEtilda[(10^x)]},{x, Log[10, ρ /. sol2], 3},Axes -> False,Frame -> True,PlotStyle -> {Blue, Green},\(PlotLabel -> \*""Δ\!\(\(E\&~\)\)""\),\(FrameLabel -> {"log(r) (km)", \*""Δ\!\(\(E\&~\)\)""}\)];]
Print["\!\(\*
StyleBox[\(Ω\_ρ\),\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\",\",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[1, 0, 0]]\)\!\(\*
StyleBox[\(Ω\_z\),\nFontColor->RGBColor[1, 0, 0]]\)"];Timing[Ωρ[ρ_] = Ωisco[ρ, 0] - (((-((((gzz[ρ, 0]))^(-1))/2) ((((gtt[ρ, 0] + gtφ[ρ, 0] Ωisco[ρ, 0]))^2) (D[gφφ[ρ, z]/(ρ^2), {ρ, 2}] /. {z -> 0}) - 2 (gtt[ρ, 0] + gtφ[ρ, 0] Ωisco[ρ, 0]) (gtφ[ρ, 0] + gφφ[ρ, 0] Ωisco[ρ, 0]) (D[gtφ[ρ, z]/(ρ^2), {ρ, 2}] /. {z -> 0}) + (((gtφ[ρ, 0] + gφφ[ρ, 0] Ωisco[ρ, 0]))^2) (D[gtt[ρ, z]/(ρ^2), {ρ, 2}] /. {z -> 0}))))^(1/2));
If[Jgeo == N[0],Null,Ωz[ρ_] = Ωisco[ρ, 0] - (((-((((gzz[ρ, 0]))^(-1))/2) ((((gtt[ρ, 0] + gtφ[ρ, 0] Ωisco[ρ, 0]))^2) (D[gφφ[ρ, z]/(ρ^2), {z, 2}] /. {z -> 0}) - 2 (gtt[ρ, 0] + gtφ[ρ, 0] Ωisco[ρ, 0]) (gtφ[ρ, 0] + gφφ[ρ, 0] Ωisco[ρ, 0]) (D[gtφ[ρ, z]/(ρ^2), {z, 2}] /. {z -> 0}) + (((gtφ[ρ, 0] + gφφ[ρ, 0] Ωisco[ρ, 0]))^2) (D[gtt[ρ, z]/(ρ^2), {z, 2}] /. {z -> 0}))))^(1/2));];]
Timing[If[Jgeo == N[0],Null,Plot[{AvNΩμInt[(10^x)] 29.979 10,Ωz[(10^x)] 29.979 10},{x, Log[10, NDEtildaTab[[1, 1]]], 3},Axes -> False,Frame -> True,PlotStyle -> {Blue, Green},PlotLabel -> "Ωz",PlotRange -> All,\(FrameLabel -> {"log(r) (km)", \*""(Subscript[Ω, z])(kHz)""}\)];]]
Timing[Plot[{AvNΩxInt[(10^x)] 29.979 10,Ωρ[(10^x)] 29.979 10},{x, Log[10, NDEtildaTab[[1, 1]]], 3},Axes -> False,Frame -> True,PlotStyle -> {Blue, Green},PlotLabel -> "Ωρ",PlotRange -> All,\(FrameLabel -> {"log(r) (km)", \*""(Subscript[Ω, ρ])(kHz)""}\)];]

ModelTypeTab // MatrixForm
DRiscoTabP // MatrixForm
DRiscoTabM // MatrixForm
ModelTypeTab3 = ModelTypeTab;
DRiscoTabP3 = DRiscoTabP;
DRiscoTabM3 = DRiscoTabM;