Clear["Global`*"];
Off[
        General::spell1];

valMessenger = {H -> 2347, R -> 6371 ,  
        GM  -> 398000,  c  -> 300000,  Vf  -> 10.38900, 
        Vinf  -> 4.0560000,  w  -> 1370,  
        di -> 31.44 Pi/180, ai -> 292.61  Pi/180, 
        do -> -31.92 Pi/180, ao -> 227.17  Pi/180, 
        dp -> 46.95  Pi/180, 
ap -> 
          107.5  Pi/180, dcmb -> -7.22 Pi/180, 
        acmb -> 167.99  Pi/180, dpriv -> -40.65 Pi/180, 
        apriv -> 128.75  Pi/180}
valCassini = {H -> 1175, R -> 6371 ,  
        GM  -> 398000,  c  -> 300000,  Vf  -> 19.02600, 
        Vinf  -> 16.010000,  w  -> 1370,  
        di -> -12.92 Pi/180, ai -> 334.31  Pi/180, 
        do -> -4.99 Pi/180, ao -> 352.54  Pi/180, 
        dp -> -23.5 Pi/180, 
ap -> 
          231.4  Pi/180, dcmb -> -7.22 Pi/180, 
        acmb -> 167.99  Pi/180, dpriv -> -40.65 Pi/180, 
        apriv -> 128.75  Pi/180}
valRosetta = {H -> 1956, R -> 6371 ,  
        GM  -> 398000,  c  -> 300000,  Vf  -> 10.51700, 
        Vinf  -> 3.863000,  w  -> 1370,  
        di -> -2.81 Pi/180, ai -> 346.12  Pi/180, 
        do -> -34.29 Pi/180, ao -> 246.51  Pi/180, 
        dp -> 20.20  .8  Pi/180, 
ap -> 
          246.8  Pi/180, dcmb -> -7.22 Pi/180, 
        acmb -> 167.99  Pi/180, dpriv -> -40.65 Pi/180, 
        apriv -> 128.75  Pi/180}
valGalileo2 = {H -> 303, R -> 6371 ,  
        GM  -> 398000,  c  -> 300000,  Vf  -> 14.08000, 
        Vinf  -> 8.877000,  w  -> 1370,  
        di -> -34.26 Pi/180, ai -> 219.35  Pi/180, 
        do -> -4.87 Pi/180, ao -> 174.35  Pi/180, 
        dp -> -33.8 Pi/180, 
ap -> 
          354.4  Pi/180, dcmb -> -7.22 Pi/180, 
        acmb -> 167.99  Pi/180, dpriv -> -40.65 Pi/180, 
        apriv -> 128.75  Pi/180}
valGalileo1 = {H -> 960, R -> 6371 ,  
        GM  -> 398000,  c  -> 300000,  Vf  -> 13.740000, 
        Vinf  -> 8.949000,  w  -> 1370,  e  -> 2.4700000, 
        di -> -12.52 Pi/180, ai -> 266.76  Pi/180, 
        do -> -34.15 Pi/180, ao -> 219.97  Pi/180, 
        dp -> 25.2  Pi/180, 
ap -> 296.5  Pi/180, 
        dcmb -> -7.22 Pi/180, acmb -> 167.99  Pi/180, 
        dpriv -> -40.65 Pi/180, 
        apriv -> 128.75  Pi/180}
valnear = {H -> 539, R -> 6371 ,  GM  -> 398000,  
        c  -> 300000,  Vf  -> 12.739000, 
        Vinf  -> 6.8510000,  w  -> 1370,  
        e  -> 1.8200000, di -> -20.76 Pi/180, 
        ai -> 261.17  Pi/180, do -> -71.96 Pi/180, 
        ao -> 183.49  Pi/180, 
        dp -> 33.  Pi/180, 
ap -> 47.2  Pi/180, 
        dcmb -> -7.22 Pi/180, acmb -> 167.99  Pi/180, 
        dpriv -> -28.93 Pi/180, 
        apriv -> 266.25  Pi/180}

val =  valGalileo1;
init  =  {r[0] == R + H,  r'[0] == 
            0,  θ[0] == 0, θ'[0] == 
            Vf/(R + H)};

ei = 
      N[{Cos [di] Cos [ai], Cos [di] Sin [ai], Sin [di]} /. 
          val]
eo  = 
      N[{Cos[ao] Cos[do], Cos[do] Sin[ao], Sin[do]} /. 
          val]
