diff --git a/pysatMissions/instruments/missions_sgp4.py b/pysatMissions/instruments/missions_sgp4.py
index 454f7cbb..cb2973cc 100644
--- a/pysatMissions/instruments/missions_sgp4.py
+++ b/pysatMissions/instruments/missions_sgp4.py
@@ -61,6 +61,9 @@ def init(self):
         'August 21–24.'))
     logger.info(self.acknowledgements)
 
+    if 'epoch' not in self.kwargs['load'].keys():
+        self.kwargs['load']['epoch'] = self.files.files.index[0]
+
     return
 
 
@@ -71,7 +74,7 @@ def init(self):
 def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None,
          alt_periapsis=None, alt_apoapsis=None,
          inclination=None, raan=0., arg_periapsis=0., mean_anomaly=0.,
-         bstar=0., one_orbit=False, num_samples=None, cadence='1S'):
+         epoch=None, bstar=0., one_orbit=False, num_samples=None, cadence='1S'):
     """Generate position of satellite in ECI co-ordinates.
 
     Note
@@ -125,6 +128,10 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None,
         orbital plane based on the orbital period.  Optional when instantiating
         via orbital elements.
         (default=0.)
+    epoch : dt.datetime or NoneType
+        The epoch used for calculating orbital propagation from Keplerians.
+        If None, then use the first date in the file list for consistency across
+        multiple days. Note that this will be set at init. (default=None)
     bstar : float
         Inverse of the ballistic coefficient. Used to model satellite drag.
         Measured in inverse distance (1 / earth radius).  Optional when
@@ -176,7 +183,7 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None,
     times, index, dates = ps_meth.generate_times(fnames, num_samples,
                                                  freq=cadence)
     # Calculate epoch for orbital propagator
-    epoch = (dates[0] - dt.datetime(1949, 12, 31)).days
+    epoch_days = (epoch - dt.datetime(1949, 12, 31)).days
     jd, _ = sapi.jday(dates[0].year, dates[0].month, dates[0].day, 0, 0, 0)
 
     if inclination is not None:
@@ -185,7 +192,7 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None,
                                                                 alt_apoapsis)
         satellite = sapi.Satrec()
         # according to module webpage, wgs72 is common
-        satellite.sgp4init(sapi.WGS72, 'i', 0, epoch, bstar, 0, 0,
+        satellite.sgp4init(sapi.WGS72, 'i', 0, epoch_days, bstar, 0, 0,
                            eccentricity, np.radians(arg_periapsis),
                            np.radians(inclination), np.radians(mean_anomaly),
                            mean_motion, np.radians(raan))
diff --git a/pysatMissions/tests/test_instruments.py b/pysatMissions/tests/test_instruments.py
index 1a023207..3894f242 100644
--- a/pysatMissions/tests/test_instruments.py
+++ b/pysatMissions/tests/test_instruments.py
@@ -90,6 +90,44 @@ def teardown_class(self):
 
     # Custom package unit tests can be added here
 
+    @pytest.mark.parametrize("kwargs",
+                             [{},
+                              {'inclination': 20, 'alt_periapsis': 400},
+                              {'inclination': 80, 'alt_periapsis': 500,
+                               'alt_apoapsis': 600,
+                               'epoch': dt.datetime(2019, 1, 1)}])
+    def test_sgp4_data_continuity(self, kwargs):
+        """Test that data is continuous for sequential days.
+
+        Parameters
+        ----------
+        kwargs : dict
+            Optional kwargs to pass through.
+        """
+
+        # Define sat with custom Keplerian inputs
+        sat = pysat.Instrument(
+            inst_module=pysatMissions.instruments.missions_sgp4,
+            **kwargs)
+
+        # Get last 10 points of day 1
+        sat.load(2018, 1)
+        day1 = sat.data[-10:]
+
+        # Get first 10 points of day 2
+        sat.load(2018, 2)
+        day2 = sat.data[:10]
+
+        average_gradient = day1.diff().mean()
+        std_gradient = day1.diff().std()
+        gradient_between_days = day2.iloc[0] - day1.iloc[-1]
+
+        # Check that the jump between days is within 3 sigma of average gradient
+        del_g = np.abs(average_gradient - gradient_between_days)
+        assert np.all(del_g < (3. * std_gradient))
+
+        return
+
     @pytest.mark.parametrize("inst_dict", [x for x in instruments['download']])
     @pytest.mark.parametrize("kwarg,output", [(None, 1), ('10s', 10)])
     def test_inst_cadence(self, inst_dict, kwarg, output):