Repository: poliastro/poliastro

Issue:
`Orbit.to_ephem` producing incorrect ephemeris time
<!--- Provide a general summary of the issue in the title above -->

The epoch time of the `Ephem` object produced by `Orbit.to_ephem` does not match the epoch of the orbit.

<!--- Provide a detailed description of the problem you have -->

```python
import astropy.units as u
from astropy.time import Time
from poliastro.twobody import Orbit
from poliastro.bodies import Earth

t = Time(123456789, format="unix")
r = [400., 500., 600.] << u.km
v = [1., 2., 8.] << u.km / u.s
orb_obj = Orbit.from_vectors(Earth, r, v, t)
ephem = orb_obj.to_ephem()
print("time", ephem.epochs[0])
print("position", ephem._coordinates[0])
```
This prints out:
```
time 123456861.18573791
position (400., 500., 600.) km
```
The time has been shifted back to the nearest periapsis time, but the initial state is the state at the original orbit epoch.

Relatedly, trying to sample the Ephem at the epoch time throws an interpolation bounds error. Instead, sampling at the orbit epoch time should return the same rv vectors as the orbit.

```python
cartrep = ephem.sample(t)
```

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 cartrep = ephem.sample(t)
      2 print("post sample", cartrep, cartrep.differentials)

File poliastro/ephem.py:349, in Ephem.sample(self, epochs, interpolator)
    346 if epochs is None or epochs.isscalar and (epochs == self.epochs).all():
    347     return self._coordinates
--> 349 coordinates = interpolator.interpolate(
    350     epochs.reshape(-1),
    351     self.epochs,
    352     self._coordinates,
    353 )
    355 return coordinates

File poliastro/ephem.py:106, in SplineInterpolator.interpolate(self, epochs, reference_epochs, coordinates)
    102 xyz_unit = coordinates.xyz.unit
    103 d_xyz_unit = coordinates.differentials["s"].d_xyz.unit
    105 result_xyz = (
--> 106     spline_interp(
    107         coordinates.xyz.value,
    108         reference_epochs.jd,
    109         epochs.jd,
    110         kind=self._kind,
    111     )
    112     << xyz_unit
    113 )
    114 result_d_xyz = (
    115     spline_interp(
    116         coordinates.differentials["s"].d_xyz.value,
   (...)
    121     << d_xyz_unit
    122 )
    124 return CartesianRepresentation(
    125     result_xyz, differentials=CartesianDifferential(result_d_xyz)
    126 )

File poliastro/_math/interpolate.py:9, in spline_interp(y, x, u, kind)
      7 def spline_interp(y, x, u, *, kind="cubic"):
      8     """Interpolates y, sampled at x instants, at u instants using `scipy.interpolate.interp1d`."""
----> 9     y_u = interp1d(x, y, kind=kind)(u)
     10     return y_u

File scipy/interpolate/_polyint.py:78, in _Interpolator1D.__call__(self, x)
     57 """
     58 Evaluate the interpolant
     59 
   (...)
     75 
     76 """
     77 x, x_shape = self._prepare_x(x)
---> 78 y = self._evaluate(x)
     79 return self._finish_y(y, x_shape)

File scipy/interpolate/_interpolate.py:695, in interp1d._evaluate(self, x_new)
    693 y_new = self._call(self, x_new)
    694 if not self._extrapolate:
--> 695     below_bounds, above_bounds = self._check_bounds(x_new)
    696     if len(y_new) > 0:
    697         # Note fill_value must be broadcast up to the proper size
    698         # and flattened to work here
    699         y_new[below_bounds] = self._fill_value_below

File scipy/interpolate/_interpolate.py:724, in interp1d._check_bounds(self, x_new)
    722 # !! Could provide more information about which values are out of bounds
    723 if self.bounds_error and below_bounds.any():
--> 724     raise ValueError("A value in x_new is below the interpolation "
    725                      "range.")
    726 if self.bounds_error and above_bounds.any():
    727     raise ValueError("A value in x_new is above the interpolation "
    728                      "range.")

ValueError: A value in x_new is below the interpolation range.
```

🖥 **Please paste the output of following commands**

<!--- If you are using conda, we need to know some extra information -->

* `conda info -a` (only if you have conda)
* `conda list` (only if you have conda)
* `pip freeze`

<!--- Paste your output of this command after the colon of "Paste your output here:" -->
<details>
  
```bash
# Paste your output here:
...
poliastro==0.17.0
...
```
</details>

💡 **Possible solutions**

I think the problem is [this line](https://github.com/poliastro/poliastro/blob/d85c9bd1e2e8c19765d34e6ee3d672cb51da22cc/src/poliastro/twobody/sampling.py#L150). I think it should be changed to:
```python
epochs = orbit.epoch + (delta_ts - orbit.t_p)
```


Discussion:

Patch:
diff --git a/src/poliastro/twobody/sampling.py b/src/poliastro/twobody/sampling.py
index 194baa55..02a3aa8e 100644
--- a/src/poliastro/twobody/sampling.py
+++ b/src/poliastro/twobody/sampling.py
@@ -147,7 +147,7 @@ class TrueAnomalyBounds(SamplingStrategy):
         delta_ts = (
             np.unwrap(delta_ts, period=orbit.period.to_value(u.s)) << u.s
         )
-        epochs = orbit.epoch + (orbit.t_p + delta_ts)
+        epochs = orbit.epoch + (delta_ts - orbit.t_p)
 
         n = nu_values.shape[0]
         rr, vv = coe2rv_many(


Test patch:
diff --git a/tests/tests_twobody/test_orbit.py b/tests/tests_twobody/test_orbit.py
index b7097d1a..d3f9d20c 100644
--- a/tests/tests_twobody/test_orbit.py
+++ b/tests/tests_twobody/test_orbit.py
@@ -1130,6 +1130,23 @@ def test_from_ephem_has_expected_properties():
     assert_quantity_allclose(ss.v, expected_v)
 
 
+def test_to_ephem_samples_correct_epochs_and_coordinates():
+    epoch = Time("2022-2-3 1:42:0")
+    orbit = Orbit.from_vectors(
+        Earth,
+        r=[6045.0, 3490.0, -2500.0] * u.km,
+        v=[-3.457, 6.618, 2.533] * u.km / u.s,
+        epoch=epoch,
+    )
+
+    ephem = orbit.to_ephem()
+    assert ephem.epochs[0].value == orbit.epoch
+
+    ephem_initial_coordinates = [coord[0] for coord in ephem.sample(epoch).xyz]
+
+    assert_quantity_allclose(ephem_initial_coordinates, orbit.r)
+
+
 def test_from_vectors_wrong_dimensions_fails():
     bad_r = [[1000, 0, 0]] * u.km
     bad_v = [[[0, 10, 0]]] * u.km / u.s