ecmb  =  
      N[{Cos [dcmb] Cos [acmb], Cos [dcmb] Sin [acmb], Sin [dcmb]} /. 
          val]
epriv = 
      N[{Cos [dpriv] Cos [apriv], Cos [dpriv] Sin [apriv], 
            Sin [dpriv]} /. val]
ep = 
      N[{Cos [dp] Cos [ap], Cos [dp] Sin [ap], Sin [dp]} /. 
          val]

eteta0 = N[Cross[Cross[ei, eo], ep] /. val]
n = N[Cross[ei, eo] /. val]
e  =  
      N[(1 + (Vinf/Vf)^2)/(1 - (Vinf/Vf)^2) /. 
          val]
epriv2  =  {-0.7411267135901881`, 0.12754849748608027`, 
        0.6586555047066535`}
scaltest  =  ep . n
scal =  epriv2 . epriv
ecmb = epriv

eq   =  {  r''[t] - r[t] (θ'[t])^2 +  
              GM/(r[t])^2 (1 - (4/
                          c^2) (Vf ((R + H)/r[t]) 
                            w ecmb . ( 
                                Cos[θ[t]] 
                                    eteta0 + (1 - 
                                        Cos[θ[t]]) (eteta0 . n) 
                                    n + Sin[θ[t]] (Cross[n, 
                                        eteta0])) + 
0 
                             w^2 (ecmb . ( 
                                    Cos[θ[t]] 
                                        ep + (1 - 
                                        Cos[θ[t]]) (ep . n) n + 
                                      Sin[θ[t]] (Cross[n, 
                                        ep])))^2 + 
0 
                             w^2 (ecmb . 
                                  n)^2)
)    ==   0, 
2  r'[t] θ'[t] + 
              r[t] θ''[
                  t] - (Sign[t]) (4/
                    c^2) GM/(r[t])^2 (Vinf w Sqrt[
                        1 + ((1 + e)/(1 - e)) ((R + H)/
                                  r[t])^2 + (2/(e - 1)) ((R + 
                                    H)/r[t])] 
                      ecmb . ( 
                          Cos[θ[t]] 
                              eteta0 + (1 - Cos[θ[t]]) (eteta0 . 
                                  n) n + 
                            Sin[θ[t]] (Cross[n, 
                                  eteta0])) + 
0 
                       w^2 (ecmb . ( 
                            Cos[θ[t]] 
                                eteta0 + (1 - 
                                    Cos[θ[t]]) (eteta0 . n) n + 
                              Sin[θ[t]] (Cross[n, 
                                    eteta0]))) (ecmb . ( 
                            Cos[θ[t]] 
                                ep + (1 - Cos[θ[t]]) (ep . n) 
                                n + Sin[θ[t]] (Cross[n, 
                                    ep])))
) == 0} ;

dsol = NDSolve[
            Join[eq, init] /. 
              val, {r[t], θ[
                t]}, 
                      {t, -200000, 200000}][[1]];

ParametricPlot[
        Evaluate[{r[t] Cos[θ[t]], r[t] Sin [θ[t]]} /. 
            dsol], {t, -20000, 20000}, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        AspectRatio -> Automatic];

Rescor = (N[r[t = -13173] /. dsol, 10] - 
            N[r[t = -13073] /. dsol, 10])/100

Rescor1 = (N[r[t = 13173] /. dsol, 10] - 
            N[r[t = 13073] /. dsol, 10])/100
Resdif  = (Rescor1 - Rescor) 10^6
Resmean = (Rescor + Rescor1)/2

Iang = N[ArcCos[n[[3]]] (180/Pi) /. val]

init0  =  {s[0] == R + H,  s'[0] == 0,  
          T[0] == 0, T'[0] == Vf/(R + H)};

eq0   =  {  s''[x] - s[x] (T'[x])^2 +  
              GM/(s[x])^2    ==   0, 
2  
                    s'[x] T'[x] + s[x] T''[x] == 
            0} ;

dsol0 = NDSolve[
            Join[eq0, init0] /. val, {s[x], 
              T[x]}, 
                      {x, -200000, 200000}][[1]];

ParametricPlot[
        Evaluate[{s[x] Cos[T[x]], s[x] Sin [T[x]]} /. 
            dsol0], {x, -20000, 20000}, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        AspectRatio -> Automatic];

Res0 = (N[s[x = 13173] /. dsol0, 10] - 
            N[s[x = 13073] /. dsol0, 10])/100

(Resmean - Res0) 10^6