Repository: JaiWillems/celest (language: python)

Issue:
Cel 26 update file saving method
The file saving is now approached in a builder style through the implementation of the `TextFileWriter` object. The new implementation allows for multiple layer files and handling of non `Quantity` data. There is also better handling of data column spacings.



Discussion / hints:

Fix patch:
diff --git a/celest/coordinates/__init__.py b/celest/coordinates/__init__.py
index 9039201..9213f0a 100644
--- a/celest/coordinates/__init__.py
+++ b/celest/coordinates/__init__.py
@@ -3,6 +3,7 @@
 
 from celest.coordinates.ground_location import GroundLocation
 
+from celest.coordinates.frames.attitude import Attitude
 from celest.coordinates.frames.azel import AzEl
 from celest.coordinates.frames.gcrs import GCRS
 from celest.coordinates.frames.itrs import ITRS
diff --git a/celest/coordinates/frames/attitude.py b/celest/coordinates/frames/attitude.py
new file mode 100644
index 0000000..365f8ee
--- /dev/null
+++ b/celest/coordinates/frames/attitude.py
@@ -0,0 +1,111 @@
+
+
+from celest.coordinates.frames.base_positions import Position3d
+from celest.coordinates.ground_location import GroundLocation
+from celest.file_save import TextFileWriter
+from celest.units.core import Unit
+from celest.units.quantity import Quantity
+from celest import units as u
+import numpy as np
+
+
+class Attitude(Position3d):
+    """Attitude(julian, roll, pitch, yaw, unit, location)
+
+    Satellite attitude.
+
+    The attitude of a satellite is defined by the roll, pitch, and yaw angles
+    that take the satellite from the lvlh frame to a ground-target-pointing
+    orientation.
+
+    Parameters
+    ----------
+    julian : np.ndarray
+        1-D array containing time in the J2000 epoch.
+    roll, pitch, yaw : np.ndarray
+        1-D array containing the coordinate data.
+    unit : Unit
+        The unit of the angular data.
+    location : GroundLocation
+        The ground location associated with the attitude data.
+
+    Attributes
+    ----------
+    roll, pitch, yaw : Quantity
+        Coordinate data.
+    location : GroundLocation
+        The ground location associated with the attitude data.
+    time : Quantity
+        Times associated with coordinate dimensions.
+
+    Methods
+    -------
+    save_text_file(file_name)
+        Save data as a pretty text file.
+
+    See Also
+    --------
+    AzEl : Azimuth-elevation coordinates.
+    GCRS : Geocentric Celestial Reference System.
+    ITRS : International Terrestrial Reference System.
+    LVLH : Local vertical local horizontal coordinates.
+    WGS84 : Geographical coordinates.
+
+    Examples
+    --------
+    Let `julian`, `roll`, `pitch`, and `yaw` be `np.ndarray` instances and
+    `location` be a `GroundLocation` instance. An `Attitude` object can be
+    initialized:
+
+    >>> attitude = Attitude(julian, roll, picth, yaw, u.deg, location)
+
+    Save `Attitude` data to a text file:
+
+    >>> attitude.save_text_file("attitude_data")
+    """
+
+    def __init__(self, julian: np.ndarray, roll: np.ndarray, pitch: np.ndarray,
+                 yaw: np.ndarray, unit: Unit, location: GroundLocation) -> None:
+
+        super().__init__(roll, unit, pitch, unit, yaw, unit, julian, u.jd2000)
+
+        self._location = location
+
+    @property
+    def roll(self) -> Quantity:
+        return self._get_x()
+
+    @property
+    def pitch(self) -> Quantity:
+        return self._get_y()
+
+    @property
+    def yaw(self) -> Quantity:
+        return self._get_z()
+
+    @property
+    def location(self) -> GroundLocation:
+        return self._location
+
+    def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
+        header = "Attitude Coordinate Data"
+        parameters = [
+            ["Location", str(self._location)]
+        ]
+        data = [
+            ["Time", self.time],
+            ["Roll", self.roll],
+            ["Pitch", self.pitch],
+            ["Yaw", self.yaw]
+        ]
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(parameters=parameters, data=data)
+        writer.save()
diff --git a/celest/coordinates/frames/azel.py b/celest/coordinates/frames/azel.py
index afdb66c..a68003d 100644
--- a/celest/coordinates/frames/azel.py
+++ b/celest/coordinates/frames/azel.py
@@ -2,7 +2,7 @@
 
 from celest.coordinates.frames.base_positions import Position2d
 from celest.coordinates.ground_location import GroundLocation
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.core import Unit
 from celest.units.quantity import Quantity
 from celest import units as u
@@ -10,7 +10,14 @@
 
 
 class AzEl(Position2d):
-    """Coordinates in the horizontal system.
+    """AzEl(julian, azimuth, elevation, unit, location)
+
+    Coordinates in the horizontal system.
+
+    The horizontal system has its origin located at a ground location and has
+    measures of azimith and elevation. Azimuth is defined as the angle of the
+    satellite in the horizontal plane clockwise from north. Elevation is defined
+    as the angle of the satellite above the horzontal plane.
 
     Parameters
     ----------
@@ -39,6 +46,7 @@ class AzEl(Position2d):
 
     See Also
     --------
+    Attitude : Satellite attitude.
     GCRS : Geocentric Celestial Reference System.
     ITRS : International Terrestrial Reference System.
     LVLH : Local vertical local horizontal coordinates.
@@ -66,24 +74,34 @@ def __init__(self, julian: np.ndarray, azimuth: np.ndarray, elevation:
 
     @property
     def azimuth(self) -> Quantity:
-        return self.x
+        return self._get_x()
 
     @property
     def elevation(self) -> Quantity:
-        return self.y
+        return self._get_y()
 
     @property
     def location(self) -> GroundLocation:
         return self._location
 
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         header = "Azimuth-Elevation Coordinate Data"
         parameters = [
             ["Location", str(self._location)]
         ]
         data = [
             ["Time", self.time],
-            ["Azimuth", self.x],
-            ["Elevation", self.y]
+            ["Azimuth", self.azimuth],
+            ["Elevation", self.elevation]
         ]
-        _save_data_as_txt(file_name, header, parameters, data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(parameters=parameters, data=data)
+        writer.save()
diff --git a/celest/coordinates/frames/base_positions.py b/celest/coordinates/frames/base_positions.py
index 55a56ae..70ea120 100644
--- a/celest/coordinates/frames/base_positions.py
+++ b/celest/coordinates/frames/base_positions.py
@@ -30,12 +30,10 @@ def __init__(self, x: np.ndarray, x_unit: Unit, y: np.ndarray, y_unit:
         self._y = Quantity(y, y_unit)
         self._time = Quantity(time, time_unit)
 
-    @property
-    def x(self) -> Quantity:
+    def _get_x(self) -> Quantity:
         return self._x
 
-    @property
-    def y(self) -> Quantity:
+    def _get_y(self) -> Quantity:
         return self._y
 
     @property
@@ -56,8 +54,8 @@ def to_numpy(self, unit: Unit) -> np.ndarray:
             The coordinate data in the specified unit.
         """
         return np.array([
-            self.x.to(unit).data,
-            self.y.to(unit).data
+            self._x.to(unit).data,
+            self._y.to(unit).data
         ]).T
 
 
@@ -88,16 +86,13 @@ def __init__(self, x: np.ndarray, x_unit: Unit, y: np.ndarray, y_unit:
         self._z = Quantity(z, z_unit)
         self._time = Quantity(time, time_unit)
 
-    @property
-    def x(self) -> Quantity:
+    def _get_x(self) -> Quantity:
         return self._x
 
-    @property
-    def y(self) -> Quantity:
+    def _get_y(self) -> Quantity:
         return self._y
 
-    @property
-    def z(self) -> Quantity:
+    def _get_z(self) -> Quantity:
         return self._z
 
     @property
@@ -118,7 +113,7 @@ def to_numpy(self, unit: Unit) -> np.ndarray:
             The coordinate data in the specified unit.
         """
         return np.array([
-            self.x.to(unit).data,
-            self.y.to(unit).data,
-            self.z.to(unit).data
+            self._x.to(unit).data,
+            self._y.to(unit).data,
+            self._z.to(unit).data
         ]).T
diff --git a/celest/coordinates/frames/gcrs.py b/celest/coordinates/frames/gcrs.py
index 0bcd6ac..dce564d 100644
--- a/celest/coordinates/frames/gcrs.py
+++ b/celest/coordinates/frames/gcrs.py
@@ -1,7 +1,7 @@
 
 
 from celest.coordinates.frames.base_positions import Position3d
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.core import Unit
 from celest.units.quantity import Quantity
 from celest import units as u
@@ -9,7 +9,9 @@
 
 
 class GCRS(Position3d):
-    """Coordinates in the Geocentric Celestial Reference System.
+    """GCRS(julian, x, y, z, unit)
+
+    Coordinates in the Geocentric Celestial Reference System.
 
     Parameters
     ----------
@@ -29,14 +31,13 @@ class GCRS(Position3d):
 
     Methods
     -------
-    convert_to(frame, **kwargs)
-        Convert current frame into a new reference frame.
     save_text_file(file_name)
         Save data as a pretty text file.
 
     See Also
     --------
-    AzEl : Azimuth elevation coordinates.
+    Attitude : Satellite attitude.
+    AzEl : Azimuth-elevation coordinates.
     ITRS : International Terrestrial Reference System.
     LVLH : Local vertical local horizontal coordinates.
     WGS84 : Geographical coordinates.
@@ -62,7 +63,27 @@ def __init__(self, julian: np.ndarray, x: np.ndarray, y: np.ndarray, z:
 
         super().__init__(x, unit, y, unit, z, unit, julian, u.jd2000)
 
+    @property
+    def x(self) -> Quantity:
+        return self._get_x()
+
+    @property
+    def y(self) -> Quantity:
+        return self._get_y()
+
+    @property
+    def z(self) -> Quantity:
+        return self._get_z()
+
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         header = "GCRS Coordinate Data"
         data = [
             ["Time", self.time],
@@ -70,4 +91,6 @@ def save_text_file(self, file_name: str) -> None:
             ["Y", self.y],
             ["Z", self.z]
         ]
-        _save_data_as_txt(file_name, header, data=data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
diff --git a/celest/coordinates/frames/itrs.py b/celest/coordinates/frames/itrs.py
index 1aa77f1..5aa91ef 100644
--- a/celest/coordinates/frames/itrs.py
+++ b/celest/coordinates/frames/itrs.py
@@ -1,7 +1,7 @@
 
 
 from celest.coordinates.frames.base_positions import Position3d
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.core import Unit
 from celest.units.quantity import Quantity
 from celest import units as u
@@ -9,7 +9,9 @@
 
 
 class ITRS(Position3d):
-    """Coordinates in the Internation Terrestrial Reference System.
+    """ITRS(julian, x, y, z, unit)
+
+    Coordinates in the Internation Terrestrial Reference System.
 
     Parameters
     ----------
@@ -29,14 +31,13 @@ class ITRS(Position3d):
 
     Methods
     -------
-    convert_to(frame, **kwargs)
-        Convert current frame into a new reference frame.
     save_text_file(file_name)
         Save data as a pretty text file.
 
     See Also
     --------
-    AzEl : Azimuth elevation coordinates.
+    Attitude : Satellite attitude.
+    AzEl : Azimuth-elevation coordinates.
     GCRS : Geocentric Celestial Reference System.
     LVLH : Local vertical local horizontal coordinates.
     WGS84 : Geographical coordinates.
@@ -62,7 +63,27 @@ def __init__(self, julian: np.ndarray, x: np.ndarray, y: np.ndarray, z:
 
         super().__init__(x, unit, y, unit, z, unit, julian, u.jd2000)
 
+    @property
+    def x(self) -> Quantity:
+        return self._get_x()
+
+    @property
+    def y(self) -> Quantity:
+        return self._get_y()
+
+    @property
+    def z(self) -> Quantity:
+        return self._get_z()
+
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         header = "ITRS Coordinate Data"
         data = [
             ["Time", self.time],
@@ -70,4 +91,6 @@ def save_text_file(self, file_name: str) -> None:
             ["Y", self.y],
             ["Z", self.z]
         ]
-        _save_data_as_txt(file_name, header, data=data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
diff --git a/celest/coordinates/frames/lvlh.py b/celest/coordinates/frames/lvlh.py
index faf33ea..9a7bc5d 100644
--- a/celest/coordinates/frames/lvlh.py
+++ b/celest/coordinates/frames/lvlh.py
@@ -2,7 +2,7 @@
 
 from celest.coordinates.frames.base_positions import Position3d
 from celest.coordinates.ground_location import GroundLocation
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.core import Unit
 from celest.units.quantity import Quantity
 from celest import units as u
@@ -10,7 +10,9 @@
 
 
 class LVLH(Position3d):
-    """Coordinates in the level-horizontal-level-vertical or Hill frame.
+    """LVLH(julian, x, y, z, unit)
+
+    Coordinates in the level-horizontal-level-vertical or Hill frame.
 
     The local-vertical local-horizontal frame (also known as the Hill
     frame) is a body frame where the z-axis is algigned with the negative
@@ -42,10 +44,10 @@ class LVLH(Position3d):
 
     See Also
     --------
-    AzEl : Azimuth elevation coordinates.
+    Attitude : Satellite attitude.
+    AzEl : Azimuth-elevation coordinates.
     GCRS : Geocentric Celestial Reference System.
     ITRS : International Terrestrial Reference System.
-    LVLH : Local vertical local horizontal coordinates.
     WGS84 : Geographical coordinates.
 
     Notes
@@ -76,7 +78,27 @@ def __init__(self, julian: np.ndarray, x: np.ndarray, y: np.ndarray,
 
         super().__init__(x, unit, y, unit, z, unit, julian, u.jd2000)
 
+    @property
+    def x(self) -> Quantity:
+        return self._get_x()
+
+    @property
+    def y(self) -> Quantity:
+        return self._get_y()
+
+    @property
+    def z(self) -> Quantity:
+        return self._get_z()
+
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         header = "Satellite LVLH Coordinate Data"
         data = [
             ["Time", self.time],
@@ -84,4 +106,6 @@ def save_text_file(self, file_name: str) -> None:
             ["Y", self.y],
             ["Z", self.z]
         ]
-        _save_data_as_txt(file_name, header, data=data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
diff --git a/celest/coordinates/frames/wgs84.py b/celest/coordinates/frames/wgs84.py
index bba9cd3..dc1a2ee 100644
--- a/celest/coordinates/frames/wgs84.py
+++ b/celest/coordinates/frames/wgs84.py
@@ -1,7 +1,7 @@
 
 
 from celest.coordinates.frames.base_positions import Position3d
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.core import Unit
 from celest.units.quantity import Quantity
 from celest import units as u
@@ -9,7 +9,9 @@
 
 
 class WGS84(Position3d):
-    """Coordinates in the WGS84 Earth ellipsoid model.
+    """WGS84(julian, latitude, longitude, height, angular_unit, length_unit)
+
+    Coordinates in the WGS84 Earth ellipsoid model.
 
     Parameters
     ----------
@@ -31,14 +33,13 @@ class WGS84(Position3d):
 
     Methods
     -------
-    convert_to(frame, **kwargs)
-        Convert current frame into a new reference frame.
     save_text_file(file_name)
         Save data as a pretty text file.
 
     See Also
     --------
-    AzEl : Azimuth elevation coordinates.
+    Attitude : Satellite attitude.
+    AzEl : Azimuth-elevation coordinates.
     GCRS : Geocentric Celestial Reference System.
     ITRS : International Terrestrial Reference System.
     LVLH : Local vertical local horizontal coordinates.
@@ -64,22 +65,32 @@ def __init__(self, julian: np.ndarray, latitude: np.ndarray, longitude:
 
     @property
     def latitude(self) -> Quantity:
-        return self.x
+        return self._get_x()
 
     @property
     def longitude(self) -> Quantity:
-        return self.y
+        return self._get_y()
 
     @property
     def height(self) -> Quantity:
-        return self.z
+        return self._get_z()
 
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         header = "WGS84 Coordinate Data"
         data = [
             ["Time", self.time],
-            ["Latitude", self.x],
-            ["Longitude", self.y],
-            ["Height", self.z]
+            ["Latitude", self.latitude],
+            ["Longitude", self.longitude],
+            ["Height", self.height]
         ]
-        _save_data_as_txt(file_name, header, data=data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
diff --git a/celest/coordinates/ground_location.py b/celest/coordinates/ground_location.py
index a13ed76..89cb62b 100644
--- a/celest/coordinates/ground_location.py
+++ b/celest/coordinates/ground_location.py
@@ -38,6 +38,11 @@ class GroundLocation:
     def __init__(self, latitude: float, longitude: float, height: float,
                  angular_unit: Unit, length_unit: Unit) -> None:
 
+        if latitude < -90 or latitude > 90:
+            raise ValueError("Latitude must be between -90 and 90 degrees.")
+        if longitude < -180 or longitude > 180:
+            raise ValueError("Longitude must be between -180 and 180 degrees.")
+
         self._latitude = Quantity(latitude, angular_unit)
         self._longitude = Quantity(longitude, angular_unit)
         self._height = Quantity(height, length_unit)
diff --git a/celest/coordinates/transforms.py b/celest/coordinates/transforms.py
index b58a393..9667da0 100644
--- a/celest/coordinates/transforms.py
+++ b/celest/coordinates/transforms.py
@@ -74,7 +74,7 @@ def _gcrs_to_itrs(gcrs: GCRS) -> ITRS:
 
 
 def _rotate_by_earth_rotation_angle(data, degree_era):
-    """Rotate vectors about the z-azis by the Earth rotation angle.
+    """Rotate vectors about the z-axis by the Earth rotation angle.
 
     Parameters
     ----------
@@ -375,7 +375,7 @@ def _azimuth(itrs: ITRS, location: GroundLocation) -> np.ndarray:
 
 
 def _get_ang(u, v) -> np.ndarray:
-    """Return the degree angle bewteen two vectors.
+    """Return the degree angle between two vectors.
 
     Parameters
     ----------
diff --git a/celest/encounter/__init__.py b/celest/encounter/__init__.py
index e31b725..262692d 100644
--- a/celest/encounter/__init__.py
+++ b/celest/encounter/__init__.py
@@ -1,3 +1,2 @@
 
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter import windows
+from celest.encounter.window_generator import generate_vtw
diff --git a/celest/encounter/_window_handling.py b/celest/encounter/_window_handling.py
deleted file mode 100644
index e4dd5c1..0000000
--- a/celest/encounter/_window_handling.py
+++ /dev/null
@@ -1,374 +0,0 @@
-
-
-from typing import Any, Literal
-import numpy as np
-import pandas as pd
-
-
-class VTW:
-    """VTW(rise_time, set_time, roll=None, pitch=None, yaw=None)
-
-    Visible window time.
-
-    A visibile window time describes an event where the satellite is in line of
-    sight of the ground station within an elevation constraint.
-
-    Parameters
-    ----------
-    rise_time : float
-        Rise time of the visible time window in julian days.
-    set_time : float
-        Set time of the visible time window in julian days.
-    roll : Stroke, optional
-        Roll angle of the satellite in radians.
-    pitch : Stroke, optional
-        Pitch angle of the satellite in radians.
-    yaw : Stroke, optional
-        Yaw angle of the satellite in radians.
-
-    Attributes
-    ----------
-    rise_time : float
-        Rise time of the visible time window in julian days.
-    set_time : float
-        Set time of the visible time window in julian days.
-    duration : float
-        Duration of the visible time window in seconds.
-    roll : Stroke
-        Roll angle of the satellite in radians.
-    pitch : Stroke
-        Pitch angle of the satellite in radians.
-    yaw : Stroke
-        Yaw angle of the satellite in radians.
-    """
-
-    def __init__(self, rise_time, set_time, roll=None, pitch=None, yaw=None) -> None:
-
-        self.rise_time = rise_time
-        self.set_time = set_time
-        self.duration = 86400 * (set_time - rise_time)
-        self.roll = roll
-        self.pitch = pitch
-        self.yaw = yaw
-
-    def __str__(self) -> str:
-
-        return f"Rise: {self.rise_time}, Set:{self.set_time}, \
-            Duration: {self.duration}s"
-
-
-class OW:
-    """OW(start_time, duration, location, deadline, roll, pitch, yaw)
-
-    Observation window.
-
-    The observation window describes a selected time window in which the
-    satellite-to-ground encounter is to be executed.
-
-    Parameters
-    ----------
-    start_time : float
-        Start time of the observing window in julian days.
-    duration : float
-        Duration of the observing window in seconds.
-    location : GroundPosition
-        Location of the satellite-ground encounter.
-    deadline : float
-        Deadline of the observing window in julian days.
-    roll : float
-        Roll angle of the satellite-ground encounter.
-    pitch : float
-        Pitch angle of the satellite-ground encounter.
-    yaw : float
-        Yaw angle of the satellite-ground encounter.
-
-    Attributes
-    ----------
-    start : float
-        Start time of the observing window in julian days.
-    duration : float
-        Duration of the observing window in seconds.
-    location : GroundPosition
-        Location of the satellite-ground encounter.
-    deadline : float
-        Deadline of the observing window in julian days.
-    roll : float
-        Roll angle of the satellite-ground encounter.
-    pitch : float
-        Pitch angle of the satellite-ground encounter.
-    yaw : float
-        Yaw angle of the satellite-ground encounter.
-    """
-
-    def __init__(self, start_time, duration, location, deadline, roll, pitch, yaw) -> None:
-
-        self.start_time = start_time
-        self.duration = duration
-        self.location = location
-        self.deadline = deadline
-        self.roll = roll
-        self.pitch = pitch
-        self.yaw = yaw
-
-    def __str__(self) -> str:
-
-        return f"Location: {(self.location.latitude, self.location.longitude)}, \
-            Start: {self.start_time}, Duration:{self.duration}s"
-
-
-class WindowHandler:
-    """WindowHandler()
-
-    Data structure to hold `VTW` or `OW` objects.
-    """
-
-    def __init__(self) -> None:
-
-        self.windows = pd.Series(dtype="object")
-        self._index_of_current_window = 0
-        self._length = 0
-
-    def __iter__(self) -> Any:
-
-        return self
-
-    def __next__(self) -> Any:
-
-        if self._index_of_current_window == self._length:
-            self._index_of_current_window = 0
-            raise StopIteration
-        else:
-            current_index = self.windows.index[self._index_of_current_window]
-            self._index_of_current_window += 1
-            return self.windows[current_index]
-
-    def __getitem__(self, key) -> Literal[VTW, OW]:
-        """Return window closest to key time.
-
-        For a scalar index, the window with the closest start time will be
-        returned. For a range index, all windows with start times within the
-        range will be returned. For a tuple index, all the windows with start
-        times closest to each tuple entry will be returned.
-        """
-
-        if isinstance(key, (int, float)):
-            return self.windows[self._get_start_time_of_min_diff(key)]
-        elif isinstance(key, slice):
-            return self.windows[key].to_numpy()
-        elif isinstance(key, tuple):
-            start_times = [self._get_start_time_of_min_diff(k) for k in key]
-            return self.windows[set(start_times)].to_numpy()
-
-    def _get_start_time_of_min_diff(self, window_start_time):
-
-        window_start_times = self.windows.index.to_numpy()
-        min_diff_index = np.abs(window_start_times - window_start_time).argmin()
-        return window_start_times[min_diff_index]
-
-    def __len__(self) -> int:
-
-        return self._length
-
-    def _window_start(self, window) -> float:
-
-        return window.rise_time if isinstance(window, VTW) else window.start_time
-
-    def _add_window_to_window_handler(self, window) -> None:
-
-        partial_window_series = pd.Series(window, index=[self._window_start(window)])
-        self.windows = pd.concat([self.windows, partial_window_series])
-        self.windows.sort_index(inplace=True)
-
-        self._length += 1
-
-    def get_window(self, julian) -> Literal[VTW, OW]:
-        """Return window closest to julian time.
-
-        Parameters
-        ----------
-        julian : float
-            Julian time in julian days.
-
-        Returns
-        -------
-        VTW or OW
-        """
-
-        return self[julian]
-
-    def get_windows_in_range(self, start, end) -> Literal[VTW, OW]:
-        """Return all windows with start times within the range.
-
-        Parameters
-        ----------
-        start : float
-            Start time of the range in julian days.
-        end : float
-            End time of the range in julian days.
-
-        Returns
-        -------
-        VTW or OW
-        """
-
-        return self[start:end]
-
-    def to_numpy(self) -> np.ndarray:
-        """Return windows within a NumPy array.
-
-        Returns
-        -------
-        np.ndarray
-            Array containing window data.
-
-        Examples
-        --------
-        Assuming `toronto_img` is a `Windows` object returned from the windows
-        generation function.
-
-        >>> toronto_img.to_numpy()
-        array([<celest.encounter._window_handling.Window object at 0x0000023FA0D8DFD0>,
-               <celest.encounter._window_handling.Window object at 0x0000023FF0757940>,
-               <celest.encounter._window_handling.Window object at 0x0000023FA10D5E50>,
-               <celest.encounter._window_handling.Window object at 0x0000023FA10D5CA0>],
-               dtype=object)
-        """
-
-        return self.windows.to_numpy()
-
-
-class VTWHandler(WindowHandler):
-    """VTWHandler()
-
-    Data structure to hold `VTW` objects.
-    """
-
-    def __init__(self) -> None:
-
-        super().__init__()
-
-    def _add_window(self, window) -> None:
-
-        if not isinstance(window, VTW):
-            raise TypeError("Expected VTW object.")
-
-        self._add_window_to_window_handler(window)
-
-    def stats(self) -> pd.DataFrame:
-        """Return visible time window statistics.
-
-        This method returns statistics on the visible time windows, including
-        the 5th, 50th, and 95th percentile of encounter durations and
-        enocounter frequency.
-
-        Return
-        ------
-        DataFrame
-            Data base of encounter statistics.
-
-        Examples
-        --------
-        Assuming `toronto_img` is a `VTWHandler` object returned from the
-        windows generation function.
-
-        >>> stats = toronto_img.stats()
-        >>> print(stats)
-                                              Statistic
-        Minimum Duration (s)                        0.0
-        Q5 Duration (s)                      108.468446
-        Q50 Duration (s)                     167.687443
-        Q75 Duration (s)                     177.506299
-        Maximum Duration (s)                      524.0
-        Encounter Frequency (per day)          1.340795
-        """
-
-        if len(self.windows) == 0:
-            return None
-
-        duration = [window.duration for window in self.windows]
-        rise_times = [window.rise_time for window in self.windows]
-        set_times = [window.set_time for window in self.windows]
-
-        data = {}
-        data['Min Duration (s)'] = np.min(duration)
-        data['Q5 Duration (s)'] = np.percentile(duration, 5)
-        data['Q50 Duration (s)'] = np.percentile(duration, 50)
-        data['Q95 Duration (s)'] = np.percentile(duration, 95)
-        data['Max Duration (s)'] = np.max(duration)
-        data['Encounter Frequency (per day)'] = self._length / (np.max(set_times) - np.min(rise_times))
-
-        df = pd.DataFrame(pd.Series(data), columns=["Statistic"])
-
-        return df
-
-    def save(self, fname, delimiter) -> None:
-        """Save visible time window information.
-
-        Parameters
-        ----------
-        fname : str
-            File name to save encounter information under.
-        delimiter : {",", "\\t"}
-            String or character separating columns.
-        """
-
-        heading = [[
-            "Rise Time (JD2000)",
-            "Set Time (JD2000)",
-            "Duration (s)"
-        ]]
-        data = [[
-            w.rise_time,
-            w.set_time,
-            w.duration
-        ] for w in self.to_numpy()]
-        window_data = np.array(heading + data)
-
-        np.savetxt(fname, window_data, delimiter=delimiter, fmt="%s,%s,%s")
-
-
-class OWHandler(WindowHandler):
-    """OWHandler()
-
-    Data structure to hold `OW` objects.
-    """
-
-    def __init__(self) -> None:
-
-        super().__init__()
-
-    def _add_window(self, window) -> None:
-
-        if not isinstance(window, OW):
-            raise TypeError("Expected OW object.")
-
-        self._add_window_to_window_handler(window)
-
-    def save(self, fname, delimiter) -> None:
-        """Save observation window information.
-
-        Parameters
-        ----------
-        fname : str
-            File name to save encounter information under.
-        delimiter : {",", "\\t"}
-            String or character separating columns.
-        """
-
-        heading = [[
-            "Latitude (deg)",
-            "Longitude (deg)",
-            "Start Time (JD2000)",
-            "Duration (s)",
-            "Deadline (JD2000)"
-        ]]
-        data = [[
-            w.location.latitude,
-            w.location.longitude,
-            w.start_time,
-            w.duration,
-            w.deadline
-        ] for w in self.to_numpy()]
-        window_data = np.array(heading + data)
-
-        np.savetxt(fname, window_data, delimiter=delimiter, fmt="%s,%s,%s,%s,%s")
diff --git a/celest/encounter/groundposition.py b/celest/encounter/groundposition.py
deleted file mode 100644
index cc72c60..0000000
--- a/celest/encounter/groundposition.py
+++ /dev/null
@@ -1,76 +0,0 @@
-
-
-import numpy as np
-
-
-WGS84_MAJOR_AXIS = 6378.137
-WGS84_MINOR_AXIS = 6356.752314245
-
-
-class GroundPosition:
-    """GroundPosition(latitude, longitude, height=0)
-
-    Localize Earth surface location information.
-
-    Parameters
-    ----------
-    latitude : float
-        Lattitude of the location in decimal degrees.
-    longitude : float
-        Longitude of the location in decimal degrees.
-    height : float, optional
-        Height of the location above the standard ellipsoid in decimal
-        kilometers.
-
-    Attributes
-    ----------
-    latitude, longitude : float
-        Latitude and longitude of the location in decimal degrees.
-    radius : float
-        Earth radius at (`latitude`, `longitude`).
-    """
-
-    def __init__(self, latitude: float, longitude: float, height: float=0) -> None:
-
-        self.latitude = latitude
-        self.longitude = longitude
-        self.height = height
-        self.radius = self._radius(latitude) + height
-
-    def __str__(self) -> str:
-
-        return f"coor={(self.latitude, self.longitude)}, radius={round(self.radius, 5)}"
-
-    def _radius(self, latitude: float) -> float:
-        """Calculate geocentric radius using WGS84.
-
-        Parameters
-        ----------
-        latitude : float
-            Location's latitude in decimal degrees.
-
-        Returns
-        -------
-        float
-            Earth's geocentric radius in decimal kilometers.
-
-        Notes
-        -----
-        The WGS84 Earth ellipsoid model is used as discussed in "Earth Radius
-        by Latitude (WGS 84)" by Timur. [1]_
-
-        References
-        ----------
-        .. [1] Timur. Earth Radius by Latitude (WGS 84). 2018. url:
-           https://planetcalc.com/7721/.
-        """
-
-        cos_latitude = np.cos(np.radians(latitude))
-        sin_latitude = np.sin(np.radians(latitude))
-
-        numerator = (WGS84_MAJOR_AXIS ** 2 * cos_latitude) ** 2 + \
-            (WGS84_MINOR_AXIS ** 2 * sin_latitude) ** 2
-        denominator = (WGS84_MAJOR_AXIS * cos_latitude) ** 2 + \
-            (WGS84_MINOR_AXIS * sin_latitude) ** 2
-
-        return np.sqrt(numerator / denominator)
diff --git a/celest/encounter/window_generator.py b/celest/encounter/window_generator.py
new file mode 100644
index 0000000..d5d37d9
--- /dev/null
+++ b/celest/encounter/window_generator.py
@@ -0,0 +1,126 @@
+
+
+from celest.coordinates.coordinate import Coordinate
+from celest.coordinates.frames.azel import AzEl
+from celest.coordinates.frames.gcrs import GCRS
+from celest.coordinates.ground_location import GroundLocation
+from celest.encounter.window_handling import VisibleTimeWindow, WindowHandler
+from celest.satellite import Satellite
+from celest.units.quantity import Quantity
+from celest import units as u
+from enum import Enum
+from jplephem.spk import SPK
+import numpy as np
+import pkg_resources
+
+
+class Lighting(Enum):
+    NIGHTTIME = -1
+    ANYTIME = 0
+    DAYTIME = 1
+
+
+def generate_vtw(satellite: Satellite, location: GroundLocation,
+                 vis_threshold: float, lighting:
+                 Lighting=Lighting.ANYTIME) -> WindowHandler:
+    """Return visible time windows for a satellite-ground encounter.
+
+    This function determines the time windows where the satellite has
+    an elevation angle greater than `vis_threshold` and where lighting
+    conditions are met.
+
+    Parameters
+    ----------
+    satellite : Satellite
+    location : GroundLocation
+        Location of ground station.
+    vis_threshold : float
+        Visibility threshold in degrees.
+        
+        The visibility threshold is the minimum elevation angle of the satellite
+        as seen from `location` where the satellite will be in visual range of
+        `location`.
+    lighting : float, optional
+        Lighting condition specifier, defaults to anytime lighting conditions.
+
+        Lighting conditions can be specified using the `Lighting` enum. The
+        different options are `Lighting.NIGHTTIME`, `Lighting.ANYTIME`, and
+        `Lighting.DAYTIME`.
+
+    Returns
+    -------
+    WindowHandler
+        The visible time windows.
+    """
+
+    if satellite.velocity is None:
+        raise Warning("Attitude data cannot be initialized because satellite "
+                      "velocity data is unavailable.")
+    if vis_threshold < 0 or 90 < vis_threshold:
+        raise ValueError("vis_threshold must be between 0 and 90 degrees.")
+
+    satellite_azel = Coordinate(satellite.position).convert_to(AzEl, location)
+    satellite_elevation = satellite_azel.elevation.to(u.deg).data
+    julian = satellite.position.time.to(u.jd2000).data
+
+    window_indices = np.where(satellite_elevation > vis_threshold)[0]
+
+    lighting_indices = _lighting_constraint_indices(julian, location, lighting)
+    window_indices = np.intersect1d(window_indices, lighting_indices)
+
+    window_indices = np.split(window_indices,
+                              np.where(np.diff(window_indices) != 1)[0] + 1)
+    window_indices = np.array(window_indices, dtype=object)
+
+    attitude = satellite.attitude(location)
+
+    windows = WindowHandler()
+    for window_index_set in window_indices:
+        rise_time = julian[window_index_set[0] - 1]
+        set_time = julian[window_index_set[-1]]
+        windows.add_window(VisibleTimeWindow(rise_time, set_time, attitude))
+
+    return windows
+
+
+def _lighting_constraint_indices(julian: np.ndarray, location: GroundLocation,
+                                 lighting: Lighting) -> np.ndarray:
+    if lighting == Lighting.NIGHTTIME:
+        return _get_night_constraint_indices(julian, location)
+    elif lighting == Lighting.ANYTIME:
+        return np.indices(julian.shape)
+    elif lighting == Lighting.DAYTIME:
+        return _get_day_constraint_indices(julian, location)
+    else:
+        raise ValueError("Invalid lighting constraint.")
+
+
+def _get_night_constraint_indices(julian: np.ndarray, location:
+                                  GroundLocation) -> np.ndarray:
+    sun_elevation = _get_sun_elevation(julian, location)
+    return np.where(sun_elevation.to(u.deg).data < 0)[0]
+
+
+def _get_sun_elevation(julian: np.ndarray, location: GroundLocation) -> Quantity:
+    sun_coordinates = Coordinate(_get_sun_gcrs(julian))
+    return sun_coordinates.convert_to(AzEl, location).elevation
+
+
+def _get_sun_gcrs(julian: np.ndarray) -> GCRS:
+    ephem = pkg_resources.resource_filename(__name__, '../data/de421.bsp')
+    kernal = SPK.open(ephem)
+
+    ssb2sunb_pos = kernal[0, 10].compute(julian)
+    ssb2eb_pos = kernal[0, 3].compute(julian)
+    eb2e_pos = kernal[3, 399].compute(julian)
+    e2sun_pos = (ssb2sunb_pos - ssb2eb_pos - eb2e_pos).T
+
+    SPK.close(kernal)
+
+    return GCRS(julian, e2sun_pos[:, 0], e2sun_pos[:, 1], e2sun_pos[:, 2], u.km)
+
+
+def _get_day_constraint_indices(julian: np.ndarray, location:
+                                GroundLocation) -> np.ndarray:
+    sun_elevation = _get_sun_elevation(julian, location)
+    return np.where(sun_elevation.to(u.deg).data > 0)[0]
diff --git a/celest/encounter/window_handling.py b/celest/encounter/window_handling.py
new file mode 100644
index 0000000..d70cfac
--- /dev/null
+++ b/celest/encounter/window_handling.py
@@ -0,0 +1,259 @@
+
+
+from celest.coordinates.frames.attitude import Attitude
+from celest.coordinates.ground_location import GroundLocation
+from celest.file_save import TextFileWriter
+from celest.units.quantity import Quantity
+from celest import units as u
+from typing import Union
+import numpy as np
+
+
+class VisibleTimeWindow:
+    """VisibleTimeWindow(rise_time, set_time, attitude)
+
+    A visible time window is the time at which a satellite is visible from a
+    ground location. It is characterized by the time at which the satellite
+    first becomes visible to when it is last visible.
+
+    Parameters
+    ----------
+    rise_time : float
+        The time at which the satellite is first visible in the jd2000 epoch.
+    set_time : float
+        The time at which the satellite is no longer visible in the jd2000 epoch.
+    attitude : Attitude
+        The attitude of the satellite.
+
+    Attributes
+    ----------
+    rise_time : Quantity
+        The time at which the satellite is visible in the J2000 frame.
+    set_time : Quantity
+        The time at which the satellite is no longer visible in the J2000 frame.
+    attitude : Attitude
+        The attitude of the satellite.
+    """
+
+    def __init__(self, rise_time: float, set_time: float, attitude:
+                 Attitude) -> None:
+        self._rise_time = Quantity(rise_time, u.jd2000)
+        self._set_time = Quantity(set_time, u.jd2000)
+        self._attitude = attitude
+
+    def __str__(self) -> str:
+        return f"Rise time: {self._rise_time}, Set time: {self._set_time}, "\
+               f"Attitude: {self._attitude}"
+
+    def __repr__(self) -> str:
+        return f"VisibleTimeWindow({self._rise_time.data}, "\
+               f"{self._set_time.data}, {self._attitude})"
+
+    @property
+    def rise_time(self) -> Quantity:
+        return self._rise_time
+
+    @property
+    def set_time(self) -> Quantity:
+        return self._set_time
+
+    @property
+    def attitude(self) -> Attitude:
+        return self._attitude
+
+
+class ObservationWindow:
+    """ObservationWindow(start_time, duration, deadline, location, attitude)
+
+    An observation window is a subset of a visible time window and represents
+    the actual encounter to fulfill some mission request.
+
+    Parameters
+    ----------
+    start_time : Quantity
+        The start time of the encounter.
+    duration : Quantity
+        The duration of the encounter.
+    location : GroundLocation
+        The ground location involved in the encounter.
+    attitude : Attitude
+        The attitude of the satellite during the encounter.
+
+    Attributes
+    ----------
+    start_time : Quantity
+        The start time of the encounter in the J2000 frame.
+    duration : Quantity
+        The duration of the encounter in seconds.
+    location : GroundLocation
+        The ground location involved in the encounter.
+    attitude : Attitude
+        The attitude of the satellite during the encounter.
+    """
+
+    def __init__(self, start_time: Quantity, duration: Quantity,
+                 location: GroundLocation, attitude: Attitude) -> None:
+        self._start_time = start_time
+        self._duration = duration
+        self._location = location
+        self._attitude = attitude
+
+    def __str__(self) -> str:
+        return f"Start time: {self._start_time}, Duration: {self._duration}, "\
+            f"Location: {self._location}, Attitude: {self._attitude}"
+
+    def __repr__(self) -> str:
+        return f"ObservationWindow({self._start_time}, {self._duration}, "\
+            f"{self.location}, {self._attitude})"
+
+    @property
+    def start_time(self) -> Quantity:
+        return self._start_time
+
+    @property
+    def duration(self) -> Quantity:
+        return self._duration
+
+    @property
+    def location(self) -> GroundLocation:
+        return self._location
+
+    @property
+    def attitude(self) -> Attitude:
+        return self._attitude
+
+
+class WindowHandler:
+    """WindowHandler()
+
+    Container to hold window data.
+
+    The `WindowHandler` class can only hold one type of window data (i.e.
+    either visible time windows or observation windows) at a time due to there
+    differing attributes and physical significance.
+
+    Methods
+    -------
+    add_window(window)
+        Add a window to the container.
+    save_text_file(filename)
+        Save the window data to a text file.
+
+    Raises
+    ------
+    TypeError
+        If attempting to add both `VisibleTimeWindow` and `ObservationWindow`
+        objects.
+    """
+
+    def __init__(self) -> None:
+        self._window_data = []
+        self._current_window_index = 0
+
+    def __str__(self) -> str:
+        return str(self._window_data)
+
+    def __repr__(self) -> str:
+        return "WindowHandler(" +\
+               ", ".join([repr(i) for i in self._window_data]) + ")"
+
+    def __len__(self) -> int:
+        return len(self._window_data)
+
+    def __iter__(self) -> Union[VisibleTimeWindow, ObservationWindow]:
+        return self
+
+    def __next__(self) -> Union[VisibleTimeWindow, ObservationWindow]:
+        if self._current_window_index < len(self._window_data):
+            self._current_window_index += 1
+            return self._window_data[self._current_window_index - 1]
+        else:
+            self._current_window_index = 0
+            raise StopIteration
+
+    def __getitem__(self, key) -> Union[list, VisibleTimeWindow,
+                                        ObservationWindow]:
+        return self._window_data[key]
+
+    def add_window(self, window: Union[VisibleTimeWindow,
+                   ObservationWindow]) -> None:
+        """Add a window to the list of windows.
+
+        Parameters
+        ----------
+        window : {VisibleTimeWindow, ObservationWindow}
+            The window to add to the list.
+
+            All windows in the window handler must be of the same type.
+        """
+
+        if not isinstance(window, (VisibleTimeWindow, ObservationWindow)):
+            raise TypeError("Window must be of type VisibleTimeWindow or "
+                            "ObservationWindow.")
+        if len(self._window_data) == 0:
+            self._window_data.append(window)
+        elif isinstance(window, self._window_data[0].__class__):
+            self._window_data.append(window)
+        else:
+            raise TypeError("All windows added must be consistent types.")
+
+    def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file to save data to.
+        """
+
+        if len(self._window_data) == 0:
+            raise Exception("No windows to save.")
+
+        header = self._window_data[0].__class__.__name__ + " Data"
+
+        if isinstance(self._window_data[0], VisibleTimeWindow):
+            data = self._get_visible_time_window_data()
+        else:
+            data = self._get_observation_window_data()
+
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
+
+    def _get_visible_time_window_data(self) -> list:
+        rise_time = []
+        set_time = []
+
+        for window in self._window_data:
+            rise_time.append(window.rise_time.data)
+            set_time.append(window.set_time.data)
+
+        return [
+            ["Rise Time", Quantity(np.array(rise_time), u.jd2000)],
+            ["Set Time", Quantity(np.array(set_time), u.jd2000)]
+        ]
+
+    def _get_observation_window_data(self) -> list:
+        start_time = []
+        duration = []
+        location = []
+        roll = []
+        pitch = []
+        yaw = []
+
+        for window in self._window_data:
+            start_time.append(window.start_time.to(u.jd2000).data)
+            duration.append(window.duration.to(u.s).data)
+            location.append(window.location)
+            roll.append(window.attitude.roll.to(u.deg).data[0])
+            pitch.append(window.attitude.pitch.to(u.deg).data[0])
+            yaw.append(window.attitude.yaw.to(u.deg).data[0])
+
+        return [
+            ["Start Time", Quantity(np.array(start_time), u.jd2000)],
+            ["Duration", Quantity(np.array(duration), u.s)],
+            ["Location", np.array(location)],
+            ["Roll", Quantity(np.array(roll), u.deg)],
+            ["Pitch", Quantity(np.array(pitch), u.deg)],
+            ["Yaw", Quantity(np.array(yaw), u.deg)]
+        ]
diff --git a/celest/encounter/windows.py b/celest/encounter/windows.py
deleted file mode 100644
index bb978ed..0000000
--- a/celest/encounter/windows.py
+++ /dev/null
@@ -1,214 +0,0 @@
-
-
-from celest.encounter._window_handling import VTW, VTWHandler
-from celest.satellite.coordinate import Coordinate
-from jplephem.spk import SPK
-import numpy as np
-import pkg_resources
-
-
-def generate_vtw(satellite, location, vis_threshold, lighting=0, tol=1e-5) -> VTWHandler:
-    """Return visible window times for a satellite-ground combination.
-
-    This function generates the visible window times where the satellite has
-    an elevation angle greater than `vis_threshold` and where lighting
-    conditions are met.
-
-    Parameters
-    ----------
-    satellite : Satellite
-        Satellite object.
-    location : Location
-        Location of ground station.
-    vis_threshold : float
-        Visibility threshold in degrees.
-    lighting : float, optional
-        Lighting condition specifier.
-
-        `-1` indicates nighttime encounters, `0` indicates anytime encounters,
-        and `1` indicates daytime encounters.
-    tol : float, optional
-        Window start and end time tolerance in days.
-
-    Returns
-    -------
-    VTWHandler
-        The visible time windows.
-    """
-
-    if vis_threshold < 0 or vis_threshold >= 90:
-        raise ValueError("Visibility threshold shall be in the range [0, 90).")
-    if lighting not in (-1, 0, 1):
-        raise ValueError("Valid lighting conditions include -1, 0, and 1.")
-
-    window_locations = _get_base_vtws(satellite, location, vis_threshold)
-    window_locations *= _lighting_constraint_factor(satellite, location, lighting)
-
-    julian = satellite._julian
-    roll, pitch, yaw = satellite.attitude(location, stroke=True)
-
-    window_indices = _get_window_indices(window_locations, julian)
-    if window_indices.size == 0:
-        return None
-
-    windows = VTWHandler()
-    for window_indices in window_indices:
-
-        rt1, rt2 = julian[window_indices[0] - 1], julian[window_indices[0]]
-        rise_time = _root_find(window_locations, rt1, rt2, tol)
-
-        st1, st2 = julian[window_indices[-1]], julian[window_indices[-1] + 1]
-        set_time = _root_find(window_locations, st1, st2, tol)
-
-        window = VTW(rise_time, set_time, roll, pitch, yaw)
-        windows._add_window(window)
-
-    return windows
-
-
-def _get_base_vtws(satellite, location, vis_threshold):
-
-    elevation = satellite._elevation(location)
-    return _constraint_highlighter(elevation, "ge", vis_threshold)
-
-
-def _constraint_highlighter(decision_var, constraint, ang):
-    """Return raw window highlighter.
-
-    This function applies the constraint between the decision variable and
-    constraint angle to generate a Stroke that is one where the constraint is
-    satisfied and zero elsewhere.
-
-    Parameters
-    ----------
-    decision_var : Stroke
-        Stroke object containing the decision variable.
-    constraint : {"gt", "ge", "lt", "le", "eq"}
-        Constraint type.
-
-        The constraint is applied with the decision variable on the left and
-        the constraint angle on the right. For example, `constraint="gt"`
-        implies `decision_var > ang`.
-    ang : float
-        Constraint angle in units equivalent to `decision_var`.
-
-    Returns
-    -------
-    Stroke
-        Stroke object containing the constraint highlighter.
-    """
-
-    if constraint == "gt":
-        return decision_var > ang
-    elif constraint == "ge":
-        return decision_var >= ang
-    elif constraint == "lt":
-        return decision_var < ang
-    elif constraint == "le":
-        return decision_var <= ang
-    elif constraint == "eq":
-        return decision_var == ang
-    else:
-        raise ValueError("Invalid constraint type.")
-
-
-def _lighting_constraint_factor(satellite, location, lighting):
-
-    if lighting == -1:
-        return _get_night_constraint_factor(satellite._julian, location)
-    elif lighting == 0:
-        return 1
-    elif lighting == 1:
-        return _get_day_constraint_factor(satellite._julian, location)
-
-
-def _get_day_constraint_factor(julian, location):
-
-    return _get_lighting_constraint_factor(julian, location, 1)
-
-
-def _get_night_constraint_factor(julian, location):
-
-    return _get_lighting_constraint_factor(julian, location, -1)
-
-
-def _get_lighting_constraint_factor(julian, location, lighting):
-
-    sun_elevation, _ = _sun_coor(julian).horizontal(location, stroke=True)
-    comparison = "gt" if lighting == 1 else "le"
-    return _constraint_highlighter(sun_elevation, comparison, 0)
-
-
-def _sun_coor(julian) -> Coordinate:
-    """Return the sun's coordinates.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    Coordinate
-        Coordinate object containing the sun's position.
-    """
-
-    ephem = pkg_resources.resource_filename(__name__, '../data/de421.bsp')
-    kernal = SPK.open(ephem)
-
-    ssb2sunb_pos, ssb2sunb_vel = kernal[0, 10].compute_and_differentiate(julian)
-    ssb2eb_pos, ssb2eb_vel = kernal[0, 3].compute_and_differentiate(julian)
-    eb2e_pos, eb2e_vel = kernal[3, 399].compute_and_differentiate(julian)
-
-    e2sun_pos = (ssb2sunb_pos - ssb2eb_pos - eb2e_pos).T
-    e2sun_vel = (ssb2sunb_vel - ssb2eb_vel - eb2e_vel).T
-
-    SPK.close(kernal)
-
-    return Coordinate(e2sun_pos, e2sun_vel, "gcrs", julian)
-
-
-def _get_window_indices(window_locations, julian):
-
-    ind = np.arange(0, julian.size, 1) * window_locations(julian).astype(int)
-    ind = ind[ind != 0]
-    ind = np.split(ind, np.where(np.diff(ind) != 1)[0] + 1)
-
-    return np.array(ind, dtype=object)
-
-
-def _root_find(func, lower_bound, upper_bound, tolerance) -> float:
-    """Return root of f(t) between lg and rg.
-
-    Parameters
-    ----------
-    func : Stroke
-        Stroke object representing the function to be evaluated.
-    tl : float
-        Lower bound of the search interval.
-    upper_bound : float
-        Upper bound of the search interval.
-    tolerance : float
-        Search tolerance.
-
-    Returns
-    -------
-    float
-        Root of `f(t)` between `lg` and `rg`.
-    """
-
-    left_value, right_value = lower_bound, upper_bound
-    left_function_value = int(func(left_value))
-
-    while abs(right_value - left_value) > tolerance:
-
-        center_value = (left_value + right_value) / 2
-        center_function_value = int(func(center_value))
-
-        if abs(left_function_value - center_function_value):
-            right_value = center_value
-        else:
-            left_value = center_value
-            left_function_value = center_function_value
-
-    return (left_value if left_function_value else right_value)
diff --git a/celest/file_save.py b/celest/file_save.py
index f9cdfe5..3032a6b 100644
--- a/celest/file_save.py
+++ b/celest/file_save.py
@@ -1,5 +1,7 @@
 
 
+from celest.units.quantity import Quantity
+from typing import IO, Any, Tuple
 import numpy as np
 
 
@@ -7,104 +9,245 @@
 BLANK_LINE_STRING = "\n\n"
 COLUMN_DELIMITER = "\t"
 
+HEADER_UNDERLINE_CHAR = "="
+COLUMN_UNDERLINE_CHAR = "-"
 
-def _save_data_as_txt(file_name, header=None, parameters=None, data=None, data_format="%.15f"):
-    """Save information as a pretty formatted text file.
+
+class TextFileWriter:
+    """TextFileWriter(file_name, header)
+
+    Build and save data to a text file.
 
     Parameters
     ----------
     file_name : str
-        File name or path string to the desired save location. The ".txt"
-        extension is not part of `file_name`.
+        File name or path to the desired save location.
     header : str
-        Header string to be placed at the top of the file.
-    parameters : List[List[str, str]]
-        List of lists containing a parameter's label-value string pair.
-
-        A label-value pair of `["Label String", "Value String"]` will be saved
-        as "Label String: Value String" after the header.
-    data : List[List[str, Quantity]]
-        List of lists containing a data's label-value pair. The label is a
-        string type and the value must be a `Quantity` object with a NumPy
-        array `data` property.
-
-        A label-value pair of `["Data Label", Quantity(np.array([1, 2]), u.m)]`
-        will have the heading "Data Label (m)".
-    data_format : str
-        Data value format string.
+        Header to be placed at the top of the file.
+
+    Methods
+    -------
+    add_layer(subheading, parameters, data, data_format="%.15f")
+        Add a data layer to the text file.
+    save()
+        Save the text file.
+
+    Examples
+    --------
+    Initialize a `TextFileWriter` object:
+
+    >>> writer = TextFileWriter("file_name.txt", "Header String")
+
+    Add a data layer to the text file:
+
+    >>> writer.add_layer(
+    ...     "Subheading String",
+    ...     [["Label String", "Value String"]],
+    ...     [["Data Label", Quantity(np.array([1, 2]), u.m)]],
+    ...     "%.15f"
+    ... )
+
+    Save the data layers to a text file:
+
+    >>> writer.save()
     """
 
-    processed_file_name = file_name + ".txt"
-    with open(processed_file_name, "w", encoding='utf-8') as file:
+    def __init__(self, file_name: str, header: str=None) -> None:
+        self._file_name = _process_file_name(file_name)
+        self._header = header
+        self._layers = []
+
+    def add_layer(self, subheading: str=None, parameters: list=None,
+                  data: list=None, data_format: str="%.15f") -> None:
+        """Add a data layer to the text file.
+
+        A layer is a collection of information to be saved to a text file
+        including a subheading, parameters, and data.
+
+        Parameters
+        ----------
+        subheading : str
+            Subheading to be placed at the top of the layer.
+        parameters : list
+            List of lists containing a parameter's label-value string pair.
+
+            Each label-value pair of `["Label String", "Value String"]` will be
+            saved as "Label String: Value String" after the subheading.
+        data : list
+            List of lists containing a data's label-value pair. The label is a
+            string type and the value must be a `Quantity` object with a NumPy
+            array `data` property or purely a NumPy array object.
+
+            If the value is a `Quantity` object, the column header will be the
+            name of the data label + " (" + the unit of the data + ")".
+            Otherwise, it will just be the data label.
+        data_format : str
+            Data value format string.
+
+        Examples
+        --------
+        Create the layer data:
+
+        >>> subheading = "Subheading String"
+        >>> parameters = [
+        ...     ["Label 1", "Value 1"],
+        ...     ["Label 2", "Value 2"]
+        ... ]
+        >>> data = [
+        ...     ["Data Label", Quantity(np.array([1, 2]), u.m)],
+        ...     ["Data Label 2", np.array([1, 2])]
+        ... ]
+        >>> data_format = "%.15f"
+
+        Add the layer data to the text file:
+
+        >>> writer = TextFileWriter("file_name.txt", "Header String")
+        >>> writer.add_layer(subheading, parameters, data, data_format)
+        """
+
+        self._layers.append([subheading, parameters, data, data_format])
+
+    def save(self) -> None:
+        """Save stored data to a text file.
+
+        This method creates a text file under the input file name and saves the
+        header and layers to the file.
+        """
 
-        if header is not None:
-            file.write(header)
+        with open(self._file_name, "w", encoding='utf-8') as file:
+            if self._header is not None:
+                _write_header_to_file(file, self._header)
 
-        if (header is not None) and (parameters is not None):
-            file.write(BLANK_LINE_STRING)
+                if len(self._layers) > 0:
+                    file.write(BLANK_LINE_STRING)
 
-        if parameters is not None:
-            for parameter_count, [label, value] in enumerate(parameters):
-                file.write(get_parameter_string(parameter_count,
-                                                len(parameters), label, value))
+            for layer_count, layer_info in enumerate(self._layers):
+                _write_layer_to_file(file, *layer_info)
 
-        if (header is not None or parameters is not None) and (data is not None):
-            file.write(BLANK_LINE_STRING)
+                if layer_count < len(self._layers) - 1:
+                    file.write(BLANK_LINE_STRING)
 
-        if data is not None:
-            column_label_lengths = []
 
-            column_header_string = ""
-            column_underline_string = ""
+def _process_file_name(file_name: str) -> str:
+    return file_name if file_name.endswith(".txt") else file_name + ".txt"
 
-            for data_count, [data_label, data_value] in enumerate(data):
 
-                value_width = get_value_width(data_format, data_value)
-                column_label = data_label + " (" + str(data_value.unit) + ")"
+def _write_header_to_file(file: IO, header: str) -> None:
+    file.write(header + NEWLINE_STRING + HEADER_UNDERLINE_CHAR * len(header))
 
-                column_width = len(column_label)
-                actual_column_width = max(value_width, column_width)
 
-                processed_column_label = column_label.ljust(actual_column_width)
-                column_label_lengths.append(len(processed_column_label))
+def _write_layer_to_file(file: IO, subheading: str, parameters: list,
+                         data: list, data_format: str) -> None:
+    if subheading is not None:
+        _write_subheading_to_file(file, subheading)
+    if (subheading is not None) and (parameters is not None):
+        file.write(BLANK_LINE_STRING)
+    if parameters is not None:
+        _write_parameters_to_file(file, parameters)
+    if (subheading is not None or parameters is not None) and (data is not None):
+        file.write(BLANK_LINE_STRING)
+    if data is not None:
+        _write_data_to_file(file, data, data_format)
+
+
+def _write_subheading_to_file(file: IO, subheading: str) -> None:
+    file.write(subheading)
+
+
+def _write_parameters_to_file(file: IO, parameters: list) -> None:
+    for parameter_count, [label, value] in enumerate(parameters):
+        parameter_string = label + ": " + value
+        if parameter_count < len(parameters) - 1:
+            parameter_string += NEWLINE_STRING
+        file.write(parameter_string)
+
+
+def _write_data_to_file(file: IO, data: list, data_format: str) -> None:
+    column_labels_string = ""
+    column_underline_string = ""
+    column_widths = []
+
+    for data_count, [data_label, data_value] in enumerate(data):
+
+        column_label, column_underline, column_width = _get_column_information(
+            data_label,
+            data_value,
+            data_format
+        )
+
+        column_labels_string += column_label
+        column_underline_string += column_underline
+        column_widths.append(column_width)
+
+        ending = COLUMN_DELIMITER if data_count < len(data) - 1 else NEWLINE_STRING
+        column_labels_string += ending
+        column_underline_string += ending
 
-                column_header_string += processed_column_label
-                column_underline_string += "-" * len(processed_column_label)
+        column_data = _get_and_reshape_column_data(data_value)
+        if data_count == 0:
+            raw_data = column_data
+        else:
+            raw_data = np.concatenate((raw_data, column_data), axis=1)
 
-                if data_count < len(data) - 1:
-                    column_header_string += COLUMN_DELIMITER
-                    column_underline_string += COLUMN_DELIMITER
-                else:
-                    column_header_string += NEWLINE_STRING
-                    column_underline_string += NEWLINE_STRING
+    file.write(column_labels_string)
+    file.write(column_underline_string)
+    _write_column_data_to_file(file, column_widths, raw_data, data_format)
 
-                reshaped_data_array = data_value.data.reshape((-1, 1))
-                if data_count == 0:
-                    raw_data = reshaped_data_array
-                else:
-                    raw_data = np.concatenate((raw_data, reshaped_data_array), axis=1)
 
-            file.write(column_header_string)
-            file.write(column_underline_string)
+def _get_column_information(label: str, value: Any, value_format:
+                            str) -> Tuple[str, str, int]:
+    column_label = label
+    if isinstance(value, Quantity):
+        column_label += " (" + str(value.unit) + ")"
 
-            number_of_rows, number_of_columns = raw_data.shape
-            for row_index in range(number_of_rows):
-                data_row_string = ""
-                for column_index in range(number_of_columns):
-                    data_value_string = data_format % raw_data[row_index][column_index]
-                    data_row_string += data_value_string.ljust(column_label_lengths[column_index])
+    value_width = _get_maximum_value_width(value, value_format)
+    column_width = len(column_label)
+    actual_width = max(value_width, column_width)
 
-                    if column_index < number_of_columns - 1:
-                        data_row_string += COLUMN_DELIMITER
-                if row_index < number_of_rows - 1:
-                    data_row_string += NEWLINE_STRING
+    processed_column_label = column_label.ljust(actual_width)
+    column_underline = COLUMN_UNDERLINE_CHAR * actual_width
 
-                file.write(data_row_string)
+    return processed_column_label, column_underline, actual_width
 
 
-def get_value_width(data_format, data_value):
-    return len(data_format % np.max(data_value.data))
+def _get_maximum_value_width(value: Any, value_format: str) -> int:
+    if isinstance(value, Quantity):
+        return len(value_format % np.max(value.data))
+    elif isinstance(value, np.ndarray):
+        if isinstance(value.dtype, (int, float)):
+            return len(value_format % np.max(value))
+        else:
+            return len(max(value.astype(str), key=len))
+    else:
+        return NotImplemented
 
 
-def get_parameter_string(parameter_count, number_of_parameters, label, value):
-    return label + ": " + value + (NEWLINE_STRING if parameter_count < number_of_parameters - 1 else "")
+def _get_and_reshape_column_data(value: Any) -> np.ndarray:
+    if isinstance(value, Quantity):
+        return value.data.reshape((-1, 1))
+    else:
+        return value.reshape((-1, 1))
+
+
+def _write_column_data_to_file(file: IO, column_widths: list,
+                               raw_data: np.ndarray, data_format: str) -> None:
+    for row_count, row in enumerate(raw_data):
+        row_string = ""
+        for column_count, column in enumerate(row):
+            row_string += _get_data_entry_string(
+                column,
+                column_widths[column_count],
+                data_format
+            )
+            if column_count < len(row) - 1:
+                row_string += COLUMN_DELIMITER
+        if row_count < len(raw_data) - 1:
+            row_string += NEWLINE_STRING
+        file.write(row_string)
+
+
+def _get_data_entry_string(value: Any, width: int, data_format: str) -> str:
+    if isinstance(value, (int, float)):
+        return (data_format % value).ljust(width)
+    else:
+        return str(value).ljust(width)
diff --git a/celest/satellite.py b/celest/satellite.py
index 7229f27..7674746 100644
--- a/celest/satellite.py
+++ b/celest/satellite.py
@@ -1,20 +1,21 @@
 
 
+from celest.coordinates.frames.attitude import Attitude
 from celest.coordinates.frames.gcrs import GCRS
 from celest.coordinates.frames.itrs import ITRS
 from celest.coordinates.ground_location import GroundLocation
 from celest.coordinates.transforms import (
-    _gcrs_to_itrs,
     _altitude,
-    _get_ang,
-    _itrs_to_gcrs,
+    _gcrs_to_itrs,
     _gcrs_to_lvlh,
-    _gcrs_to_lvlh_matrix
+    _gcrs_to_lvlh_matrix,
+    _get_ang,
+    _itrs_to_gcrs
 )
-from celest.file_save import _save_data_as_txt
+from celest.file_save import TextFileWriter
 from celest.units.quantity import Quantity
 from celest import units as u
-from typing import Union, Tuple
+from typing import Tuple, Union
 import numpy as np
 
 
@@ -69,8 +70,7 @@ def __init__(self, position: Union[GCRS, ITRS], velocity:
         else:
             self.velocity = None
 
-    def attitude(self, location: GroundLocation) -> Tuple[Quantity, Quantity,
-                                                          Quantity]:
+    def attitude(self, location: GroundLocation) -> Attitude:
         """Return satellite roll, pitch, and yaw angles.
 
         This method returns the attitude angles required to align the
@@ -115,9 +115,9 @@ def attitude(self, location: GroundLocation) -> Tuple[Quantity, Quantity,
             self.position.time.data,
             np.full((len(lvlh_position.x.data),),
                     location.itrs_x.to(u.km).data),
-            np.full((len(lvlh_position.x.data),),
+            np.full((len(lvlh_position.y.data),),
                     location.itrs_y.to(u.km).data),
-            np.full((len(lvlh_position.x.data),),
+            np.full((len(lvlh_position.z.data),),
                     location.itrs_z.to(u.km).data),
             u.km
         )
@@ -149,11 +149,12 @@ def attitude(self, location: GroundLocation) -> Tuple[Quantity, Quantity,
         a12 = - v[:, 2] + v[:, 0] * v[:, 1] / (1 + c)
         a22 = 1 - (v[:, 0] ** 2 + v[:, 2] ** 2) / (1 + c)
 
-        roll = Quantity(-np.arcsin(a32), u.rad)
-        pitch = Quantity(np.arctan2(a31, a33), u.rad)
-        yaw = Quantity(np.arctan2(a12, a22), u.rad)
+        roll = -np.arcsin(a32)
+        pitch = np.arctan2(a31, a33)
+        yaw = np.arctan2(a12, a22)
 
-        return roll, pitch, yaw
+        return Attitude(self.position.time.data, roll, pitch, yaw, u.rad,
+                        location)
 
     def look_angle(self, location: GroundLocation) -> Quantity:
         """Return look-angles to a ground location.
@@ -233,23 +234,33 @@ def distance(self, location: GroundLocation) -> Quantity:
         return Quantity(distance, u.km)
 
     def save_text_file(self, file_name: str) -> None:
+        """Save data as a pretty text file.
+
+        Parameters
+        ----------
+        file_name : str
+            Name of the text file for the saved data.
+        """
+
         gcrs_position = _itrs_to_gcrs(self.position)
         gcrs_velocity = _itrs_to_gcrs(self.velocity)
 
-        header = "Satellite position and velocity."
+        header = "Satellite Position and Velocity"
         data = [
             ["Time", self.position.time],
-            ["ITRS X", self.position.x],
-            ["ITRS Y", self.position.y],
-            ["ITRS Z", self.position.z],
-            ["ITRS VX", self.velocity.x],
-            ["ITRS VY", self.velocity.y],
-            ["ITRS VZ", self.velocity.z],
-            ["GCRS X", gcrs_position.x],
-            ["GCRS Y", gcrs_position.y],
-            ["GCRS Z", gcrs_position.z],
-            ["GCRS VX", gcrs_velocity.x],
-            ["GCRS VY", gcrs_velocity.y],
-            ["GCRS VZ", gcrs_velocity.z]
+            ["Itrs X", self.position.x],
+            ["Itrs Y", self.position.y],
+            ["Itrs Z", self.position.z],
+            ["Itrs VX", self.velocity.x],
+            ["Itrs VY", self.velocity.y],
+            ["Itrs VZ", self.velocity.z],
+            ["Gcrs X", gcrs_position.x],
+            ["Gcrs Y", gcrs_position.y],
+            ["Gcrs Z", gcrs_position.z],
+            ["Gcrs VX", gcrs_velocity.x],
+            ["Gcrs VY", gcrs_velocity.y],
+            ["Gcrs VZ", gcrs_velocity.z]
         ]
-        _save_data_as_txt(file_name, header, data=data)
+        writer = TextFileWriter(file_name, header)
+        writer.add_layer(data=data)
+        writer.save()
diff --git a/celest/satellite_old/__init__.py b/celest/satellite_old/__init__.py
deleted file mode 100644
index fabbc70..0000000
--- a/celest/satellite_old/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-
-from celest.satellite.time import Time
-from celest.satellite.coordinate import Coordinate
-from celest.satellite.satellite import Satellite
diff --git a/celest/satellite_old/_angle_representations.py b/celest/satellite_old/_angle_representations.py
deleted file mode 100644
index ce9666f..0000000
--- a/celest/satellite_old/_angle_representations.py
+++ /dev/null
@@ -1,102 +0,0 @@
-
-
-import numpy as np
-
-
-DEGREE_SYMBOL = "\u00B0"
-MINUTE_SYMBOL = "\u2032"
-SECOND_SYMBOL = "\u2033"
-
-
-def sexagesimal(angles: np.ndarray) -> np.ndarray:
-    """Convert decimal angles into sexagesimal angles.
-
-    Parameters
-    ----------
-    angles : np.ndarray
-        1-D array containing angles in decimal degrees.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing sexagesimal angles as strings.
-
-    Examples
-    --------
-    >>> angles = np.array([43.6532, -79.3832, -33.2833, 149.1000])
-    >>> coor = Coordinate(...)
-    >>> coor.sexagesimal(angle=angles)
-    np.array(['+43°39′11.52″',
-              '-79°22′59.52″',
-              '-33°16′59.88″',
-              '+149°06′00.00″'])
-    """
-
-    sexagesimal_angles = []
-    for angle in angles:
-        angle_sign = "+" if angle >= 0 else "-"
-        sexagesimal_string = _get_sexagesimal_string(abs(angle))
-        sexagesimal_angles.append(angle_sign + sexagesimal_string)
-
-    return np.array(sexagesimal_angles)
-
-
-def _get_sexagesimal_string(absolute_angle: float) -> str:
-
-    degree_string = "%.3s" % str(int(absolute_angle)).zfill(2)
-    minute_string = "%.2s" % str(int(absolute_angle * 60 % 60)).zfill(2)
-    second_string = str("%.2f" % (absolute_angle * 3600 % 60)).zfill(5)
-
-    return degree_string + DEGREE_SYMBOL + minute_string + MINUTE_SYMBOL + \
-        second_string + SECOND_SYMBOL
-
-
-def _ISO6709_representation(latitude: np.ndarray, longitude: np.ndarray,
-                            altitude: np.ndarray) -> np.ndarray:
-    """Format geographical data in the ISO6709 standard.
-
-    Parameters
-    ----------
-    latitude : np.ndarray
-        1-D array containing latitude angles in decimal degrees.
-    longitude : np.ndarray
-        1-D array containing longitude angles in decimal degrees.
-    altitude : np.ndarray
-        1-D array containing altitude angles in km.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing ISO6709 standardized position strings.
-    """
-
-    iso_strings = []
-    for lat, lon, alt in zip(latitude, longitude, altitude):
-        latitude_string = _get_latitude_string(lat)
-        longitude_string = _get_longitude_string(lon)
-        altitude_string = _get_altitude_string(alt)
-
-        iso_strings.append(latitude_string + " " + longitude_string + " " + altitude_string)
-
-    return np.array(iso_strings)
-
-
-def _get_latitude_string(latitude: float) -> str:
-
-    sexigesimal_latitude = _get_sexagesimal_string(abs(latitude))
-    direction = "N" if latitude >= 0 else "S"
-
-    return sexigesimal_latitude + direction
-
-
-def _get_longitude_string(longitude: float) -> str:
-
-    sexigesimal_longitude = _get_sexagesimal_string(abs(longitude))
-    direction = "E" if longitude >= 0 else "W"
-
-    return sexigesimal_longitude + direction
-
-
-def _get_altitude_string(altitude: float) -> str:
-
-    return ("%.2f" % altitude) + "km"
diff --git a/celest/satellite_old/_astronomical_quantities.py b/celest/satellite_old/_astronomical_quantities.py
deleted file mode 100644
index 38e932e..0000000
--- a/celest/satellite_old/_astronomical_quantities.py
+++ /dev/null
@@ -1,564 +0,0 @@
-
-
-from celest.units.quantity import Quantity
-from typing import Tuple
-import numpy as np
-
-
-def nutation_angles(julian: np.ndarray) -> Tuple:
-    """Return five Earth nutation angles.
-
-    This method calculates the five nutation angles for each time in the
-    `julian` array. The calculated angles include the mean elongation of
-    the Moon from the Sun (D), mean anomaly of the Sun (M), mean anomaly of
-    the Moon (N), Moon's argument of latitude (F), and the longitude of the
-    ascending node of the Moon's mean orbit on the ecliptic measured from the
-    mean equinox date (O).
-
-    Parameters
-    ----------
-    jullian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    tuple
-        Tuple of the form, `(D, M, N, F, O)`, where each item are 1-D NumPy
-        arrays containing the Earth nutation angles in decimal degrees.
-
-    Notes
-    -----
-    The time in Julian centeries since J2000, :math:`T`, can be calculated
-    from the Julian day, :math:`JD`, from the following:
-
-    .. math:: T = \\frac{JD - 2451545}{36525}
-
-    The nutation angles can then be calculated in decimal degrees. [Mee98b]_
-
-    .. math:: D = 297.85036 + 445267.111480T - 0.0019142T^2 + T^3 / 189474
-
-    .. math:: M = 357.52772 + 35999.050340T - 0.0001603T^2 - T^3 / 300000
-
-    .. math:: N = 134.96298 + 477198.867398T + 0.0086972T^2 + T^3 / 56250
-
-    .. math:: F = 93.27191 + 483202.017538T - 0.0036825T^2 + T^3 / 327270
-
-    .. math:: O = 125.04452 - 1934.136261T + 0.0020708T^2 + T^3 / 450000
-
-    References
-    ----------
-    .. [Mee98b] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 143-144. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2446895.5])
-    >>> D, M, N, F, O = nutation_angles(julian=julian)
-    """
-
-    T1, T2, T3 = _calculate_elapsed_JD_century_powers(julian, 3)
-
-    D = _get_mean_elongation_moon_from_sun_deg(T1, T2, T3)
-    M = _get_mean_anomaly_of_sun_deg(T1, T2, T3)
-    N = _get_mean_anomaly_of_moon_deg(T1, T2, T3)
-    F = _get_moons_argument_of_latitude_deg(T1, T2, T3)
-    O = _get_longitude_of_ascending_node_deg(T1, T2, T3)
-
-    return D, M, N, F, O
-
-
-def _calculate_elapsed_JD_century_powers(julian, highest_power):
-
-    elapsd_centuries = (julian - 2451545) / 36525
-    return [elapsd_centuries ** i for i in range(1, highest_power + 1)]
-
-
-def _get_mean_elongation_moon_from_sun_deg(T, T2, T3):
-
-    return (297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474) % 360
-
-
-def _get_mean_anomaly_of_sun_deg(T, T2, T3):
-
-    return (357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000) % 360
-
-
-def _get_mean_anomaly_of_moon_deg(T, T2, T3):
-
-    return (134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250) % 360
-
-
-def _get_moons_argument_of_latitude_deg(T, T2, T3):
-
-    return (93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270) % 360
-
-
-def _get_longitude_of_ascending_node_deg(T, T2, T3):
-
-    return (125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000) % 360
-
-
-def nutation_components(julian: np.ndarray) -> Tuple:
-    """Return the nutations in longitude and obliquity.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    tuple
-        Tuple of the form, `(longitude, obliquity)`, containing 1-D NumPy
-        arrays with the nutation components of longitude and obliquity in
-        decimal seconds.
-
-    Notes
-    -----
-    The methods to calculate the nutations in longitude and obliquity are taken
-    from "Astronomical Algorithms" by Jean Meeus. [Mee98c]_
-
-    References
-    ----------
-    .. [Mee98c] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 143-144. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2446895.5])
-    >>> delta_psi, delta_epsilon = nutation_componenets(julian=julian)
-    """
-
-    T1, T2, T3 = _calculate_elapsed_JD_century_powers(julian, 3)
-
-    ascending_node_longitude = _get_longitude_of_ascending_node_deg(T1, T2, T3)
-    sun_mean_lontitude = _get_mean_longitude_of_sun_deg(T1)
-    moon_mean_longitude = _get_mean_longitude_of_moon_deg(T1)
-
-    ascending_node_longitude = np.radians(ascending_node_longitude)
-    sun_mean_lontitude = np.radians(sun_mean_lontitude)
-    moon_mean_longitude = np.radians(moon_mean_longitude)
-
-    nutation_in_longitude = _get_nutation_in_longitude(ascending_node_longitude,
-                                                       sun_mean_lontitude,
-                                                       moon_mean_longitude)
-    nutation_in_obliquity = _get_nutation_in_obliquity(ascending_node_longitude,
-                                                       sun_mean_lontitude,
-                                                       moon_mean_longitude)
-
-    return nutation_in_longitude, nutation_in_obliquity
-
-
-def _get_mean_longitude_of_sun_deg(T):
-
-    return (280.46645 + 36000.76983 * T)
-
-
-def _get_mean_longitude_of_moon_deg(T):
-
-    return 218.3165 + 481267.8813 * T
-
-
-def _get_nutation_in_longitude(ascending_node_longitude, sun_mean_longitude,
-                               moon_mean_longitude):
-
-    P1 = -17.20 * np.sin(ascending_node_longitude)
-    P2 = -1.32 * np.sin(2 * sun_mean_longitude)
-    P3 = -0.23 * np.sin(2 * moon_mean_longitude)
-    P4 = 0.21 * np.sin(2 * ascending_node_longitude)
-
-    return P1 + P2 + P3 + P4
-
-
-def _get_nutation_in_obliquity(ascending_node_longitude, sun_mean_longitude,
-                               moon_mean_longitude):
-
-    E1 = 9.20 * np.cos(ascending_node_longitude)
-    E2 = 0.57 * np.cos(2 * sun_mean_longitude)
-    E3 = 0.10 * np.cos(2 * moon_mean_longitude)
-    E4 = -0.09 * np.cos(2 * ascending_node_longitude)
-
-    return E1 + E2 + E3 + E4
-
-
-def conventional_precession_angles(julian: np.ndarray) -> Tuple:
-    """Return conventional precession angles.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    Tuple
-        Tuple containing precession angles in decimal arcseconds.
-
-    Notes
-    -----
-    The convetional precession angles are those derived using the IAU 2000A
-    model. [SL13c]_
-
-    References
-    ----------
-    .. [SL13c] M. Soffel and R. Langhans. Space-Time Reference Systems.
-       Astronomy and Astrophysics Library. Springer-Verlag, 2013, pp. 219.
-    """
-
-    T1, T2, T3, T4, T5 = _calculate_elapsed_JD_century_powers(julian, 5)
-
-    zeta = _calculate_zeta_precession_angle_arcseconds(T1, T2, T3, T4, T5)
-    theta = _calculate_theta_precession_angle_arcseconds(T1, T2, T3, T4, T5)
-    z = _calculate_z_precession_angle_arcseconds(T1, T2, T3, T4, T5)
-
-    return zeta, theta, z
-
-
-def _calculate_zeta_precession_angle_arcseconds(T1, T2, T3, T4, T5):
-
-    return 2.59796176 + 2306.0809506 * T1 + 0.3019015 * T2 + 0.0179663 * T3 - \
-        0.0000327 * T4 - 0.0000002 * T5
-
-
-def _calculate_theta_precession_angle_arcseconds(T1, T2, T3, T4, T5):
-
-    return 2004.1917476 * T1 - 0.4269353 * T2 - 0.0418251 * T3 - \
-        0.0000601 * T4 - 0.0000001 * T5
-
-
-def _calculate_z_precession_angle_arcseconds(T1, T2, T3, T4, T5):
-
-    return - 2.5976176 + 2306.0803226 * T1 + 1.0947790 * T2 + 0.0182273 * T3 + \
-        0.0000470 * T4 - 0.0000003 * T5
-
-
-def mean_obliquity(julian: np.ndarray) -> np.ndarray:
-    """Return the ecliptic's mean obliquity in the J2000 epoch.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing the ecliptic's mean obliquity in degrees and
-        decimals.
-
-    See Also
-    --------
-    apparent_obliquity :
-        Return the ecliptic's apparent obliquity in the J2000 epoch.
-
-    Notes
-    -----
-    The methods to calculate the mean obliquity of the ecliptic are given in
-    "Astronomical Algorithms" by Jean Meeus and are only valid for 10000 years
-    on either side of J2000. [Mee98d]_
-
-    References
-    ----------
-    .. [Mee98d] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 147. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2446895.5])
-    >>> epsilon_0 = mean_obliquity(julian=julian)
-    """
-
-    T = _calculate_elapsed_JD_century_powers(julian, 1)[0]
-    U = T / 100
-
-    mean_obliquity = (84381.448 - 4680.93 * U - 1.55 * U ** 2 +
-                      1999.25 * U ** 3 - 51.38 * U ** 4 - 249.67 * U ** 5 -
-                      39.05 * U ** 6 + 7.12 * U ** 7 + 27.87 * U ** 8 +
-                      5.79 * U ** 9 + 2.45 * U ** 10) / 3600
-
-    return mean_obliquity
-
-
-def apparent_obliquity(julian: np.ndarray) -> np.ndarray:
-    """Return the ecliptic's apparent obliquity in the J2000 epoch.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing the ecliptic's apparent obliquity in degrees and
-        decimals.
-
-    See Also
-    --------
-    mean_obliquity :
-        Return the ecliptic's mean obliquity in the J2000 epoch.
-
-    Notes
-    -----
-    The methods to calculate the apparent obliquity of the ecliptic are given
-    in "Astronomical Algorithms" by Jean Meeus and are only valid for 10000
-    years on either side of J2000. [Mee98e]_
-
-    References
-    ----------
-    .. [Mee98e] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 147. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2446895.5])
-    >>> epsilon = apparent_obliquity(julian=julian)
-    """
-
-    average_obliquity = mean_obliquity(julian)
-    _, nutation_in_obliquity = nutation_components(julian)
-
-    return average_obliquity + nutation_in_obliquity / 3600
-
-
-def from_julian(julian: np.ndarray) -> Tuple:
-    """Return year, month, day from Julian times.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    Tuple
-        Returns a tuple, `(year, month, day)`, containing 1-D NumPy arrays
-        representing the year, month, and decimal days of the input times.
-
-    Notes
-    -----
-    Refer to page 63 of Astronomical Algorithms by Jean Meeus for a detailed
-    implementation of the algorithm. [Mee98f]_
-
-    References
-    ----------
-    .. [Mee98f] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 63. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2436116.31])
-    >>> apparent_obliquity(julian=julian)
-    (1957, 10, 4.81)
-    """
-
-    JD = julian + 0.5
-    A, F = JD.astype(int), JD % 1
-
-    ind = np.where(A >= 2291161)[0]
-    alpha = ((A[ind] - 1867216.25) / 36524.25).astype(int)
-    A[ind] = A + 1 + alpha - (alpha / 4).astype(int)
-
-    B = A + 1524
-    C = ((B - 122.1) / 365.25).astype(int)
-    D = (365.25 * C).astype(int)
-    E = ((B - D) / 30.6001).astype(int)
-
-    day = B - D - (30.6001 * E).astype(int) + F
-
-    month = E
-    idx1 = np.where(month < 14)[0]
-    idx2 = np.where((month == 14) | (month == 15))[0]
-    month[idx1] = month[idx1] - 1
-    month[idx2] = month[idx2] - 13
-
-    year = C
-    idx1 = np.where(month > 2)[0]
-    idx2 = np.where((month == 1) | (month == 2))[0]
-    year[idx1] = year[idx1] - 4716
-    year[idx2] = year[idx2] - 4715
-
-    return year, month, day
-
-
-def day_of_year(julian: np.ndarray) -> np.ndarray:
-    """Return the day of the year.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing the day of the year.
-
-    Notes
-    -----
-    The day of year is calculated using the methods in "Astronomical
-    Algorithms" by Jean Meeus. [Mee98g]_
-
-    References
-    ----------
-    .. [Mee98g] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 65. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julData = np.array([2447273.5])
-    >>> AstronomicalQuantities().day_of_year(julData=julData)
-    113
-    """
-
-    year, month, day = from_julian(julian)
-
-    K = np.full(day.shape, 2)
-    K[np.where(year % 4 == 0)] = 1
-
-    N = (275 * month / 9).astype(int)
-    N = N - K * ((month + 9) / 12).astype(int) + day - 30
-
-    return N.astype(int)
-
-
-def equation_of_time(julian: np.ndarray) -> np.ndarray:
-    """Return the equation of time in decimal degrees.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing the Equation of Time in decimal degrees.
-
-    Notes
-    -----
-    The equation of time is calculated using the methods described in
-    "Astronomical Algorithms" by Jean Meeus. [Mee98h]_
-
-    References
-    ----------
-    .. [Mee98h] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 184 - 185. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2448908.5])
-    >>> equation_of_time(julian=julian)
-    3.427351
-    """
-
-    T1, T2, T3 = _calculate_elapsed_JD_century_powers(julian, 3)
-
-    sun_mean_longitude = np.radians(_get_mean_longitude_of_sun_deg(T1))
-    sun_mean_anomaly = np.radians(_get_mean_anomaly_of_sun_deg(T1, T2, T3))
-    earth_eccentricity = _get_earth_eccentricity(T1)
-
-    sun_apparent_obliquity = np.radians(apparent_obliquity(julian))
-    y = np.tan(sun_apparent_obliquity / 2) ** 2
-
-    equation_of_time = y * np.sin(2 * sun_mean_longitude) - \
-        2 * earth_eccentricity * np.sin(sun_mean_anomaly) + \
-        4 * earth_eccentricity * y * np.sin(sun_mean_anomaly) * np.cos(2 * sun_mean_longitude) - \
-        0.5 * y ** 2 * np.sin(4 * sun_mean_longitude) - \
-        1.25 * earth_eccentricity ** 2 * np.sin(2 * sun_mean_anomaly)
-
-    return np.degrees(equation_of_time)
-
-
-def _get_earth_eccentricity(T):
-
-    return 0.016708634 - 0.000042037 * T - 0.0000001267 * T ** 2
-
-
-def equation_of_equinoxes(julian: np.ndarray) -> np.ndarray:
-    """Return the equation of the equinoxes in decimal arcseconds.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing the equation of equinoxes in decimal arcseconds.
-
-    Notes
-    -----
-    The equation of the equinoxes is given by
-    :math:`\frac{1}{15}\Delta\psi\cos\epsilon` where :math:`\Delta\psi` is the
-    nutation in longitude represented in arcseconds and, :math:`\epsilon` is
-    the true obliquity of the ecliptic. [Mee98i]_
-
-    References
-    ----------
-    .. [Mee98i] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 88. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2446895.5]])
-    >>> equation_of_equinoxes(julian=julian)
-    -0.2317
-    """
-
-    nutation_in_longitude, _ = nutation_components(julian)
-    obliquity = np.radians(apparent_obliquity(julian))
-
-    return nutation_in_longitude * np.cos(obliquity) / 15
-
-
-def sun_right_ascension(julian: np.ndarray) -> np.ndarray:
-    """Return the right ascension of the mean sun position.
-
-    This method calculates the Sun's right ascension to an accuracy of 0.01 of
-    a degree.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        1-D array containing containing the right ascension of the mean sun in
-        decimal degrees.
-
-    Notes
-    -----
-    The right ascension of the Sun is calculated using the methods discussed in
-    "Astronomical Algorithms" by Jean Meeus. [Mee98j]_
-
-    References
-    ----------
-    .. [Mee98j] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 163 - 165. isbn: 9780943396613.
-
-    Examples
-    --------
-    >>> julian = np.array([2448908.5])
-    >>> sun_right_ascension(julian=julian)
-    np.array([198.38166])
-    """
-
-    T1, T2, T3 = _calculate_elapsed_JD_century_powers(julian, 3)
-
-    mean_sun_longitude = _get_mean_longitude_of_sun_deg(T1)
-    mean_sun_anomaly = np.radians(_get_mean_anomaly_of_sun_deg(T1, T2, T3))
-
-    sun_center = (1.914602 - 0.004817 * T1 - 0.000014 * T2) * np.sin(mean_sun_anomaly) + \
-        (0.019993 - 0.000101 * T1) * np.sin(2 * mean_sun_anomaly) + \
-        0.000289 * np.sin(3 * mean_sun_anomaly)
-
-    sun_true_longitude = np.radians(mean_sun_longitude + sun_center)
-    obliquity = np.radians(mean_obliquity(julian))
-
-    alpha = np.arctan2(np.cos(obliquity) * np.sin(sun_true_longitude), np.cos(sun_true_longitude))
-
-    return np.degrees(alpha) % 360
diff --git a/celest/satellite_old/coordinate.py b/celest/satellite_old/coordinate.py
deleted file mode 100644
index a215d51..0000000
--- a/celest/satellite_old/coordinate.py
+++ /dev/null
@@ -1,963 +0,0 @@
-
-
-from celest.satellite._angle_representations import _ISO6709_representation
-from celest.satellite.nutation_precession import (
-    bias_matrix,
-    precession_matrix,
-    nutation_matrix
-)
-from celest.satellite.time import Time
-from polare import Stroke
-from typing import Tuple, Union
-import numpy as np
-
-
-WGS84_MAJOR_AXIS = 6378.137
-WGS84_MINOR_AXIS = 6356.752314245
-
-EARTH_ROTATION_RATE = 360.9856123035484
-ERA_AT_JD2000 = 280.46
-
-
-class Coordinate(Time):
-    """Coordinate(position, velocity, frame, julian, offset=0)
-
-    Satellite coordinate transformations.
-
-    The input position, velocity, and time can be converted into different
-    reference frames useful for satellite applications. Here, `julian + offset`
-    is the Julian time in J2000 epoch.
-
-    The `position` and `velocity` inputs require three columns corresponding to
-    the x, y, z cartesian coordinates.
-
-    Parameters
-    ----------
-    position : array_like
-        2-D array containing position coordinates in km. Axis 0 is the time
-        series axis.
-
-        Supported coordinate frames include the Geocentric Celestial Reference
-        System (gcrs) and International Terrestrial Reference System (itrs).
-    velocity : array_like
-        2-D array containing velocity coordinates in m/s. Axis 0 is the time
-        series axis.
-
-        See `position` input for supported coordinate frames.
-    frame : {"gcrs", "itrs"}
-        Frame specifier for the input position.
-    julian : array_like
-        1-D array containing time data in Julian days.
-
-        If times are not in the J2000 epoch, a non-zero offset must be pased
-        in to add to the julian times.
-    offset : float, optional
-        Offset to convert input time data to the J2000 epoch, default is zero.
-
-    Examples
-    --------
-    Initialize `Coordinate` using gcrs data:
-
-    >>> julian = [30462.50000, 30462.50069]
-    >>> position = [[-4681.50824, -5030.09119, 000.00000]
-    ...             [-4714.35352, -4978.74326, 454.41493]]
-    >>> velocity = [[-0.72067, 0.67072, 7.57919]
-    ...             [-0.37378, 1.04024, 7.56237]]
-    >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-    ...                julian=julian, offset=2430000)
-
-    Get horizontal coordinates for the ground location at (43.65, -79.38):
-
-    >>> location = GroundPosition(latitude=43.65, longitude=-79.38, height=0.076)
-    >>> elevation, azimuth = c.horizontal(location=location)
-    """
-
-    def __init__(self, position, velocity, frame, julian, offset=0) -> None:
-
-        super().__init__(julian, offset)
-
-        position = np.array(position)
-        velocity = np.array(velocity)
-        julian = np.array(julian)
-
-        if frame not in ["gcrs", "itrs"]:
-            raise ValueError(f"{frame} is not a valid frame.")
-        if position.ndim != 2:
-            raise ValueError("position input must be 2-D.")
-        if velocity.ndim != 2:
-            raise ValueError("velocity input must be 2-D.")
-        if julian.ndim != 1:
-            raise ValueError("time input must be 1-D.")
-        if position.shape[0] != julian.size:
-            raise ValueError("position and time lengths are mismatched")
-        if velocity.shape[0] != julian.size:
-            raise ValueError("velocity and time lengths are mismatched")
-
-        self._set_interpolataion_kind(julian.size)
-        self._set_base_data(position, velocity, frame)
-
-    def __len__(self) -> int:
-
-        return self._length
-
-    def _set_interpolataion_kind(self, data_length):
-
-        length_to_kind = {1: "zero", 2: "linear", 3: "quadratic"}
-        self._kind = "cubic" if data_length > 3 else length_to_kind[data_length]
-
-    def _set_base_data(self, position, velocity, frame) -> None:
-        """Initialize base data.
-
-        The `position` and `velocity` inputs require three columns
-        corresponding to the x, y, z cartesian coordinates.
-
-        Parameters
-        ----------
-        position : array_like
-            2-D array containing position coordinates.
-
-            Supported coordinate frames include the Geocentric Celestial
-            Reference System (gcrs) and International Terrestrial Reference
-            System (itrs).
-        velocity : array_like
-            2-D array containing velocity coordinates.
-
-            See `position` input for supported coordinate frames.
-        frame : {"gcrs", "itrs"}
-            Frame specifier for the input position.
-        """
-
-        self._length = len(position)
-
-        base_coor = np.concatenate((position, velocity), axis=1)
-
-        converted_position = self._gcrs_and_itrs(position, frame)
-        converted_velocity = self._gcrs_and_itrs(velocity, frame)
-        converted_coor = np.concatenate((converted_position, converted_velocity),
-                                        axis=1)
-
-        if frame == "gcrs":
-            self._GCRS = self._stroke_array_from_2D_array_cols(base_coor)
-            self._ITRS = self._stroke_array_from_2D_array_cols(converted_coor)
-        elif frame == "itrs":
-            self._GCRS = self._stroke_array_from_2D_array_cols(converted_coor)
-            self._ITRS = self._stroke_array_from_2D_array_cols(base_coor)
-
-    def _stroke_array_from_2D_array_cols(self, data):
-
-        return np.array([self._stroke_init(i) for i in data.T])
-
-    def _stroke_init(self, data):
-
-        return Stroke(self._julian, data, self._kind)
-
-    def geo(self, iso=False, stroke=False) -> Union[Tuple, np.ndarray]:
-        """Return geographical coordinates.
-
-        The method cannot return data as both a stroke and iso
-        formatted strings. If both `iso=True` and `stroke=True`, the Stroke
-        return will take precedence.
-
-        Parameters
-        ----------
-        iso : bool, optional
-            Formats output as ISO6709 sexagesimal position strings if true.
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        Union[Tuple, np.ndarray]
-            Tuple of latitude, longitude, and altitude given in decimal
-            degrees and km.
-
-            If `iso=True`, the output is a 1-D array containing formatted
-            strings.
-
-            If `stroke=True`, the output is a tuple of stroke objects.
-
-        See Also
-        --------
-        itrs : Return itrs coordinates.
-        gcrs : Return gcrs coordinates.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Get geographical data:
-
-        >>> c.geo()
-        (array([-0.09540, 3.69543]),
-         array([-22.60391, -23.35431]),
-         array([493.42306, 493.59707]))
-
-        Return formatted output strings by setting `iso=True`:
-
-        >>> c.geo(iso=True)
-        array(['00°05′37.86″S 22°36′05.44″W 493.42km'
-               '03°42′10.85″N 23°21′06.73″W 493.60km'])
-        """
-
-        latitude, longitude = self._itrs_to_geo(self._ITRS[:3])
-        altitude = self.altitude(True)
-
-        if stroke:
-            return latitude, longitude, altitude
-
-        latitude = latitude(self._julian)
-        longitude = longitude(self._julian)
-        altitude = altitude(self._julian)
-
-        if iso:
-            return _ISO6709_representation(latitude, longitude, altitude)
-        else:
-            return latitude, longitude, altitude
-
-    def _itrs_to_geo(self, position) -> np.ndarray:
-        """Itrs to geographical transformation.
-
-        Parameters
-        ----------
-        position : np.ndarray
-            2-D array with columns of x, y, z cartesian coordinates in ITRS.
-
-        Returns
-        -------
-        np.ndarray
-            2-D array with columns of latitude, longitude, and altitude given
-            in decimal degrees and kilometres.
-
-        See Also
-        --------
-        _geo_to_itrs : Geographical to itrs transformation.
-        """
-
-        x, y, z = position
-        latitude = np.arctan(z / np.sqrt(x ** 2 + y ** 2))
-        longitude = np.arctan2(y, x)
-
-        return np.degrees(latitude), np.degrees(longitude)
-
-    def _geo_to_itrs(self, position) -> np.ndarray:
-        """Geographical to itrs transformation.
-
-        Parameters
-        ----------
-        position : np.ndarray
-            2-D array with columns of geodetic latitude, terrestrial longitude,
-            and geodetic altitude given in decimal degrees and kilometres.
-
-        Returns
-        -------
-        np.ndarray
-            2-D array with columns of x, y, z cartesian coordinates in ITRS.
-
-        See Also
-        --------
-        _itrs_to_geo : Itrs to geographical transformation.
-
-        Notes
-        -----
-        An Earth ellipsoid model is used for the geographical to itrs
-        conversion using the methods described in "Coordinate Systems in
-        Geodesy" by E. J. Krakiwsky and D.E. Wells as presented by Christopher
-        Lum. [KW98a]_ [Lum20]_
-
-        References
-        ----------
-        .. [KW98a] E. J. Krakiwsky and D. E. Wells. Coordinate Systems in
-           Geodesy. Jan. 1998.
-        .. [Lum20] Christopher Lum. Geodetic Coordinates: Computing Latitude and
-           Longitude. June 2020.url:https://www.youtube.com/watch?v=4BJ-GpYbZlU.
-        """
-
-        latitude, longitude, height = position.T
-        latitude, longitude = np.radians(latitude), np.radians(longitude)
-
-        e = np.sqrt(1 - WGS84_MINOR_AXIS ** 2 / WGS84_MAJOR_AXIS ** 2)
-        n = WGS84_MAJOR_AXIS / np.sqrt(1 - e ** 2 * np.sin(latitude) ** 2)
-
-        a = n + height
-        b = np.cos(latitude)
-        c = np.cos(longitude)
-        x = (n + height) * np.cos(latitude) * np.cos(longitude)
-        y = (n + height) * np.cos(latitude) * np.sin(longitude)
-        z = (n * (1 - e ** 2) + height) * np.sin(latitude)
-
-        return np.concatenate([i.reshape((-1, 1)) for i in [x, y, z]], axis=1)
-
-    def era(self, stroke=False) -> Union[np.ndarray, Stroke]:
-        """Return Earth rotation angle in decimal degrees.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a stroke objects if true. Otherwise, the output
-            is a 1-D array.
-
-        Returns
-        -------
-        Union[np.ndarray, Stroke]
-            Earth rotation angle in decimal degrees.
-
-            If `stroke=True`, the output is a stroke object. Otherwise, the
-            output is a 1-D array.
-
-        Notes
-        -----
-        The Earth rotation angle is calculated using the following:
-
-        .. math:: \gamma^\circ = 360.9856123035484\Delta T + 280.46
-
-        where :math:`\Delta T=JD-2451545` is the elapsed days since the J2000
-        epoch where :math:`JD` is the Julian day. [Kok17a]_
-
-        References
-        ----------
-        .. [Kok17a] Don Koks. Changing Coordinates in the Context of Orbital
-           Mechanics. Cyber and Electronic Warfare Division, Defence Science,
-           and Technology Group, Jan.2017, p. 12 - 13.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Get era data:
-
-        >>> c.era()
-        array([249.65772, 249.90840])
-        """
-
-        time_since_JD2000 = self._julian - 2451545
-        rotation_angle = (EARTH_ROTATION_RATE * time_since_JD2000 +
-                          ERA_AT_JD2000) % 360
-
-        return self._stroke_init(rotation_angle) if stroke else rotation_angle
-
-    def _gcrs_and_itrs(self, data, frame) -> np.ndarray:
-        """Transform between gcrs and itrs coordinates.
-
-        Parameters
-        ----------
-        data : np.ndarray
-            2-D array containing x, y, z cartesian coordinates.
-        frame : {"itrs", "gcrs"}
-            Input data frame.
-
-        Returns
-        -------
-        np.ndarray
-            2-D array containing the transformed x, y, z cartesian data.
-
-        Notes
-        -----
-        The conversion between the gcrs and itrs frames accounts for Earth
-        rotation, precession, and nutation. Polar motion effects are ignored
-        due to their poor predictability.
-        """
-
-        return self._gcrs_to_itrs(data) if frame == "gcrs" else self._itrs_to_gcrs(data)
-
-    def _gcrs_to_itrs(self, gcrs) -> np.ndarray:
-
-        gcrs = np.copy(gcrs)
-
-        earth_rotation_angle = self.era()
-        bias_mtrx = bias_matrix()
-        precession_mtrx = precession_matrix(self._julian)
-        nutation_mtrx = nutation_matrix(self._julian)
-
-        gcrs = np.einsum('ij, kj -> ki', bias_mtrx, gcrs)
-        gcrs = np.einsum('ijk, ik -> ij', precession_mtrx, gcrs)
-        gcrs = np.einsum('ijk, ik -> ij', nutation_mtrx, gcrs)
-
-        gcrs = self._rotate_by_earth_rotation_angle(gcrs, -earth_rotation_angle)
-
-        return gcrs
-
-    def _rotate_by_earth_rotation_angle(self, data, degree_era):
-
-        radian_era = np.radians(degree_era)
-        c, s = np.cos(radian_era), np.sin(radian_era)
-        c1, c2, c3 = np.copy(data).T
-
-        data[:, 0] = c * c1 - s * c2
-        data[:, 1] = s * c1 + c * c2
-        data[:, 2] = c3
-
-        return data
-
-    def _itrs_to_gcrs(self, itrs) -> np.ndarray:
-
-        itrs = np.copy(itrs)
-
-        earth_rotation_angle = self.era()
-        bias_mtrx = bias_matrix()
-        precession_mtrx = precession_matrix(self._julian)
-        nutation_mtrx = nutation_matrix(self._julian)
-
-        itrs = self._rotate_by_earth_rotation_angle(itrs, earth_rotation_angle)
-
-        itrs = np.einsum('ijk, ik -> ij', np.linalg.inv(nutation_mtrx), itrs)
-        itrs = np.einsum('ijk, ik -> ij', np.linalg.inv(precession_mtrx), itrs)
-        itrs = np.einsum('ij, kj -> ki', np.linalg.inv(bias_mtrx), itrs)
-
-        return itrs
-
-    def gcrs(self, stroke=False) -> Tuple:
-        """Return gcrs coordinates.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        Tuple
-            Tuple containing x, y, z cartesian positions and velocities in the
-            gcrs frame.
-
-            The coordinates are stroke objects if `stroke=True`. Otherwise,
-            1-D arrays are returned.
-
-        See Also
-        --------
-        itrs : Return itrs coordinates.
-        geo : Return geographical coordinates.
-
-        Examples
-        --------
-        Initialize `Coordinate` using itrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[6343.70553, -2641.13728, -11.44111]
-        ...             [6295.54131, -2718.36512, 442.89640]]
-        >>> velocity = [[-0.37176, -0.92574, 7.57747]
-        ...             [-0.84200, -0.72521, 7.56150]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Get era data:
-
-        >>> c.gcrs()
-        (array([-4681.50824, -4714.35352]),
-         array([-5030.09119, -4978.74326]),
-         array([0, 454.41493]),
-         array([-0.72067, -0.37378]))
-         array([0.67072, 1.04024]),
-         array([7.57919, 7.56237]))
-        """
-
-        return tuple(self._GCRS if stroke else [i(self._julian) for i in self._GCRS])
-
-    def itrs(self, stroke=False) -> np.ndarray:
-        """Return itrs coordinates.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        tuple
-            Tuple containing x, y, z cartesian positions and velocities in the
-            itrs frame.
-
-            The coordinates are stroke objects if `stroke=True`. Otherwise,
-            1-D arrays are returned.
-
-        See Also
-        --------
-        gcrs : Return gcrs coordinates.
-        geo : Return geographical coordinates.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Get itrs data:
-
-        >>> c.itrs()
-        (array([[6343.70553, 6295.54131]),
-         array([-2641.13728, -2718.36512]),
-         array([-11.44111, 442.89640]),
-         array([-0.37176, -0.84200]),
-         array([-0.92574, -0.72521]),
-         array([7.57747, 7.56150]))
-        """
-
-        return tuple(self._ITRS if stroke else [i(self._julian) for i in self._ITRS])
-
-    def lvlh(self, stroke=False) -> Tuple:
-        """Return the LVLH (Hill frame) coordinates.
-
-        The local-vertical local-horizontal frame (also known as the Hill
-        frame) is a body frame where the z-axis is algigned with the negative
-        of the geocentric position vector, the y-axis is aligned with the
-        negative orbit normal, and the x-axis completes the right handed triad.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        Tuple
-            Tuple containing x, y, z cartesian positions and velocities in the
-            VLVH frame.
-
-            The coordinates are stroke objects if `stroke=True`. Otherwise,
-            1-D arrays are returned.
-
-        Notes
-        -----
-        The LVLH frame definition was taken from NASA's technical memorandum
-        on coordinate frames for the space shuttle program. [NASA1974]_
-
-        References
-        ----------
-        .. [NASA1974] Coordinate Systems for the Space Shuttle Program, Lyndon
-           B. Johnson Space Center, Houston, Texas 77058, Oct 1974, no. NASA
-           TM X-58153.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Get lvlh data:
-
-        >>> c.lvlh()
-        (array([0, -1.13687e-13]),
-         array([-4.54747e-13, 1.06581e-13]),
-         array([-6871.56000, -6871.64511]),
-         array([7.64286, 7.64272]),
-         array([-2.22045e-16, -1.11022e-16]),
-         array([5.55112e-17, -2.83622e-03]))
-        """
-
-        gcrs_position, gcrs_velocity = self._GCRS[:3], self._GCRS[3:]
-
-        Aoi = self._gcrs_to_lvlh_matrix(stroke=True)
-        lvlh_position = np.dot(Aoi, gcrs_position)
-        lvlh_velocity = np.dot(Aoi, gcrs_velocity)
-
-        lvlh = np.concatenate((lvlh_position, lvlh_velocity), axis=0)
-
-        return tuple(lvlh if stroke else [i(self._julian) for i in lvlh.T])
-
-    def _gcrs_to_lvlh_matrix(self, stroke=False) -> np.ndarray:
-        """Return the gcrs to lvlh rotation matrix.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a 2-D array of Stroke objects if true.
-            Otherwise, the output is a 3-D array.
-
-        Returns
-        -------
-        np.ndarray
-            2-D rotation matrix containing Stroke objects if `stroke=True`.
-            Otherwise, 3-D rotation matrix.
-        """
-
-        position, velocity = self._GCRS[:3], self._GCRS[3:]
-
-        norm_position = np.linalg.norm(position)
-        position_cross_velocity = np.cross(position, velocity)
-        n = np.linalg.norm(position_cross_velocity)
-
-        lvlh_z = - np.array([i / norm_position for i in position])
-        lvlh_y = - np.array([i / n for i in position_cross_velocity])
-        lvlh_x = np.cross(lvlh_y, lvlh_z)
-
-        Aoi = np.array([lvlh_x, lvlh_y, lvlh_z])
-
-        if stroke:
-            return Aoi
-        else:
-            a11 = Aoi[0, 0](self._julian)
-            a12 = Aoi[0, 1](self._julian)
-            a13 = Aoi[0, 2](self._julian)
-            a21 = Aoi[1, 0](self._julian)
-            a22 = Aoi[1, 1](self._julian)
-            a23 = Aoi[1, 2](self._julian)
-            a31 = Aoi[2, 0](self._julian)
-            a32 = Aoi[2, 1](self._julian)
-            a33 = Aoi[2, 2](self._julian)
-            return np.array([[a11, a12, a13],
-                             [a21, a22, a23],
-                             [a31, a32, a33]])
-
-    def _get_ang(self, u, v) -> np.ndarray:
-        """Return the degree angle bewteen two vectors.
-
-        Parameters
-        ----------
-        u, v : np.ndarray
-            1-D or 2-D arrays containing row vectors.
-
-            If the dimensions of `u` and `v` do not match, the 1-D array will
-            be broadcast to have the same number of rows as the 2-D array.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing the degree angle between to vector arrays.
-        """
-
-        ua = None if u.ndim == 1 else 1
-        va = None if v.ndim == 1 else 1
-
-        numerator = np.sum(u * v, axis=(ua or va))
-        denominator = np.linalg.norm(u, axis=ua) * np.linalg.norm(v, axis=va)
-        ang = np.degrees(np.arccos(numerator / denominator))
-
-        return ang
-
-    def horizontal(self, location, stroke=False) -> Tuple:
-        """Return horizontal coordinates in decimal degrees.
-
-        The azimuth angle ranges from 0 to 360 degrees, measured clockwise
-        from North. The elevation angle ranges from 0 to 90 degrees measured
-        above the local horizon.
-
-        Parameters
-        ----------
-        location : GroundPosition
-            Eath bound origin of the horizontal system.
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        tuple
-            Tuple containing the elevation and azimuth angles.
-
-            The angles are stroke objects if `stroke=True`. Otherwise, 1-D
-            arrays are returned.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Calculate the horizontal coordinates at (52.16, -106.67):
-
-        >>> location = GroundPosition(latitude=52.16, longitude=-106.67,
-        ...                           height=0.482)
-        >>> c.horizontal(location=location)
-        (array([-40.88076, -39.01006]),
-         array([94.73900, 92.99101]))
-        """
-
-        el_az = [self._elevation(location), self._azimuth(location)]
-
-        return tuple(el_az if stroke else [i(self._julian) for i in el_az])
-
-    def _elevation(self, location) -> Stroke:
-        """Return the elevation angle of a satellite.
-
-        Parameters
-        ----------
-        location : GroundPosition
-            Eath bound origin of the horizontal system.
-
-        Returns
-        -------
-        Stroke
-            Stroke containing the elevation of the satellite.
-        """
-
-        ground_itrs = self._calculate_ground_location_itrs(location)
-        return 90 - self._get_ang(self._ITRS[:3] - ground_itrs, ground_itrs)
-    
-    def _calculate_ground_location_itrs(self, location):
-
-        ground_geographic = np.array([[location.latitude, location.longitude,
-                                       location.height]])
-
-        return self._geo_to_itrs(ground_geographic).reshape((3,))
-
-
-    def _azimuth(self, location) -> Stroke:
-        """Return the azimuth angle of the satellite.
-
-        Parameters
-        ----------
-        location : GroundPosition
-            Eath bound origin of the horizontal system.
-
-        Returns
-        -------
-        Stroke
-            Stroke containing the azimuth of the satellite.
-        """
-
-        latitude, radius = location.latitude, location.radius
-        ground_itrs = self._calculate_ground_location_itrs(location)
-
-        ground_to_satellite = self._ITRS[:3] - ground_itrs
-
-        z_axis_vector = [0, 0, radius / np.sin(np.radians(latitude))]
-        surface_tangent = z_axis_vector - ground_itrs
-
-        sat_to_ground_on_normal = np.sum(ground_to_satellite * ground_itrs)
-        misc_vector = np.array([sat_to_ground_on_normal * ground_itrs[i] for i in range(3)]) / radius ** 2
-        sat_to_ground_on_plane = ground_to_satellite - misc_vector
-
-        direction = np.cross(surface_tangent, ground_itrs)
-        negative_indices = np.sum(sat_to_ground_on_plane * direction) < 0
-        azimuth = self._get_ang(surface_tangent, sat_to_ground_on_plane)
-        azimuth = 360 * negative_indices - 1 * (2 * negative_indices - 1) * azimuth
-
-        return azimuth
-
-    def off_nadir(self, location, stroke=False) -> Union[Stroke, np.ndarray]:
-        """Return off-nadir angle in decimal degrees.
-
-        The off-nadir angle is the angular distance of a ground location from
-        the satellites nadir.
-
-        Parameters
-        ----------
-        location : GroundPosition
-            Ground location of interest.
-        stroke : bool, optional
-            Formats output as a stroke objects if true. Otherwise, the output
-            is a 1-D arrays.
-
-        Returns
-        -------
-        Union[Stroke, np.ndarray]
-            Off-nadir angles in decimal degrees.
-
-            The off-nadir angles are returned as a stroke object if
-            `stroke=True`. Otherwise, a 1-D array is returned.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Calculate the off-nadir data at (52.16, -106.67):
-
-        >>> location = GroundPosition(latitude=52.16, longitude=-106.67,
-        ...                           height=0.482)
-        >>> c.off_nadir(location=location)
-        array([67.10070, 65.29457])
-        """
-
-        ground_itrs = self._calculate_ground_location_itrs(location)
-        ground_to_satellite = (self._ITRS[:3] - ground_itrs).reshape((3,))
-        off_nadir = self._get_ang(ground_to_satellite, self._ITRS[:3])
-
-        return off_nadir if stroke else off_nadir(self._julian)
-
-    def _WGS84_radius(self, latitude) -> np.ndarray:
-        """Return Earth's geocentric radius using WGS84.
-
-        Parameters
-        ----------
-        latitude : np.ndarray
-            1-D array containing latitude in decimal degrees.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing the geocentric radius in kilometres.
-
-        Notes
-        -----
-        The radius is calculated using the WGS84 Earth ellipsoid model as
-        detailed by "Earth Radius by Latitude" by Timur. [Tim18]_
-
-        References
-        ----------
-        .. [Tim18] Timur. Earth Radius by Latitude (WGS 84). 2018. url:
-           https://planetcalc.com/7721/.
-        """
-
-        c, s = np.cos(np.radians(latitude)), np.sin(np.radians(latitude))
-
-        num = (WGS84_MAJOR_AXIS ** 2 * c) ** 2 + (WGS84_MINOR_AXIS ** 2 * s) ** 2
-        denom = (WGS84_MAJOR_AXIS * c) ** 2 + (WGS84_MINOR_AXIS * s) ** 2
-        radius = np.sqrt(num / denom)
-
-        return radius
-
-    def altitude(self, stroke=False) -> Union[Stroke, np.ndarray]:
-        """Return the geodetic altitude in kilometers.
-
-        This method uses the WGS84 reference ellipsoid to calculate the
-        geodetic altitude above the Earth's surface.
-
-        Parameters
-        ----------
-        stroke : bool, optional
-            Formats output as a stroke object if true. Otherwise, the output
-            is a 1-D array.
-
-        Returns
-        -------
-        Union[Stroke, np.ndarray]
-            Geodetic altitude in kilometres.
-
-            The geodetic altitude is returned as a stroke object if
-            `stroke=True`. Otherwise, a 1-D array is returned.
-
-        Notes
-        -----
-        This method uses an ellipsoid based model of the Earth to calculate
-        the ellipsoid height in an iterative manner described in "Coordinate
-        Systems in Geodesy" by E. J. Krakiwsky and D.E. Wells. [KW98b]_
-
-        References
-        ----------
-        .. [KW98b] E. J. Krakiwsky and D. E. Wells. Coordinate Systems in
-           Geodesy. Jan. 1998, pp. 31–33.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Calculate altitude data:
-
-        >>> c.altitude()
-        array([493.42306, 493.59714])
-        """
-
-        a, b = WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS
-
-        tol = 10e-10
-
-        x, y, z, _, _, _ = self._ITRS
-
-        e = np.sqrt(1 - b ** 2 / a ** 2)
-        p = np.sqrt(x ** 2 + y ** 2)
-
-        h = np.sqrt(x ** 2 + y ** 2 + z ** 2) - np.sqrt(a * b)
-        phi = np.arctan((z / p) * (1 - (e ** 2 * a) / (a + h)) ** -1)
-
-        def new_vals(a, b, e, p, N, h, phi, z):
-
-            N = a / np.sqrt(np.cos(phi) ** 2 + b ** 2 / a ** 2 * np.sin(phi) ** 2)
-            h = p / np.cos(phi) - N
-            phi = np.arctan((z / p) * (1 - (e ** 2 * N) / (N + h)) ** -1)
-
-            return N, h, phi
-
-        Np, hp, phip = new_vals(a, b, e, p, a, h, phi, z)
-
-        while (np.mean((hp - h)(self._julian)) > a * tol) or (np.mean((phip - phi)(self._julian)) > tol):
-
-            n, h, phi = Np, hp, phip
-            Np, hp, phip = new_vals(a, b, e, p, n, h, phi, z)
-
-        return h if stroke else h(self._julian)
-
-    def distance(self, location, stroke=False) -> Union[Stroke, np.ndarray]:
-        """Return distance to the ground location.
-
-        Parameters
-        ----------
-        location : GroundPosition
-        stroke : bool, optional
-            Formats output as a stroke object if true. Otherwise, the output
-            is a 1-D array.
-
-        Returns
-        -------
-        Union[Stroke, np.ndarray]
-            Distance in kilometres.
-
-            The distance is returned as a stroke object if `stroke=True`.
-            Otherwise, a 1-D array is returned.
-
-        Examples
-        --------
-        Initialize `Coordinate` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> c = Coordinate(position=position, velocity=velocity, frame="gcrs",
-        ...                julian=julian, offset=2430000)
-
-        Calculate the satellite distances from position (52.16, -106.67):
-
-        >>> location = GroundPosition(latitude=52.16, longitude=-106.67,
-        ...                           height=0.482)
-        >>> c.distance(location=location)
-        array([9070.80824, 8777.02141])
-        """
-
-        ground_itrs = self._calculate_ground_location_itrs(location)
-        distance = np.linalg.norm(self._ITRS[:3] - ground_itrs)
-
-        return distance if stroke else distance(self._julian)
diff --git a/celest/satellite_old/nutation_precession.py b/celest/satellite_old/nutation_precession.py
deleted file mode 100644
index ca59270..0000000
--- a/celest/satellite_old/nutation_precession.py
+++ /dev/null
@@ -1,154 +0,0 @@
-
-
-from celest.satellite._astronomical_quantities import (
-    _calculate_elapsed_JD_century_powers,
-    mean_obliquity,
-    nutation_components,
-    conventional_precession_angles
-)
-import numpy as np
-
-
-def bias_matrix() -> np.ndarray:
-    r"""Generate bias matrix for GCRS and ITRS conversions.
-
-    Returns
-    -------
-    np.ndarray
-        2-D array representing the bias matrix.
-
-    Notes
-    -----
-    The bias matrixs is a 3-2-1 set of Euler angle rotations about
-    :math:`d\\alpha_0=-0.01460`, :math:`\xi_0=-0.0166170`, and
-    :math:`-\eta_0=0.0068192`.[SL13c]_
-
-    References
-    ----------
-    .. [SL13c] M. Soffel and R. Langhans. Space-Time Reference Systems.
-       Astronomy and Astrophysics Library. Springer-Verlag, 2013, pp. 197-233.
-    """
-
-    xi, eta, dalpha = -0.0166170, -0.0068192, -0.01460
-
-    a1 = np.radians(dalpha / 3600)
-    a2 = np.radians(xi / 3600)
-    a3 = np.radians(-eta / 3600)
-
-    s1, c1 = np.sin(a1), np.cos(a1)
-    s2, c2 = np.sin(a2), np.cos(a2)
-    s3, c3 = np.sin(a3), np.cos(a3)
-
-    matrix = np.zeros((3, 3))
-    matrix[0, 0] = c1 * c2
-    matrix[0, 1] = s1 * c2
-    matrix[0, 2] = - s2
-    matrix[1, 0] = - s1 * c3 + c1 * s2 * s3
-    matrix[1, 1] = c1 * c3 + s1 * s2 * s3
-    matrix[1, 2] = c2 * s3
-    matrix[2, 0] = s1 * s3 + c1 * s2 * c3
-    matrix[2, 1] = - c1 * s3 + s1 * s2 * c3
-    matrix[2, 2] = c2 * c3
-
-    return matrix
-
-
-def precession_matrix(julian: np.ndarray) -> np.ndarray:
-    """Return precession tensor.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        3-D array of shape (n,3,3) containing 3x3 precession matrices
-        associated with each time along the 0 axis.
-
-    Notes
-    -----
-    The precession matrix is implemented as defined in "Space-Time Reference
-    Systems". [SL13c]_
-
-    References
-    ----------
-    .. [SL13c] M. Soffel and R. Langhans. Space-Time Reference Systems.
-       Astronomy and Astrophysics Library. Springer-Verlag, 2013, pp. 197-233.
-    """
-
-    zeta, theta, z = conventional_precession_angles(julian=julian)
-
-    a1 = np.radians(zeta / 3600)
-    a2 = np.radians(theta / 3600)
-    a3 = - np.radians(z / 3600)
-
-    s1, c1 = np.sin(a1), np.cos(a1)
-    s2, c2 = np.sin(a2), np.cos(a2)
-    s3, c3 = np.sin(a3), np.cos(a3)
-
-    matrix = np.zeros((julian.size, 3, 3))
-    matrix[:, 0, 0] = - s1 * s3 + c1 * c2 * c3
-    matrix[:, 0, 1] = c1 * s3 + s1 * c2 * c3
-    matrix[:, 0, 2] = - s2 * c3
-    matrix[:, 1, 0] = - s1 * c3 - c1 * c2 * s3
-    matrix[:, 1, 1] = c1 * c3 - s1 * c2 * s3
-    matrix[:, 1, 2] = s2 * s3
-    matrix[:, 2, 0] = c1 * s2
-    matrix[:, 2, 1] = s1 * s2
-    matrix[:, 2, 2] = c2
-
-    return matrix
-
-
-def nutation_matrix(julian: np.ndarray) -> np.ndarray:
-    """Return nutation tensor.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        1-D array containing Julian times in the J2000 epoch.
-
-    Returns
-    -------
-    np.ndarray
-        3-D array of shape (n,3,3) containing 3x3 nutation matrices associated
-        with each time along the 0 axis.
-
-    Notes
-    -----
-    The nutation matrix is implemented as defined in "Space-Time Reference
-    Systems". [SL13c]_
-
-    References
-    ----------
-    .. [SL13c] M. Soffel and R. Langhans. Space-Time Reference Systems.
-       Astronomy and Astrophysics Library. Springer-Verlag, 2013, pp. 197-233.
-    """
-
-    T1, T2, T3 = _calculate_elapsed_JD_century_powers(julian, 3)
-
-    nutation_in_longitude, nutation_in_obliquity = nutation_components(julian)
-    eps_A = 3600 * mean_obliquity(julian) - 46.84024 * T1 - 0.00059 * T2 + 0.001813 * T3
-
-    a1 = np.radians(eps_A / 3600)
-    a2 = - np.radians(nutation_in_longitude / 3600)
-    a3 = - np.radians(eps_A / 3600 + nutation_in_obliquity / 3600)
-
-    s1, c1 = np.sin(a1), np.cos(a1)
-    s2, c2 = np.sin(a2), np.cos(a2)
-    s3, c3 = np.sin(a3), np.cos(a3)
-
-    matrix = np.zeros((julian.size, 3, 3))
-    matrix[:, 0, 0] = c2
-    matrix[:, 0, 1] = s1 * s2
-    matrix[:, 0, 2] = - c1 * s2
-    matrix[:, 1, 0] = s2 * s3
-    matrix[:, 1, 1] = c1 * c3 - s1 * c2 * s3
-    matrix[:, 1, 2] = s1 * c3 + c1 * c2 * s3
-    matrix[:, 2, 0] = s2 * c3
-    matrix[:, 2, 1] = - c1 * s3 - s1 * c2 * c3
-    matrix[:, 2, 2] = - s1 * s3 + c1 * c2 * c3
-
-    return matrix
diff --git a/celest/satellite_old/satellite.py b/celest/satellite_old/satellite.py
deleted file mode 100644
index ec9963e..0000000
--- a/celest/satellite_old/satellite.py
+++ /dev/null
@@ -1,214 +0,0 @@
-
-
-from celest.satellite.coordinate import Coordinate
-from typing import Tuple
-import numpy as np
-import pandas as pd
-
-
-class Satellite(Coordinate):
-    """Satellite(position, velocity, frame, julian, offset=0)
-
-    Satellite abstraction for satellite-ground encounters.
-
-    `julian + offset` is the Julian time in the J2000 epoch associated with
-    input positions.
-
-    Parameters
-    ----------
-    position : array_like
-        2-D array containing position coordinates.
-
-        Supported coordinate frames include the Geocentric Celestial Reference
-        System (gcrs) and International Terrestrial Reference System (itrs).
-    velocity : array_like
-        2-D array containing velocity coordinates.
-
-        Refer to the `position` parameter for supported coordinate frames.
-    frame : {"gcrs", "geo", "itrs"}
-        Frame specifier for the input position.
-    julian : array_like
-        1-D array containing time data in Julian days.
-
-        If times are not in the J2000 epoch, a non-zero offset must be passed
-        in to add to the julian times.
-    offset : float, optional
-        Offset to convert input time data to the J2000 epoch, default is zero.
-
-    Examples
-    --------
-    Initialize `Satellite` using gcrs data:
-
-    >>> julian = [30462.50000, 30462.50069]
-    >>> position = [[-4681.50824, -5030.09119, 000.00000]
-    ...             [-4714.35352, -4978.74326, 454.41493]]
-    >>> velocity = [[-0.72067, 0.67072, 7.57919]
-    ...             [-0.37378, 1.04024, 7.56237]]
-    >>> s = Satellite(position=position, velocity=velocity, frame="gcrs",
-    ...               julian=julian, offset=2430000)
-    """
-
-    def __init__(self, position, velocity, frame, julian, offset=0) -> None:
-
-        super().__init__(position, velocity, frame, julian, offset)
-
-    def __len__(self) -> int:
-
-        return self._length
-
-    def attitude(self, location, stroke=False) -> Tuple:
-        """Return satellite roll, pitch, and yaw angles in decimal degrees.
-
-        This method returns the attitude angles required to align the
-        satellite's nadir with a ground location for ground imaging.
-
-        Parameters
-        ----------
-        location : GroundPosition
-            Location for satellite attitude pointing.
-        stroke : bool, optional
-            Formats output as a tuple of stroke objects if true. Otherwise,
-            the output is a tuple of 1-D arrays.
-
-        Returns
-        -------
-        Tuple
-            Tuple containing roll, pitch, and yaw angles in decimal degrees.
-
-        Notes
-        -----
-        The methods used were taken from [adcs1]_.
-
-        References
-        ----------
-        .. [adcs1] G. H. J. van Vuuren, “The design and simulation analysis of
-           an attitude determination and control system for a small earth
-           observation satellite,” Master of Engineering, Stellenbosch
-           University, Stellenbosch, South Africa, Mar 2015.
-
-        Examples
-        --------
-        Initialize `Satellite` using gcrs data:
-
-        >>> julian = [30462.50000, 30462.50069]
-        >>> position = [[-4681.50824, -5030.09119, 000.00000]
-        ...             [-4714.35352, -4978.74326, 454.41493]]
-        >>> velocity = [[-0.72067, 0.67072, 7.57919]
-        ...             [-0.37378, 1.04024, 7.56237]]
-        >>> s = Satellite(position=position, velocity=velocity, frame="gcrs",
-        ...               julian=julian, offset=2430000)
-
-        Get attitude data:
-
-        >>> location = GroundPosition(latitude=43.65, longitude=-79.38,
-        ...                           height=0.076)
-        >>> roll, pitch, yaw = s.attitude(location)
-        """
-
-        x, y, z, _, _, _ = [i.reshape((-1, 1)) for i in self.lvlh()]
-        satellite_position = np.concatenate((x, y, z), axis=1)
-
-        ground_itrs = self._calculate_ground_location_itrs(location).reshape((1, -1))
-        repeated_ground_itrs = np.repeat(ground_itrs, self._length, axis=0)
-        ground_gcrs = self._itrs_to_gcrs(repeated_ground_itrs)
-
-        Aoi = self._gcrs_to_lvlh_matrix()
-        ground_lvlh = np.einsum('jki, ik -> ij', Aoi, ground_gcrs)
-
-        satellite_position_norm = np.linalg.norm(satellite_position, axis=1)
-        s = - satellite_position / satellite_position_norm[:, None]
-
-        ground_norm = np.linalg.norm(ground_lvlh - satellite_position, axis=1)
-        g = (ground_lvlh - satellite_position) / ground_norm[:, None]
-
-        v, c = np.cross(s, g, axis=1), np.sum(s * g, axis=1)
-
-        a32 = v[:, 0] + v[:, 1] * v[:, 2] / (1 + c)
-        a31 = - v[:, 1] + v[:, 0] * v[:, 2] / (1 + c)
-        a33 = 1 - (v[:, 0] ** 2 + v[:, 1] ** 2) / (1 + c)
-        a12 = - v[:, 2] + v[:, 0] * v[:, 1] / (1 + c)
-        a22 = 1 - (v[:, 0] ** 2 + v[:, 2] ** 2) / (1 + c)
-
-        roll = - np.degrees(np.arcsin(a32))
-        pitch = np.degrees(np.arctan2(a31, a33))
-        yaw = np.degrees(np.arctan2(a12, a22))
-
-        if stroke:
-            return tuple([self._stroke_init(i) for i in (roll, pitch, yaw)])
-        else:
-            return roll, pitch, yaw
-
-    def save(self, times=("julian",), positions=("gcrs",), path=None,
-             sep=",", float_format=None) -> None:
-        """Export satellite data.
-
-        Parameters
-        ----------
-        times : Tuple, optional
-            Tuple containing time representations to save.
-
-            Available times include "julian", "ut1", "gmst", and "gast".
-        positions : Tuple, optional
-            Tuple containing position representations to save.
-
-            Available positions include "gcrs", "itrs", and "geo".
-        path : str, optional
-            String or path object to save data to. If not provided, data is
-            returned as a string.
-        sep : str, optional
-            String of length 1 representing the field delimiter for the
-            output file.
-        float_format : str, optional
-            Format string for floating point numbers.
-
-        Examples
-        --------
-        >>> positions = ["gcrs", "itrs"]
-        >>> finch.save(fname="data.csv", positions=positions, delimiter=",")
-        """
-
-        key_mapping = {
-            "julian": (self.julian, "Julian.J2000"),
-            "ut1": (self.ut1, "UT1.hr"),
-            "gmst": (self.gmst, "GMST.hr"),
-            "gast": (self.gast, "GAST.hr"),
-            "itrs": (self.itrs, "ITRS.x.km", "ITRS.y.km", "ITRS.z.km",
-                     "ITRS.vx.m/s", "ITRS.vy.m/s", "ITRS.vz.m/s"),
-            "gcrs": (self.gcrs, "GCRS.x.km", "GCRS.y.km", "GCRS.z.km",
-                     "GCRS.vx.m/s", "GCRS.vy.m/s", "GCRS.vz.m/s"),
-            "geo": (self.geo, "GEO.lat.deg", "GEO.lon.deg", "GEO.height.km")
-        }
-
-        for time in times:
-            if time not in key_mapping:
-                raise ValueError(f"{time} is not a valid time type.")
-
-        for position in positions:
-            if position not in key_mapping:
-                raise ValueError(f"{position} is not a valid position type.")
-
-        data = {}
-
-        if isinstance(times, str):
-            times = (times,)
-
-        for time in times:
-            time_information = key_mapping[time]
-            time_data = time_information[0]()
-            column = time_information[1]
-
-            data[column] = time_data
-
-        if isinstance(positions, str):
-            positions = (positions,)
-
-        for position in positions:
-            position_information = key_mapping[position]
-            position_data = position_information[0]()
-            columns = position_information[1:]
-
-            for i, column in enumerate(columns):
-                data[column] = position_data[i]
-
-        df = pd.DataFrame(data)
-        df.to_csv(path_or_buf=path, sep=sep, float_format=float_format)
diff --git a/celest/satellite_old/time.py b/celest/satellite_old/time.py
deleted file mode 100644
index 73d14bd..0000000
--- a/celest/satellite_old/time.py
+++ /dev/null
@@ -1,499 +0,0 @@
-
-
-from celest.satellite._astronomical_quantities import (
-    equation_of_time,
-    from_julian,
-    _calculate_elapsed_JD_century_powers,
-    sun_right_ascension,
-    equation_of_equinoxes
-)
-from datetime import datetime, timedelta
-import numpy as np
-
-
-class Time:
-    """Time(julian, offset=0)
-
-    Time transformations.
-
-    `julian + offset` is the Julian time in J2000 epoch which can be converted
-    into different time representations useful for astronomical calculations.
-
-    Parameters
-    ----------
-    julian : array_like
-        1-D array containing time data in Julian days.
-
-        If time data is not in the J2000 epoch, an offset must be pased in to
-        be added to the julian times.
-    offset : float, optional
-        Offset to convert input time data to the J2000 epoch, default to zero.
-
-    Examples
-    --------
-    Construct a `Time` object using J2000 input times:
-
-    >>> time = [2460462.50, 2460462.91, 2460463.02]
-    >>> t = Time(julian=time, offset=0)
-    >>> gmst = t.gmst()
-
-    For modified Julian input, we add an offset `2430000`:
-
-    >>> time = [30462.50, 30462.91, 30463.02]
-    >>> t = Time(julian=time, offset=2430000)
-
-    We can then convert to various time representations:
-
-    >>> gmst = t.gmst()
-    >>> last = t.last(longitude=145)
-    """
-
-    def __init__(self, julian, offset=0) -> None:
-
-        time = np.array(julian)
-
-        if time.ndim == 0:
-            time = time.reshape((1))
-        elif time.ndim != 1:
-            raise ValueError("Time array must be 1-D.")
-
-        self._julian = time + offset
-        self._length = self._julian.size
-
-    def __len__(self):
-
-        return self._length
-
-    def true_solar_time(self, longitude) -> np.ndarray:
-        """Return true solar time in decimal hours.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing true solar time in decimal hours.
-
-        See Also
-        --------
-        mean_solar_time :
-            Return mean solar time in decimal hours.
-
-        Notes
-        -----
-        The true solar time, :math:`TT_s` can be calculated by converting
-        Julian times, :math:`JD`, into UTC as follows:
-
-        .. math:: UTC = 24\left(JD\%1\\right)^h + \\alpha^h
-
-        The true solar time can then be calculated using the Equation of Time,
-        :math:`EoT`, definition and the mean solar time, :math:`MT_s`. [SL13a]_
-
-        .. math:: TT_s = MT_s + EoT
-
-        .. math:: TT_s = (UTC^h + (EoT^\circ + \phi^\circ)/15)\%24
-
-        References
-        ----------
-        .. [SL13a] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 203.
-
-        Examples
-        --------
-        Calculate the true solar time for a constant longitude:
-
-        >>> Time(julian=[2456293.54167]).true_solar_time(longitude=147.46)
-        np.array([10.773109])
-        """
-
-        eqn_of_time = equation_of_time(julian=self._julian)
-        utc_time = 24 * (self._julian % 1) - 12
-        true_solar_time = (utc_time + (np.array(longitude) + eqn_of_time) / 15) % 24
-
-        return true_solar_time
-
-    def mean_solar_time(self, longitude) -> np.ndarray:
-        """Return mean solar time in decimal hours.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing mean solar time in decimal hours.
-
-        See Also
-        --------
-        true_solar_time :
-            Return true solar time in decimal hours.
-
-        Notes
-        -----
-        The mean solar time, :math:`MT_s` can be calculated by converting
-        Julian times, :math:`JD`, into UTC as follows:
-
-        .. math:: UTC = 24\left(JD\%1\\right)^h + \\alpha^h
-
-        The mean solar time can then be calculated using the Equation of Time,
-        :math:`EoT`, and the local meridians longitude, :math:`\phi` as
-        follows:
-
-        .. math:: MT_s = (UTC^h + \phi^\circ/15)\%24
-
-        Examples
-        --------
-        Calculate the mean solar time for a constant longitude:
-
-        >>> Time(julian=2456293.54167).mean_solar_time(longitude=147.46)
-        np.array([10.830667])
-        """
-
-        utc_time = 24 * (self._julian % 1) - 12
-        mean_solar_time = (utc_time + np.array(longitude) / 15) % 24
-
-        return mean_solar_time
-
-    def true_hour_angle(self, longitude) -> np.ndarray:
-        """Return true hour angle in decimal hours.
-
-        The true hour angle is the angle between the Sun's apparent position at
-        a given time and its position at local solar noon. The value falls
-        between `0` and `360`, measured in increasing degrees past local solar
-        noon.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing true solar hour angles in decimal hours.
-
-        See Also
-        --------
-        mean_hour_angle :
-            Return mean hour angle in decimal hours.
-
-        Notes
-        -----
-        The true hour angle is zero at local solar noon and increases with the
-        true solar time by a rate of :math:`15^\circ` per hour. Thus, the true
-        hour angle can be calculated from the following:
-
-        .. math:: h_{Sun}^\circ = TT_s^\circ - 180^\circ
-
-        where :math:`TT_s^\circ` is the true solar time in degrees. [SL13b]_
-
-        References
-        ----------
-        .. [SL13b] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 203.
-
-        Examples
-        --------
-        Calculate the true hour angle at different longitudes:
-
-        >>> julian = [2455368.75, 2459450.85]
-        >>> lon = [-105, -118.24]
-        >>> Time(julian=julian).true_hour_angle(longitude=lon)
-        np.array([10.97157662,
-                  12.47807655])
-        """
-
-        return (self.true_solar_time(longitude) - 12) % 24
-
-    def mean_hour_angle(self, longitude) -> np.ndarray:
-        """Return mean hour angle in decimal hours.
-
-        The mean hour angle is the angle between the Sun's mean position at
-        a given time and its position at local solar noon. The value falls
-        between `0` and `360`, measured in increasing degrees past local solar
-        noon.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D containing mean solar hour angles in decimal hours.
-
-        See Also
-        --------
-        true_hour_angle :
-            Return true hour angle in decimal hours.
-
-        Notes
-        -----
-        The mean hour angle is zero at local solar noon and increases with the
-        mean solar time by a rate of :math:`15^\circ` per hour. Thus, the mean
-        hour angle can be calculated from the following:
-
-        .. math:: h_{hSun}^\circ = MT_s^\circ - 180^\circ
-
-        where :math:`MT_s^\circ` is the mean solar time in degrees. [SL13c]_
-
-        References
-        ----------
-        .. [SL13c] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 203.
-
-        Examples
-        --------
-        Calculate the meann hour angle for a single longitude:
-
-        >>> Time(julian=2456293.5416666665).mean_hour_angle(longitude=147.46)
-        np.array([22.83066666])
-        """
-
-        return (self.mean_solar_time(longitude) - 12) % 24
-
-    def ut1(self) -> np.ndarray:
-        """Return the universal time (same as GMT) in decimal hours.
-
-        Due to approximations in the mean solar time calculations, the DUT1
-        time correction is not accounted for which will introduce an error of
-        at most 1.8 seconds.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing universal time in decimal hours.
-
-        Notes
-        -----
-        The universal time is equal to the mean solar time at the Greenwich
-        meridian. [SL13d]_ It can be calculated using
-        `Time.mean_solar_time(longitude=0)`.
-
-        References
-        ----------
-        .. [SL13d] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 203.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).ut1()
-        np.array([6.00000,
-                  8.40000,
-                  1.00000])
-        """
-
-        return self.mean_solar_time(0)
-
-    def julian(self) -> np.ndarray:
-        """Convenience method to return Julian times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing the Julian times in decimal days.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).julian()
-        np.array([2455368.75,
-                  2459450.85,
-                  2456293.5416666665])
-        """
-
-        return self._julian
-
-    def datetime(self) -> np.ndarray:
-        """Return `datetime.datetime` object array.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing `datetime.datetime` objects.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).datetime()
-        np.array([datetime.datetime(2010, 6, 21, 6, 0)
-                  datetime.datetime(2021, 8, 24, 8, 24, 0, 8)
-                  datetime.datetime(2013, 1, 1, 0, 59, 59, 999987)])
-        """
-
-        year, month, day = from_julian(self._julian)
-        full_days, fractional_day = day.astype(int), day % 1
-
-        datetime_list = []
-        for y, m, d, f in zip(year, month, full_days, fractional_day):
-            integer_date = datetime(year=y, month=m, day=d)
-            extra_date = timedelta(days=f)
-            datetime_list.append(integer_date + extra_date)
-
-        return np.array(datetime_list)
-
-    def gmst(self) -> np.ndarray:
-        """Return Greenwich Mean Sidereal Time in decimal hours.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing Greenwich Mean Sidereal Time in hours and
-            decimals.
-
-        Notes
-        -----
-        The Greenwich Mean Sidereal Time is calculated using the methods
-        described in Astronomical Algorithms by Jean Meeus. [Mee98a]_
-
-        References
-        ----------
-        .. [Mee98a] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 87 - 88. isbn: 9780943396613.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).gmst()
-        np.array([23.955316,
-                  6.589391,
-                  7.723214])
-        """
-
-        T1, T2, T3 = _calculate_elapsed_JD_century_powers(self._julian, 3)
-
-        gmst = 280.46061837 + 360.98564736629 * T1 * 36525 + \
-            0.000387933 * T2 - T3 / 38710000
-
-        return gmst % 360 / 15
-
-    def lmst(self, longitude) -> np.ndarray:
-        """Return Local Mean Sidereal Time in decimal hours.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing Local Mean Sidereal Time in decimal hours.
-
-        Notes
-        -----
-        The Local Mean Sidereal Time can be calculated using the following
-        formulation:
-
-        .. math:: LMST^h = h^h_{mSun} + \\alpha^h_{mSun}
-
-        where :math:`h^h_{mSun}` is the mean solar hour angle at the observer's
-        longitude and :math:`\\alpha^h_{mSun}` is the right ascension of the
-        mean Sun position. [SL13e]_
-
-        References
-        ----------
-        .. [SL13e] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 205.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).lmst(longitude=150)
-        np.array([9.98419514,
-                  16.62892539,
-                  17.78123885])
-        """
-
-        hour_angle = self.mean_hour_angle(longitude)
-        right_ascension_of_sun = sun_right_ascension(self._julian) / 15
-
-        return (hour_angle + right_ascension_of_sun) % 24
-
-    def gast(self) -> np.ndarray:
-        """Return Greenwich Apparent Sidereal Time in decimal hours.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing Greenwich Apparent Sidereal Time in hours and
-            degrees.
-
-        Notes
-        -----
-        The difference between apparent and mean solar position is the equation
-        of time. Thus, the Greenwich Apparent Sidereal Time can be calculated
-        as follows:
-
-        .. math:: GAST^h = GMST^h + EoT^h
-
-        where :math:`EoT` is the equation of time in decimal hours.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).gast()
-        np.array([23.955596,
-                  6.589146,
-                  7.723465])
-        """
-
-        gmst = self.gmst()
-        eqn_of_equinoxes = equation_of_equinoxes(self._julian) / 3600
-
-        return gmst + eqn_of_equinoxes
-
-    def last(self, longitude) -> np.ndarray:
-        """Return Local Apparent Sidereal Time in hours and degrees.
-
-        Parameters
-        ----------
-        longitude : array_like
-            1-D array containing longitude in decimal degrees. If a scalar
-            value is given, it will be applied for all times.
-
-        Returns
-        -------
-        np.ndarray
-            1-D array containing Local Apparent Sidereal Time in hours and
-            degrees.
-
-        Notes
-        -----
-        The Local Apparent Sidereal Time is calculated using the methods
-        described in Space-Time Reference Systems by M. Soffel and R.
-        Langhans. [SL13f]_
-
-        References
-        ----------
-        .. [SL13f] M. Soffel and R. Langhans. Space-Time Reference Systems.
-           Astronomy and Astrophysics Library. Springer-Verlag, 2013, p. 205.
-
-        Examples
-        --------
-        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        >>> Time(julian=julian).last(longitude=150)
-        np.array([9.98447567,
-                  16.6286799,
-                  17.78148932])
-        """
-
-        ut1 = self.ut1()
-        right_ascendion_of_sun = sun_right_ascension(self._julian) / 15
-        eqn_of_equinoxes = equation_of_equinoxes(self._julian) / 3600
-
-        last = ut1 - 12 + right_ascendion_of_sun + eqn_of_equinoxes + np.array(longitude) / 15
-
-        return last % 24
diff --git a/celest/schedule/alns.py b/celest/schedule/alns.py
index 2429732..0a4ab71 100644
--- a/celest/schedule/alns.py
+++ b/celest/schedule/alns.py
@@ -214,7 +214,7 @@ def _get_repair_function(self, repair_function_index) -> Callable:
 
     def _accept(self, temporary_solution: Any, current_solution: Any,
                 T: float) -> bool:
-        """Accpetance criterion.
+        """Acceptance criterion.
 
         This method implements a simulated annealing acceptance criterion.
 
diff --git a/celest/schedule/insertion_operators.py b/celest/schedule/insertion_operators.py
index 9fcbc95..61457cb 100644
--- a/celest/schedule/insertion_operators.py
+++ b/celest/schedule/insertion_operators.py
@@ -1,16 +1,17 @@
 
 
+from celest.units.quantity import Quantity
+from celest import units as u
 import math
+import numpy as np
 
 
 def greedy_insertion(request_handler, number_to_insert):
-
     request_handler.sort_by_decreasing_priority()
     insert_first_n_requests(request_handler, number_to_insert)
 
 
 def insert_first_n_requests(request_handler, number_to_insert):
-
     for request_index in range(request_handler.number_of_requests):
 
         if not number_to_insert:
@@ -25,12 +26,18 @@ def insert_first_n_requests(request_handler, number_to_insert):
         look_angle = request_handler.look_angle(request_index)
         request_vtw_list = request_handler.vtw_list(request_index)
 
-        duration_in_days = duration / 86400
+        duration_in_days = Quantity(duration.to(u.s).data / 86400, u.jd2000)
 
         for vtw_index, vtw in enumerate(request_vtw_list):
 
             if look_angle is not None:
-                start = _look_angle_time(vtw.pitch, look_angle, vtw.rise_time, vtw.set_time)
+                start = _look_angle_time(
+                    vtw.attitude.pitch,
+                    look_angle,
+                    vtw.attitude.time,
+                    vtw.rise_time,
+                    vtw.set_time
+                )
                 if start is None:
                     continue
             else:
@@ -42,7 +49,7 @@ def insert_first_n_requests(request_handler, number_to_insert):
                 continue
             if not image_quality_is_met(vtw, start, minimum_image_quality):
                 continue
-            if does_conflict_exists(request_handler, start, duration):
+            if does_conflict_exists(request_handler, start, duration_in_days):
                 continue
 
             request_handler.schedule_request(request_index, vtw_index, start, duration)
@@ -50,69 +57,66 @@ def insert_first_n_requests(request_handler, number_to_insert):
             break
 
 
-def _look_angle_time(look_angle, desired_look_angle, start, end, tol=1e-6):
+def _look_angle_time(look_angle, desired_look_angle, julian, start, end):
+    julian = julian.to(u.jd2000).data
+    look_angle = look_angle.to(u.deg).data - desired_look_angle.to(u.deg).data
 
-    look_angle = look_angle - desired_look_angle
+    left_value_index = np.where(julian == start.to(u.jd2000).data)[0][0]
+    right_value_index = np.where(julian == end.to(u.jd2000).data)[0][0]
 
-    left_value = start
-    right_value = end
-
-    left_look_angle = look_angle(left_value)
-    right_look_angle = look_angle(right_value)
+    left_look_angle = look_angle[left_value_index]
+    right_look_angle = look_angle[right_value_index]
 
     if left_look_angle * right_look_angle > 0:
         return None
 
-    while (right_value - left_value > tol):
-        mid_value = (left_value + right_value) / 2
-        mid_look_angle = look_angle(mid_value)
+    while right_value_index - left_value_index > 1:
+        mid_value_index = (left_value_index + right_value_index) // 2
+        mid_look_angle = look_angle[mid_value_index]
         if left_look_angle * mid_look_angle > 0:
-            left_value = mid_value
+            left_value_index = mid_value_index
             left_look_angle = mid_look_angle
         else:
-            right_value = mid_value
-            right_look_angle = mid_look_angle
+            right_value_index = mid_value_index
 
-    return (left_value + right_value) / 2
+    return Quantity(julian[left_value_index], u.jd2000)
 
 
 def image_quality_is_met(vtw, start, minimum_image_quality):
-
     return image_quality(vtw, start) >= minimum_image_quality
 
 
 def image_quality(vtw, start):
-
-    nadir_time = (vtw.rise_time + vtw.set_time) / 2
-    return math.floor(10 - 9 * abs(start - nadir_time) / (nadir_time - vtw.rise_time))
+    nadir_time = ((vtw.rise_time + vtw.set_time) / 2).to(u.jd2000).data
+    start_time = start.to(u.jd2000).data
+    ride_time = vtw.rise_time.to(u.jd2000).data
+    return math.floor(10 - 9 * abs(start_time - nadir_time) / (nadir_time - ride_time))
 
 
 def does_conflict_exists(request_handler, start, duration_in_days):
-
     for request in request_handler:
         if not request[0]:
             continue
         scheduled_start = request[2]
-        scheduled_end = scheduled_start + request[3] / 86400
+        scheduled_duration = Quantity(request[3].to(u.s).data / 86400, u.jd2000)
+        scheduled_end = scheduled_start + scheduled_duration
 
-        if start < scheduled_start and scheduled_start < start + duration_in_days:
+        if start <= scheduled_start < start + duration_in_days:
             return True
-        if start < scheduled_end and scheduled_end < start + duration_in_days:
+        if start < scheduled_end <= start + duration_in_days:
             return True
-        if scheduled_start < start and start + duration_in_days < scheduled_end:
+        if scheduled_start <= start and start + duration_in_days <= scheduled_end:
             return True
 
     return False
 
 
 def minimum_opportunity_insertion(request_handler, number_to_insert):
-
     request_handler.sort_by_decreasing_opportunity()
     insert_first_n_requests(request_handler, number_to_insert)
 
 
 def minimum_conflict_insertion(request_handler, number_to_insert):
-
     request_handler.sort_vtws_by_increasing_conflict_degree()
     insert_first_n_requests(request_handler, number_to_insert)
 
diff --git a/celest/schedule/request_handler.py b/celest/schedule/request_handler.py
index 86036c4..04244a3 100644
--- a/celest/schedule/request_handler.py
+++ b/celest/schedule/request_handler.py
@@ -1,3 +1,5 @@
+
+
 from enum import IntEnum
 
 
@@ -19,22 +21,18 @@ class RequestIndices(IntEnum):
 class RequestHandler:
 
     def __init__(self):
-
         self.requests = []
         self.number_of_requests = 0
         self.number_of_scheduled_requests = 0
 
     def __getitem__(self, index):
-
         return self.requests[index]
 
     def __iter__(self):
-
         self.iteration_index = 0
         return self
 
     def __next__(self):
-
         current_index = self.iteration_index
         self.iteration_index += 1
         if self.iteration_index < self.number_of_requests:
@@ -42,13 +40,11 @@ def __next__(self):
         raise StopIteration
 
     def add_request(self, location, deadline, duration, priority, quality, look_ang, vtws):
-
         self.number_of_requests += 1
         self.requests.append([False, None, None, None, location, deadline,
                               duration, priority, quality, look_ang, vtws])
 
     def schedule_request(self, request_index, vtw_index, start, duration):
-
         self.requests[request_index][RequestIndices.is_scheduled] = True
         self.requests[request_index][RequestIndices.vtw_index] = vtw_index
         self.requests[request_index][RequestIndices.scheduled_start_time] = start
@@ -56,14 +52,12 @@ def schedule_request(self, request_index, vtw_index, start, duration):
         self.number_of_scheduled_requests += 1
 
     def unschedule_all_requests(self):
-
         for index in range(self.number_of_requests):
             self.unschedule_request(index)
 
         self.number_of_scheduled_requests = 0
 
     def unschedule_request(self, request_index):
-
         self.requests[request_index][RequestIndices.is_scheduled] = False
         self.requests[request_index][RequestIndices.vtw_index] = None
         self.requests[request_index][RequestIndices.scheduled_start_time] = None
@@ -72,57 +66,44 @@ def unschedule_request(self, request_index):
         self.number_of_scheduled_requests -= 1
 
     def is_request_scheduled(self, request_index):
-
         return self.requests[request_index][RequestIndices.is_scheduled]
 
     def vtw_index(self, request_index):
-
         return self.requests[request_index][RequestIndices.vtw_index]
 
     def start_time(self, request_index):
-
         return self.requests[request_index][RequestIndices.scheduled_start_time]
 
     def scheduled_duration(self, request_index):
-
         return self.requests[request_index][RequestIndices.scheduled_duration]
 
     def deadline(self, request_index):
-
         return self.requests[request_index][RequestIndices.deadline]
 
     def duration(self, request_index):
-
         return self.requests[request_index][RequestIndices.duration]
 
     def priority(self, request_index):
-
         return self.requests[request_index][RequestIndices.priority]
 
     def image_quality(self, request_index):
-
         return self.requests[request_index][RequestIndices.quality]
 
     def look_angle(self, request_index):
-
         return self.requests[request_index][RequestIndices.look_angle]
 
     def vtw_list(self, request_index):
-
         return self.requests[request_index][RequestIndices.vtw_list]
 
     def sort_by_decreasing_priority(self):
-
         key = lambda x: x[RequestIndices.priority]
         self.requests = sorted(self.requests, key=key, reverse=True)
 
     def sort_by_decreasing_opportunity(self):
-
         key = lambda x: len(x[RequestIndices.vtw_list])
         self.requests = sorted(self.requests, key=key, reverse=True)
 
     def sort_by_decreasing_conflict_degree(self):
-
         conflict_degree = []
         conflict_degree_index = []
         for request_index, request in enumerate(self.requests):
@@ -135,23 +116,23 @@ def sort_by_decreasing_conflict_degree(self):
 
             conflict_degree_index.append(request_index)
 
-        indices = [index for _, index in sorted(zip(conflict_degree, conflict_degree_index))]
+        indices = [index for _, index in
+                   sorted(zip(conflict_degree, conflict_degree_index))]
         indices.reverse()
 
         self.requests = [self.requests[index] for index in indices]
 
     def _conflict_degree(self, vtw) -> float:
-
         over_lapping_vtws = self._set_of_overlapping_vtws(vtw)
         total_overlapping_time = 0
 
         for other_vtw in over_lapping_vtws:
-            total_overlapping_time += self._julian_day_overlap_between_vtws(vtw, other_vtw)
+            total_overlapping_time +=\
+                self._julian_day_overlap_between_vtws(vtw, other_vtw)
 
         return total_overlapping_time / len(over_lapping_vtws)
 
     def _set_of_overlapping_vtws(self, current_vtw) -> list:
-
         overlap_vtws = []
 
         for request in self.requests:
@@ -163,14 +144,9 @@ def _set_of_overlapping_vtws(self, current_vtw) -> list:
         return overlap_vtws
 
     def _are_vtws_overlaping(self, vtw1, vtw2) -> bool:
-
-        if (vtw1.set_time < vtw2.rise_time) | (vtw2.set_time < vtw1.rise_time):
-            return False
-        else:
-            return True
+        return not (vtw1.set_time < vtw2.rise_time) | (vtw2.set_time < vtw1.rise_time)
 
     def _julian_day_overlap_between_vtws(self, vtw1, vtw2) -> float:
-
         s1, s2 = vtw1.rise_time, vtw2.rise_time
         e1, e2 = vtw1.set_time, vtw2.set_time
 
@@ -186,17 +162,17 @@ def _julian_day_overlap_between_vtws(self, vtw1, vtw2) -> float:
             return e2 - s1
 
     def sort_vtws_by_increasing_conflict_degree(self):
-
         for request_index in range(self.number_of_requests):
             conflict_degree = []
             vtw_indices = []
-            for vtw_index, vtw in enumerate(self.requests[request_index][RequestIndices.vtw_list]):
+            vtw_list = self.requests[request_index][RequestIndices.vtw_list]
+            for vtw_index, vtw in enumerate(vtw_list):
                 conflict_degree.append(self._conflict_degree(vtw))
                 vtw_indices.append(vtw_index)
 
-            self.requests[request_index][RequestIndices.vtw_list] = [vtw for _, vtw in sorted(zip(conflict_degree, self.requests[request_index][RequestIndices.vtw_list]))]
+            self.requests[request_index][RequestIndices.vtw_list] = \
+                [vtw for _, vtw in sorted(zip(conflict_degree, vtw_list))]
 
     def sort_vtws_by_increasing_rise_time(self):
-
         for request in self.requests:
-            request[RequestIndices.vtw_list] = sorted(request[RequestIndices.vtw_list], key=lambda x: x.rise_time)
+            request[RequestIndices.vtw_list] = sorted(request[RequestIndices.vtw_list], key=lambda x: x.rise_time.data)
diff --git a/celest/schedule/schedule.py b/celest/schedule/schedule.py
index 72b5f7f..f67db93 100644
--- a/celest/schedule/schedule.py
+++ b/celest/schedule/schedule.py
@@ -1,40 +1,114 @@
 
 
-from celest.satellite.satellite import Satellite
+from celest.coordinates.frames.attitude import Attitude
+from celest.coordinates.ground_location import GroundLocation
+from celest.satellite import Satellite
 from celest.schedule.alns import ALNS
 from celest.schedule.insertion_operators import _INSERTION_FUNCTIONS
 from celest.schedule.removal_operators import _REMOVAL_FUNCTIONS
-from celest.schedule.request_handler import RequestHandler
+from celest.schedule.request_handler import RequestHandler, RequestIndices
 from celest.schedule.scheduling_utils import (
     cost,
     initialize_solution,
     is_complete
 )
-from celest.encounter.windows import generate_vtw
-from celest.encounter._window_handling import OW, OWHandler
+from celest.encounter.window_generator import generate_vtw, Lighting
+from celest.encounter.window_handling import ObservationWindow, WindowHandler
+from celest import units as u
 import math
+import numpy as np
 
 
 class Schedule(ALNS):
+    """Schedule(satellite, visibility_threshold)
 
-    def __init__(self, satellite, visibility_threshold):
+    Creates a scheduling object for the given satellite and visibility threshold.
 
+    This class implements an adaptive large neighborhood search metaheuristic
+    for scheduling satellite observations.
+
+    Parameters
+    ----------
+    satellite : Satellite
+        The satellite executing the scheduled observations.
+    vis_threshold : float
+        Visibility threshold in degrees.
+
+        The visibility threshold is the minimum elevation angle of the satellite
+        as seen from `location` where the satellite will be in visual range of
+        `location`.
+
+    Methods
+    -------
+    add_request(location, deadline, duration, priority, quality, look_ang, lighting)
+        Adds an encounter request to the schedule.
+    generate(max_iter, annealing_coeff, react_factor)
+        Generates a schedule for the given satellite.
+    """
+
+    def __init__(self, satellite: Satellite, vis_threshold: float) -> None:
         if not isinstance(satellite, Satellite):
             raise TypeError("satellite must be of type Satellite.")
-        if not isinstance(visibility_threshold, (float, int)):
-            raise TypeError("visibility_threshold must be a float or integer.")
+        if not isinstance(vis_threshold, (float, int)):
+            raise TypeError("vis_threshold must be a float or integer.")
 
         self.satellite = satellite
-        self.visibility_threshold = visibility_threshold
+        self.visibility_threshold = vis_threshold
         self.request_handler = RequestHandler()
 
-    def add_request(self, location, deadline, duration, priority, quality, look_ang):
-
-        vtws = generate_vtw(self.satellite, location, self.visibility_threshold, 1)
+    def add_request(self, location: GroundLocation, deadline: float, duration:
+                    float, priority: int, quality: int, look_ang: float=None,
+                    lighting: Lighting=Lighting.DAYTIME) -> None:
+        """Add a request to be scheduled.
+        
+        Parameters
+        ----------
+        location : GroundLocation
+            The location associated with the encounter.
+        deadline : float
+            The deadline for the request to be scheduled by in the jd2000 epoch.
+        duration : float
+            The desired duration of the encounter in seconds.
+        priority: int
+            The priority of the request on a scale of 1 to 10 where higher
+            priority requests are scheduled first.
+        quality : int
+            The quality of the request on a scale of 1 to 10. A higher quality
+            request will occur with a smaller look_angle. For imaging requests,
+            the imaging will be taken closer to nadir for a higher quality.
+        look_ang : float, optional
+            The desired look angle of the encounter in degrees, default is None.
+        lighting : Lighting, optional
+            The lighting condition of the encounter, default is DAYTIME.
+        """
+
+        vtws = generate_vtw(self.satellite, location, self.visibility_threshold, lighting)
         self.request_handler.add_request(location, deadline, duration, priority, quality, look_ang, vtws)
 
-    def generate(self, max_iter, annealing_coeff, react_factor):
-
+    def generate(self, max_iter: int, annealing_coeff: float, react_factor: float) -> None:
+        """Generate the schedule for the given satellite and requests.
+
+        Parameters
+        ----------
+        max_iter : int
+            The maximum number of iterations to perform in search for an
+            improved solution.
+        annealing_coeff : float
+            The annealing coefficient between 0 and 1 for the simulated
+            annealing process.
+        react_factor : float
+            Decay parameter within the range [0, 1].
+
+            The decay parameter dictates the effect of previous iterations on
+            the current iteration. A `react_factor` of 1 will cause the
+            algorithm to ignore previous iterations and a `react_factor` of 0
+            will cause the algorithm to only consider the previous iterations.
+
+        Returns
+        -------
+        WindowHandler
+            A container holding the scheduled windows.
+        """
         if not self.request_handler.number_of_requests:
             raise Exception("No requests to schedule.")
 
@@ -55,19 +129,31 @@ def generate(self, max_iter, annealing_coeff, react_factor):
         return self._generate_OWHandler_from_request_list(best_solution)
 
     def _generate_OWHandler_from_request_list(self, request_list):
-
-        ow_handler = OWHandler()
+        ow_handler = WindowHandler()
         for request in request_list:
 
-            if request[0]:
-
-                idx = request[1]
-                roll = request[10][idx].roll(request[2])
-                pitch = request[10][idx].pitch(request[2])
-                yaw = request[10][idx].yaw(request[2])
-
-                ow = OW(request[2], request[3],
-                        request[4], request[5], roll, pitch, yaw)
-                ow_handler._add_window(ow)
+            if request[RequestIndices.is_scheduled]:
+
+                idx = request[RequestIndices.vtw_index]
+                attitude = request[RequestIndices.vtw_list][idx].attitude
+                attitude_idx = np.where(
+                    attitude.time.to(u.jd2000).data ==
+                    request[RequestIndices.scheduled_start_time])[0]
+                time = attitude.time.data[attitude_idx]
+                roll = attitude.roll.data[attitude_idx]
+                pitch = attitude.pitch.data[attitude_idx]
+                yaw = attitude.yaw.data[attitude_idx]
+                unit = attitude.roll.unit
+                location = attitude.location
+                new_attitude = Attitude(time, roll, pitch, yaw, unit, location)
+
+                ow = ObservationWindow(
+                    request[RequestIndices.scheduled_start_time],
+                    request[RequestIndices.scheduled_duration],
+                    request[RequestIndices.deadline],
+                    request[RequestIndices.location],
+                    new_attitude
+                )
+                ow_handler.add_window(ow)
 
         return ow_handler
diff --git a/celest/units/core.py b/celest/units/core.py
index 6f9db32..53f2362 100644
--- a/celest/units/core.py
+++ b/celest/units/core.py
@@ -117,7 +117,6 @@ def insert_into_namespace(self, namespace):
             namespace[self._short_name] = self
 
 
-# TODO: Is there ever a need to not pass in a namespace?
 class Unit(NamedUnit, BaseUnit):
     """Unit(short_name, long_name, namespace=None, base_units=None)
 
diff --git a/celest/units/display.py b/celest/units/display.py
index 1407012..c08645f 100644
--- a/celest/units/display.py
+++ b/celest/units/display.py
@@ -3,9 +3,7 @@
 from celest.units import core  # Import core to prevent circular import errors.
 
 
-# TODO: There is a lot of repeated code here, try to refactor.
 def _to_string(unit):
-
     if isinstance(unit, core.CompoundUnit):
         return generate_compound_unit_string(unit)
     elif isinstance(unit, core.NamedUnit):
@@ -15,7 +13,6 @@ def _to_string(unit):
 
 
 def generate_compound_unit_string(unit):
-
     numerator_strings = []
     denominator_strings = []
     for base, power in zip(unit.bases, unit.powers):
diff --git a/celest/units/quantity.py b/celest/units/quantity.py
index ba93d66..c25429c 100644
--- a/celest/units/quantity.py
+++ b/celest/units/quantity.py
@@ -40,7 +40,8 @@ def __add__(self, other):
             return Quantity(data, self._unit)
         elif isinstance(other, Quantity):
             if repr(self._unit) != repr(other._unit):
-                raise ArithmeticError("Units must match for Quantity addition.")
+                raise ArithmeticError(
+                    "Units must match for Quantity addition.")
             data = self._data + other._data
             return Quantity(data, self._unit)
         else:
@@ -52,7 +53,8 @@ def __radd__(self, other):
             return Quantity(data, self._unit)
         elif isinstance(other, Quantity):
             if repr(self._unit) != repr(other._unit):
-                raise ArithmeticError("Units must match for Quantity addition.")
+                raise ArithmeticError(
+                    "Units must match for Quantity addition.")
             data = self._data + other._data
             return Quantity(data, self._unit)
         else:
@@ -64,7 +66,8 @@ def __sub__(self, other):
             return Quantity(data, self._unit)
         elif isinstance(other, Quantity):
             if repr(self._unit) != repr(other._unit):
-                raise ArithmeticError("Units must match for Quantity subtraction.")
+                raise ArithmeticError(
+                    "Units must match for Quantity subtraction.")
             data = self._data - other._data
             return Quantity(data, self._unit)
         else:
@@ -76,7 +79,8 @@ def __rsub__(self, other):
             return Quantity(data, self._unit)
         elif isinstance(other, Quantity):
             if repr(self._unit) != repr(other._unit):
-                raise ArithmeticError("Units must match for Quantity subtraction.")
+                raise ArithmeticError(
+                    "Units must match for Quantity subtraction.")
             data = other._data - self._data
             return Quantity(data, self._unit)
         else:
@@ -127,6 +131,54 @@ def __rtruediv__(self, other):
         else:
             return NotImplemented
 
+    def __eq__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data == other._data
+        else:
+            return NotImplemented
+
+    def __ne__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data != other._data
+        else:
+            return NotImplemented
+
+    def __lt__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data < other._data
+        else:
+            return NotImplemented
+
+    def __le__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data <= other._data
+        else:
+            return NotImplemented
+
+    def __gt__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data > other._data
+        else:
+            return NotImplemented
+
+    def __ge__(self, other):
+        if isinstance(other, Quantity):
+            if repr(self._unit) != repr(other._unit):
+                raise ArithmeticError("Units must match for comparison.")
+            return self._data >= other._data
+        else:
+            return NotImplemented
+
     @property
     def data(self):
         return self._data
@@ -182,8 +234,8 @@ def convert_to(self, new_unit):
         """
 
         if self._unit.dimension != new_unit.dimension:
-            raise ValueError("New unit has a different dimension then the \
-                             current unit.")
+            raise ValueError(
+                "New unit has a different dimension then the current unit.")
 
         def _get_unit_scale(unit):
             if isinstance(unit, CompoundUnit):
diff --git a/celest/units/unit_initialization.py b/celest/units/unit_initialization.py
index 1f8c6ef..88dc14d 100644
--- a/celest/units/unit_initialization.py
+++ b/celest/units/unit_initialization.py
@@ -8,7 +8,7 @@
 
 # Si length measures.
 setup_unit("m", "meter", namespace)
-setup_unit("mm", "milimeter", namespace, 0.001 * m)
+setup_unit("mm", "millimeter", namespace, 0.001 * m)
 setup_unit("cm", "centimeter", namespace, 0.01 * m)
 setup_unit("km", "kilometer", namespace, 1000 * m)
 


Test patch:
diff --git a/tests/test_coordinates/test_frames/test_attitude.py b/tests/test_coordinates/test_frames/test_attitude.py
new file mode 100644
index 0000000..972f692
--- /dev/null
+++ b/tests/test_coordinates/test_frames/test_attitude.py
@@ -0,0 +1,60 @@
+
+
+from celest.coordinates.frames.attitude import Attitude
+from celest.coordinates.ground_location import GroundLocation
+from celest import units as u
+from unittest import TestCase
+import numpy as np
+import os
+
+
+class TestAttitude(TestCase):
+
+    def setUp(self):
+        self.julian = np.random.rand(5)
+        self.roll = np.random.rand(5)
+        self.pitch = np.random.rand(5)
+        self.yaw = np.random.rand(5)
+        self.unit = u.deg
+        self.location = GroundLocation(52.1579, -106.6702, 0.482, u.deg, u.km)
+        self.valid_attitude = Attitude(self.julian, self.roll, self.pitch,
+                                       self.yaw, self.unit, self.location)
+
+    def test_initialization(self):
+        self.assertIsInstance(self.valid_attitude, Attitude)
+
+    def test_parent_class_is_initialized(self):
+        self.assertTrue(np.array_equal(self.julian, self.valid_attitude.time.data))
+        self.assertEqual(u.jd2000, self.valid_attitude.time.unit)
+
+        self.assertTrue(np.array_equal(self.roll, self.valid_attitude.roll.data))
+        self.assertEqual(self.unit, self.valid_attitude.roll.unit)
+
+        self.assertTrue(np.array_equal(self.pitch, self.valid_attitude.pitch.data))
+        self.assertEqual(self.unit, self.valid_attitude.pitch.unit)
+
+        self.assertTrue(np.array_equal(self.yaw, self.valid_attitude.yaw.data))
+        self.assertEqual(self.unit, self.valid_attitude.yaw.unit)
+
+    def test_error_raised_with_different_length_inputs(self):
+        julian = np.random.rand(1)
+        x = np.random.rand(2)
+        y = np.random.rand(3)
+        z = np.random.rand(4)
+        self.assertRaises(ValueError, Attitude, julian, x, y, z, self.unit,
+                          self.location)
+
+    def test_error_raised_when_wrong_dimension_input_passed_in(self):
+        julian = np.random.rand(1, 1)
+        x = np.random.rand(2, 2)
+        y = np.random.rand(3, 3)
+        z = np.random.rand(4, 4)
+        self.assertRaises(ValueError, Attitude, julian, x, y, z, self.unit,
+                          self.location)
+
+    def test_save_to_text_file_saves_file(self):
+        file_name = "attitude_test_file"
+        self.valid_attitude.save_text_file(file_name)
+        self.assertTrue(os.path.exists(file_name + ".txt"))
+        if os.path.exists(file_name + ".txt"):
+            os.remove(file_name + ".txt")
diff --git a/tests/test_coordinates/test_frames/test_azel.py b/tests/test_coordinates/test_frames/test_azel.py
index ea1b157..0365e5d 100644
--- a/tests/test_coordinates/test_frames/test_azel.py
+++ b/tests/test_coordinates/test_frames/test_azel.py
@@ -26,11 +26,11 @@ def test_parent_class_is_initialized(self):
         self.assertTrue(np.array_equal(self.julian, self.valid_azel.time.data))
         self.assertEqual(u.jd2000, self.valid_azel.time.unit)
 
-        self.assertTrue(np.array_equal(self.azimuth, self.valid_azel.x.data))
-        self.assertEqual(self.unit, self.valid_azel.x.unit)
+        self.assertTrue(np.array_equal(self.azimuth, self.valid_azel.azimuth.data))
+        self.assertEqual(self.unit, self.valid_azel.azimuth.unit)
 
-        self.assertTrue(np.array_equal(self.elevation, self.valid_azel.y.data))
-        self.assertEqual(self.unit, self.valid_azel.y.unit)
+        self.assertTrue(np.array_equal(self.elevation, self.valid_azel.elevation.data))
+        self.assertEqual(self.unit, self.valid_azel.elevation.unit)
 
         self.assertEqual(self.location, self.valid_azel.location)
 
diff --git a/tests/test_coordinates/test_frames/test_base_positions.py b/tests/test_coordinates/test_frames/test_base_positions.py
index a2e7b19..937475b 100644
--- a/tests/test_coordinates/test_frames/test_base_positions.py
+++ b/tests/test_coordinates/test_frames/test_base_positions.py
@@ -36,13 +36,13 @@ def test_initialization_with_non_arrays_of_different_lengths(self):
         self.assertRaises(ValueError, Position2d, x, self.x_unit, y,
                           self.y_unit, time, self.time_unit)
 
-    def test_x(self):
-        self.assertTrue(np.array_equal(self.x, self.position.x.data))
-        self.assertEqual(self.x_unit, self.position.x.unit)
+    def test_get_x(self):
+        self.assertTrue(np.array_equal(self.x, self.position._get_x().data))
+        self.assertEqual(self.x_unit, self.position._get_x().unit)
 
-    def test_y(self):
-        self.assertTrue(np.array_equal(self.y, self.position.y.data))
-        self.assertEqual(self.y_unit, self.position.y.unit)
+    def test_get_y(self):
+        self.assertTrue(np.array_equal(self.y, self.position._get_y().data))
+        self.assertEqual(self.y_unit, self.position._get_y().unit)
 
     def test_time(self):
         self.assertTrue(np.array_equal(self.time, self.position.time.data))
@@ -91,17 +91,17 @@ def test_initialization_with_non_arrays_of_different_lengths(self):
         self.assertRaises(ValueError, Position3d, x, self.x_unit, y,
                           self.y_unit, z, self.z_unit, time, self.time_unit)
 
-    def test_x(self):
-        self.assertTrue(np.array_equal(self.x, self.position.x.data))
-        self.assertEqual(self.x_unit, self.position.x.unit)
+    def test_get_x(self):
+        self.assertTrue(np.array_equal(self.x, self.position._get_x().data))
+        self.assertEqual(self.x_unit, self.position._get_x().unit)
 
-    def test_y(self):
-        self.assertTrue(np.array_equal(self.y, self.position.y.data))
-        self.assertEqual(self.y_unit, self.position.y.unit)
+    def test_get_y(self):
+        self.assertTrue(np.array_equal(self.y, self.position._get_y().data))
+        self.assertEqual(self.y_unit, self.position._get_y().unit)
 
-    def test_z(self):
-        self.assertTrue(np.array_equal(self.z, self.position.z.data))
-        self.assertEqual(self.z_unit, self.position.z.unit)
+    def test_get_z(self):
+        self.assertTrue(np.array_equal(self.z, self.position._get_z().data))
+        self.assertEqual(self.z_unit, self.position._get_z().unit)
 
     def test_time(self):
         self.assertTrue(np.array_equal(self.time, self.position.time.data))
diff --git a/tests/test_coordinates/test_ground_location.py b/tests/test_coordinates/test_ground_location.py
index e19e012..bf5bab5 100644
--- a/tests/test_coordinates/test_ground_location.py
+++ b/tests/test_coordinates/test_ground_location.py
@@ -20,6 +20,14 @@ def setUp(self):
     def test_initialization(self):
         self.assertIsInstance(self.location, GroundLocation)
 
+    def test_initialization_with_invalid_latitude_raises_value_error(self):
+        self.assertRaises(ValueError, GroundLocation, -91, 0, 0, u.deg, u.km)
+        self.assertRaises(ValueError, GroundLocation, 91, 0, 0, u.deg, u.km)
+
+    def test_initialization_with_invalid_longitude_raises_value_error(self):
+        self.assertRaises(ValueError, GroundLocation, 0, -181, 0, u.deg, u.km)
+        self.assertRaises(ValueError, GroundLocation, 0, 181, 0, u.deg, u.km)
+
     def test_str(self):
         self.assertEqual(str(self.location),
                          "00°05′39.84″S 22°36′05.40″W 493.42km")
diff --git a/tests/test_encounter/test_groundposition.py b/tests/test_encounter/test_groundposition.py
deleted file mode 100644
index 58f87b8..0000000
--- a/tests/test_encounter/test_groundposition.py
+++ /dev/null
@@ -1,49 +0,0 @@
-
-
-from celest.encounter.groundposition import GroundPosition
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-WGS84_MAJOR_AXIS = 6378.137
-WGS84_MINOR_AXIS = 6356.752314245
-
-
-class TestGroundPosition(TestCase):
-
-    def setUp(self):
-
-        self.latitude = 43.6532
-        self.longitude = -79.3832
-        self.elevation = 0.076
-
-    def test_radius(self):
-        """Test `GroundPosition._WGS84_radius`.
-
-        Notes
-        -----
-        Validation against online calculator methodology.[1]_
-
-        References
-        ----------
-        .. [1] Timur. Earth Radius by Latitude (WGS 84). 2018. url:
-           https://planetcalc.com/7721/.
-        """
-
-        location = GroundPosition(self.latitude, self.longitude, self.elevation)
-
-        cos_latitude = np.cos(np.radians(self.latitude))
-        sin_latitude = np.sin(np.radians(self.latitude))
-
-        numinator = (WGS84_MAJOR_AXIS ** 2 * cos_latitude) ** 2 + \
-            (WGS84_MINOR_AXIS ** 2 * sin_latitude) ** 2
-        denominator = (WGS84_MAJOR_AXIS * cos_latitude) ** 2 + \
-            (WGS84_MINOR_AXIS * sin_latitude) ** 2
-        radius = np.sqrt(numinator / denominator) + self.elevation
-
-        self.assertAlmostEqual(radius, location.radius, delta=0.001)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_encounter/test_window_generator.py b/tests/test_encounter/test_window_generator.py
new file mode 100644
index 0000000..ee53bec
--- /dev/null
+++ b/tests/test_encounter/test_window_generator.py
@@ -0,0 +1,142 @@
+
+
+from celest.coordinates.coordinate import Coordinate
+from celest.coordinates.frames.azel import AzEl
+from celest.coordinates.frames.gcrs import GCRS
+from celest.coordinates.ground_location import GroundLocation
+from celest.encounter.window_generator import (
+    generate_vtw,
+    _get_day_constraint_indices,
+    _get_night_constraint_indices,
+    Lighting,
+    _get_sun_gcrs
+)
+from celest.satellite import Satellite
+from celest import units as u
+from unittest import TestCase
+import numpy as np
+
+
+class TestGenerateVTW(TestCase):
+
+    def setUp(self):
+        data = np.loadtxt(
+            fname="tests/test_data/coordinate_validation_long.txt",
+            usecols=(0, 11, 12, 13, 14, 15, 16),
+            skiprows=1,
+            max_rows=5000
+        )
+
+        self.julian = data[:, 0] + 2430000
+        gcrs_position = data[:, 1:4]
+        gcrs_velocity = data[:, 4:]
+
+        self.gcrs_position = GCRS(
+            self.julian,
+            gcrs_position[:, 0],
+            gcrs_position[:, 1],
+            gcrs_position[:, 2],
+            u.km
+        )
+        self.gcrs_velocity = GCRS(
+            self.julian,
+            gcrs_velocity[:, 0],
+            gcrs_velocity[:, 1],
+            gcrs_velocity[:, 2],
+            u.m / u.s
+        )
+        self.satellite = Satellite(self.gcrs_position, self.gcrs_velocity)
+        self.location = GroundLocation(52.1579, -106.6702, 0.482, u.deg, u.km)
+
+    def test_generate_vtw_raises_warning_when_satellite_has_no_velocity(self):
+        self.assertRaises(Warning, generate_vtw, Satellite(self.gcrs_position),
+                          self.location, 10)
+
+    def test_generate_vtw_raises_value_error_for_negative_threshold(self):
+        self.assertRaises(ValueError, generate_vtw, self.satellite,
+                          self.location, -10)
+
+    def test_generate_vtw_raises_value_error_for_large_threshold(self):
+        self.assertRaises(ValueError, generate_vtw, self.satellite,
+                          self.location, 100)
+
+    def test_window_constraints_valid_at_day(self):
+        vis_threshold = 10
+        windows = generate_vtw(
+            self.satellite,
+            self.location,
+            vis_threshold,
+            Lighting.DAYTIME
+        )
+        satellite_elevation = self._get_satellite_elevation().to(u.deg).data
+        sun_elevation = self._get_sun_elevation().to(u.deg).data
+
+        for window in windows:
+            rise_time, set_time = window.rise_time.data, window.set_time.data
+            window_indices = np.where((rise_time <= self.julian) &
+                                      (self.julian <= set_time))[0]
+
+            # The first and last indices are always invalid.
+            self.assertTrue(np.alltrue(
+                satellite_elevation[window_indices][1:-1] > vis_threshold))
+            self.assertTrue(np.alltrue(sun_elevation[window_indices][1:-1] > 0))
+
+    def _get_satellite_elevation(self):
+        satellite_coordinates = Coordinate(self.satellite.position)
+        return satellite_coordinates.convert_to(AzEl, self.location).elevation
+
+    def _get_sun_elevation(self):
+        sun_gcrs = _get_sun_gcrs(self.julian)
+        sun_azel = Coordinate(sun_gcrs).convert_to(AzEl, self.location)
+        return sun_azel.elevation
+
+    def test_window_constraints_valid_at_night(self):
+        vis_threshold = 10
+        windows = generate_vtw(
+            self.satellite,
+            self.location,
+            vis_threshold,
+            Lighting.NIGHTTIME
+        )
+        satellite_elevation = self._get_satellite_elevation().to(u.deg).data
+        sun_elevation = self._get_sun_elevation().to(u.deg).data
+
+        for window in windows:
+            rise_time, set_time = window.rise_time.data, window.set_time.data
+            window_indices = np.where((rise_time <= self.julian) &
+                                      (self.julian <= set_time))[0]
+
+            # The first and last indices are always invalid.
+            self.assertTrue(np.alltrue(
+                satellite_elevation[window_indices][1:-1] > vis_threshold))
+            self.assertTrue(np.alltrue(sun_elevation[window_indices][1:-1] < 0))
+
+    def test_get_night_constraint_indices(self):
+        indices = _get_night_constraint_indices(self.julian, self.location)
+        sun_elevation = self._get_sun_elevation()
+        self.assertTrue(np.alltrue(sun_elevation.to(u.deg).data[indices] < 0))
+
+    def test_get_sun_gcrs(self):
+        from astropy.coordinates import get_body
+        from astropy import time, units
+
+        julian = np.array([2455368.75000, 2456293.54167, 2459450.85000])
+        astropy_time = time.Time(julian, format="jd")
+
+        sun_position = get_body("sun", astropy_time)
+        sun_position.representation_type = "cartesian"
+
+        expected_gcrs = sun_position.gcrs
+        actual_gcrs = _get_sun_gcrs(julian)
+
+        self.assertTrue(np.allclose(expected_gcrs.x.to(units.km).value,
+                                    actual_gcrs.x.to(u.km).data, rtol=0.05))
+        self.assertTrue(np.allclose(expected_gcrs.y.to(units.km).value,
+                                    actual_gcrs.y.to(u.km).data, rtol=0.05))
+        self.assertTrue(np.allclose(expected_gcrs.z.to(units.km).value,
+                                    actual_gcrs.z.to(u.km).data, rtol=0.05))
+
+    def test_get_day_constraint_indices(self):
+        sun_elevation = self._get_sun_elevation()
+        indices = _get_day_constraint_indices(self.julian, self.location)
+        self.assertTrue(np.alltrue(sun_elevation.to(u.deg).data[indices] > 0))
diff --git a/tests/test_encounter/test_window_handling.py b/tests/test_encounter/test_window_handling.py
index 3e4cfb5..6d3d0a9 100644
--- a/tests/test_encounter/test_window_handling.py
+++ b/tests/test_encounter/test_window_handling.py
@@ -1,177 +1,250 @@
 
 
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._window_handling import (
-    VTW,
-    OW,
-    WindowHandler,
-    VTWHandler,
-    OWHandler
+from celest.coordinates.frames.attitude import Attitude
+from celest.coordinates.ground_location import GroundLocation
+from celest.encounter.window_handling import (
+    ObservationWindow,
+    VisibleTimeWindow,
+    WindowHandler
 )
+from celest.units.quantity import Quantity
+from celest import units as u
 from unittest import TestCase
 import numpy as np
-import unittest
 import os
 
 
-class TestVTW(TestCase):
+class TestVisibleTimeWindow(TestCase):
 
-    def test_vtw_initialization(self):
-
-        rise_time = 21282.4
-        set_time = 21282.5
-        self.window = VTW(rise_time, set_time, None, None, None)
-
-
-class TestOW(TestCase):
-
-    def test_ow_initialization(self):
-
-        start_time = 21282.4
-        duration = 30
-        location = GroundPosition(43.6532, 79.3832, 0.076)
-        deadline = 21284.3
+    def setUp(self):
+        self.rise_time = 0
+        self.set_time = 1
+        self.attitude = Attitude(
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            u.deg,
+            GroundLocation(0, 0, 0, u.deg, u.km)
+        )
+        self.vtw = VisibleTimeWindow(
+            self.rise_time,
+            self.set_time,
+            self.attitude
+        )
+
+    def test_initialization(self):
+        self.assertIsInstance(self.vtw, VisibleTimeWindow)
+
+    def test_str(self):
+        self.assertEqual(f"Rise time: {self.rise_time} {u.jd2000}, Set time: "
+                         f"{self.set_time} {u.jd2000}, Attitude: "
+                         f"{self.attitude}", str(self.vtw))
+
+    def test_repr(self):
+        self.assertEqual(f"VisibleTimeWindow({self.rise_time}, "
+                         f"{self.set_time}, {self.attitude})", repr(self.vtw))
+
+    def test_rise_time_property(self):
+        self.assertEqual(self.rise_time, self.vtw.rise_time.data)
+
+    def test_set_time_property(self):
+        self.assertEqual(self.set_time, self.vtw.set_time.data)
+
+    def test_attitude_property(self):
+        self.assertEqual(self.attitude, self.vtw.attitude)
+
+
+class TestObservationWindow(TestCase):
 
-        self.window = OW(start_time, duration, location, deadline, None, None, None)
+    def setUp(self):
+        self.start_time = Quantity(0, u.jd2000)
+        self.duration = Quantity(1, u.s)
+        self.location = GroundLocation(0, 0, 0, u.deg, u.km)
+        self.attitude = Attitude(
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            u.deg,
+            self.location
+        )
+        self.obw = ObservationWindow(
+            self.start_time,
+            self.duration,
+            self.location,
+            self.attitude
+        )
+
+    def test_initialization(self):
+        self.assertIsInstance(self.obw, ObservationWindow)
+
+    def test_str(self):
+        self.assertEqual(f"Start time: {self.start_time}, Duration: "
+                         f"{self.duration}, Location: {self.location}, "
+                         f"Attitude: {self.attitude}", str(self.obw))
+
+    def test_repr(self):
+        self.assertEqual(f"ObservationWindow({self.start_time}, "
+                         f"{self.duration}, {self.location}, {self.attitude})",
+                         repr(self.obw))
+
+    def test_start_time_property(self):
+        self.assertEqual(self.start_time, self.obw.start_time)
+
+    def test_duration_property(self):
+        self.assertEqual(self.duration, self.obw.duration)
+
+    def test_location_property(self):
+        self.assertEqual(self.location, self.obw.location)
+
+    def test_attitude_property(self):
+        self.assertEqual(self.attitude, self.obw.attitude)
 
 
 class TestWindowHandler(TestCase):
 
     def setUp(self):
-
-        self.windows = WindowHandler()
-
-        self.W1 = VTW(1319, 1319.2, None, None, None)
-        self.W2 = VTW(1322.3, 1322.31, None, None, None)
-        self.W3 = VTW(1324, 1324.2, None, None, None)
-        self.W4 = VTW(1324.3, 1325.3, None, None, None)
-        self.W5 = VTW(1325, 1326.1, None, None, None)
-
-        self.windows._add_window_to_window_handler(self.W1)
-        self.windows._add_window_to_window_handler(self.W2)
-        self.windows._add_window_to_window_handler(self.W3)
-        self.windows._add_window_to_window_handler(self.W4)
-        self.windows._add_window_to_window_handler(self.W5)
-
-    def test_getitem(self):
-
-        self.assertEqual(self.windows[1310], self.W1)
-        self.assertEqual(self.windows[1319], self.W1)
-        self.assertEqual(self.windows[1323], self.W2)
-        self.assertEqual(self.windows[1324.4], self.W4)
-        self.assertEqual(self.windows[1330], self.W5)
-
-        test_windows = self.windows[1322.5, 1325]
-        true_windows = np.array([self.W2, self.W5], dtype=object)
-        self.assertTrue(np.array_equiv(test_windows, true_windows))
-
-        self.assertEqual(self.windows[1323.5, 1324.1], self.W3)
-
-        test_windows = self.windows[1322.5:1326]
-        true_windows = np.array([self.W3, self.W4, self.W5], dtype=object)
-        self.assertTrue(np.array_equiv(test_windows, true_windows))
-
-        test_windows = self.windows[1322.3:1325]
-        true_windows = np.array([self.W2, self.W3, self.W4, self.W5], dtype=object)
-        self.assertTrue(np.array_equiv(test_windows, true_windows))
-
-    def test_add_window(self):
-
-        test_val = 0
-        for w in self.windows:
-            window_start = w.rise_time if isinstance(w, VTW) else w.start_time
-            self.assertLessEqual(test_val, window_start)
-            test_val = window_start
-
-    def test_get_window(self):
-
-        self.assertEqual(self.windows.get_window(1320), self.W1)
-        self.assertEqual(self.windows.get_window(1322.5), self.W2)
-        self.assertEqual(self.windows.get_window(1324.5), self.W4)
-        self.assertEqual(self.windows.get_window(1325.5), self.W5)
-
-    def test_get_windows_in_range(self):
-
-        test_windows = self.windows.get_windows_in_range(1320, 1326)
-        true_windows = np.array([self.W2, self.W3, self.W4, self.W5], dtype=object)
-        self.assertTrue(np.array_equiv(test_windows, true_windows))
-
-    def test_to_numpy(self):
-
-        self.assertIsInstance(self.windows.to_numpy(), np.ndarray)
-
-
-class TestVTWHandler(TestCase):
-
-    def setUp(self) -> None:
-
-        self.windows = VTWHandler()
-
-        self.W1 = VTW(1319, 1319.2, None, None, None)
-        self.W2 = VTW(1322.3, 1322.31, None, None, None)
-        self.W3 = VTW(1324, 1324.2, None, None, None)
-        self.W4 = VTW(1324.3, 1325.3, None, None, None)
-        self.W5 = VTW(1325, 1326.1, None, None, None)
-
-        self.windows._add_window(self.W1)
-        self.windows._add_window(self.W2)
-        self.windows._add_window(self.W3)
-        self.windows._add_window(self.W4)
-        self.windows._add_window(self.W5)
-
-    def test_stats(self):
-
-        self.assertIsNotNone(self.windows.stats())
-
-    def test_save(self):
-
-        self.windows.save('test_vtw.csv', ',')
-        self.assertTrue(os.path.isfile('test_vtw.csv'))
-
-        loaded_data = np.loadtxt('test_vtw.csv', delimiter=',', skiprows=1)
-        for i, window in enumerate(self.windows):
-            self.assertEqual(window.rise_time, loaded_data[i, 0])
-            self.assertEqual(window.set_time, loaded_data[i, 1])
-            self.assertEqual(window.duration, loaded_data[i, 2])
-
-        os.remove('test_vtw.csv')
-
-
-class TestOWHandler(TestCase):
-
-    def setUp(self) -> None:
-
-        self.windows = OWHandler()
-        location = GroundPosition(43.6532, 79.3832, 0.076)
-
-        self.W1 = OW(1319, 30, location, 1319.2, None, None, None)
-        self.W2 = OW(1322.3, 20, location, 1322.31, None, None, None)
-        self.W3 = OW(1324, 32, location, 1324.2, None, None, None)
-        self.W4 = OW(1324.3, 28, location, 1325.3, None, None, None)
-        self.W5 = OW(1325, 10, location, 1326.1, None, None, None)
-
-        self.windows._add_window(self.W1)
-        self.windows._add_window(self.W2)
-        self.windows._add_window(self.W3)
-        self.windows._add_window(self.W4)
-        self.windows._add_window(self.W5)
-
-    def test_save(self):
-
-        self.windows.save('test_ow.csv', ',')
-        self.assertTrue(os.path.isfile('test_ow.csv'))
-
-        loaded_data = np.loadtxt('test_ow.csv', delimiter=',', skiprows=1)
-        for i, window in enumerate(self.windows):
-
-            self.assertEqual(window.location.latitude, loaded_data[i, 0])
-            self.assertEqual(window.location.longitude, loaded_data[i, 1])
-            self.assertEqual(window.start_time, loaded_data[i, 2])
-            self.assertEqual(window.duration, loaded_data[i, 3])
-            self.assertEqual(window.deadline, loaded_data[i, 4])
-
-        os.remove('test_ow.csv')
-
-
-if __name__ == "__main__":
-    unittest.main()
+        self.window_handler = WindowHandler()
+        self.ground_location = GroundLocation(0, 0, 0, u.deg, u.km)
+        self.vtw_attitude = Attitude(
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            np.random.rand(5),
+            u.deg,
+            self.ground_location
+        )
+        self.ow_attitude = Attitude(
+            np.random.rand(1),
+            np.random.rand(1),
+            np.random.rand(1),
+            np.random.rand(1),
+            u.deg,
+            self.ground_location
+        )
+        self.test_vtw = VisibleTimeWindow(0, 1, self.vtw_attitude)
+        self.test_obw = ObservationWindow(
+            Quantity(0, u.jd2000),
+            Quantity(1, u.s),
+            self.ground_location,
+            self.ow_attitude
+        )
+
+    def test_initialization(self):
+        self.assertIsInstance(self.window_handler, WindowHandler)
+        self.assertListEqual(self.window_handler._window_data, [])
+
+    def test_add_window_handler_data_raises_error_for_non_window(self):
+        self.assertRaises(TypeError, self.window_handler.add_window,
+                          "not a window")
+
+    def test_add_window_handler_data_alters_window_data_with_vtw(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.assertListEqual(self.window_handler._window_data, [self.test_vtw])
+
+    def test_add_window_handler_data_alters_window_data_with_ow(self):
+        self.window_handler.add_window(self.test_obw)
+        self.assertListEqual(self.window_handler._window_data, [self.test_obw])
+
+    def test_add_window_handler_data_enforces_same_window_type(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.assertRaises(TypeError, self.window_handler.add_window,
+                          self.test_obw)
+
+    def test_str(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.window_handler.add_window(self.test_vtw)
+        self.assertEqual(str([self.test_vtw, self.test_vtw]),
+                         str(self.window_handler))
+
+    def test_repr(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.window_handler.add_window(self.test_vtw)
+        self.assertEqual(f"WindowHandler({repr(self.test_vtw)}, {repr(self.test_vtw)})",
+                         repr(self.window_handler))
+
+    def test_len_returns_correct_length(self):
+        self.assertEqual(len(self.window_handler), 0)
+        self.window_handler.add_window(self.test_vtw)
+        self.assertEqual(len(self.window_handler), 1)
+        self.window_handler.add_window(self.test_vtw)
+        self.assertEqual(len(self.window_handler), 2)
+
+    def test_iteration(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.window_handler.add_window(self.test_vtw)
+        self.assertListEqual([window for window in self.window_handler],
+                             self.window_handler._window_data)
+
+    def test_indexing_with_integers(self):
+        test_ow_1 = ObservationWindow(Quantity(1, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_2 = ObservationWindow(Quantity(2, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_3 = ObservationWindow(Quantity(3, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_4 = ObservationWindow(Quantity(4, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_5 = ObservationWindow(Quantity(5, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+
+        self.window_handler.add_window(test_ow_1)
+        self.window_handler.add_window(test_ow_2)
+        self.window_handler.add_window(test_ow_3)
+        self.window_handler.add_window(test_ow_4)
+        self.window_handler.add_window(test_ow_5)
+
+        self.assertEqual(self.window_handler[0], test_ow_1)
+        self.assertEqual(self.window_handler[1], test_ow_2)
+        self.assertEqual(self.window_handler[2], test_ow_3)
+        self.assertEqual(self.window_handler[3], test_ow_4)
+        self.assertEqual(self.window_handler[4], test_ow_5)
+
+    def test_indexing_with_range(self):
+        test_ow_1 = ObservationWindow(Quantity(1, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_2 = ObservationWindow(Quantity(2, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_3 = ObservationWindow(Quantity(3, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_4 = ObservationWindow(Quantity(4, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+        test_ow_5 = ObservationWindow(Quantity(5, u.jd2000), Quantity(1, u.s),
+                                      self.ground_location, self.ow_attitude)
+
+        self.window_handler.add_window(test_ow_1)
+        self.window_handler.add_window(test_ow_2)
+        self.window_handler.add_window(test_ow_3)
+        self.window_handler.add_window(test_ow_4)
+        self.window_handler.add_window(test_ow_5)
+
+        self.assertListEqual(self.window_handler[:2], [test_ow_1, test_ow_2])
+        self.assertListEqual(self.window_handler[1:3], [test_ow_2, test_ow_3])
+        self.assertListEqual(self.window_handler[3:], [test_ow_4, test_ow_5])
+        self.assertListEqual(self.window_handler[:], [test_ow_1, test_ow_2,
+                                                      test_ow_3, test_ow_4,
+                                                      test_ow_5])
+
+    def test_save_text_file_raises_error_for_no_windows(self):
+        self.assertRaises(Exception, self.window_handler.save_text_file)
+
+    def test_save_text_file_for_vtws(self):
+        self.window_handler.add_window(self.test_vtw)
+        self.window_handler.add_window(self.test_vtw)
+
+        file_name = "vtw_test_file"
+        self.window_handler.save_text_file(file_name)
+        self.assertTrue(os.path.exists(file_name + ".txt"))
+        if os.path.exists(file_name + ".txt"):
+            os.remove(file_name + ".txt")
+
+    def test_save_text_file_for_ows(self):
+        self.window_handler.add_window(self.test_obw)
+        self.window_handler.add_window(self.test_obw)
+
+        file_name = "ow_test_file"
+        self.window_handler.save_text_file(file_name)
+        self.assertTrue(os.path.exists(file_name + ".txt"))
+        if os.path.exists(file_name + ".txt"):
+            os.remove(file_name + ".txt")
diff --git a/tests/test_encounter/test_windows.py b/tests/test_encounter/test_windows.py
deleted file mode 100644
index c914481..0000000
--- a/tests/test_encounter/test_windows.py
+++ /dev/null
@@ -1,93 +0,0 @@
-
-
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter.windows import _sun_coor, generate_vtw
-from celest.satellite.satellite import Satellite
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-class TestWindows(TestCase):
-
-    def setUp(self):
-
-        fname = "tests/test_data/coordinate_validation_long.txt"
-        cols = (0, 11, 12, 13, 14, 15, 16)
-        skiprows = 1
-        max_rows = 5000
-        data = np.loadtxt(fname=fname, usecols=cols, skiprows=skiprows,
-                          max_rows=max_rows)
-
-        times = data[:, 0]
-        gcrs_position = data[:, 1:4]
-        gcrs_velocity = data[:, 4:7]
-
-        offset = 2430000
-        self.satellite = Satellite(gcrs_position, gcrs_velocity, "gcrs", times,
-                                   offset)
-
-    def test_sun_coor(self):
-
-        from astropy.coordinates import get_body, ITRS
-        from astropy import time, units
-
-        julian = np.array([2455368.75000, 2456293.54167, 2459450.85000])
-        astropy_time = time.Time(julian, format="jd")
-
-        sun_position = get_body("sun", astropy_time)
-        sun_position.representation_type = "cartesian"
-
-        itrsCoor = sun_position.transform_to(ITRS(obstime=astropy_time))
-
-        true_x = itrsCoor.x.to(units.km).value
-        true_y = itrsCoor.y.to(units.km).value
-        true_z = itrsCoor.z.to(units.km).value
-
-        test_x, test_y, test_z, _, _, _ = _sun_coor(julian).itrs(False)
-        self.assertTrue(np.allclose(true_x, test_x, rtol=0.05))
-        self.assertTrue(np.allclose(true_y, test_y, rtol=0.05))
-        self.assertTrue(np.allclose(true_z, test_z, rtol=0.05))
-
-    def test_windows(self):
-
-        location = GroundPosition(43.6532, -79.3832)
-        elevation, _ = self.satellite.horizontal(location, True)
-
-        sun_coordinates = _sun_coor(self.satellite._julian)
-        sun_elevation, _ = sun_coordinates.horizontal(location, True)
-
-        vis_threshold, lighting, tol = 10, 0, 1e-5
-        windows = generate_vtw(self.satellite, location, vis_threshold, lighting, tol)
-
-        for window in windows:
-            rise_time, set_time = window.rise_time, window.set_time
-            window_times = np.linspace(rise_time, set_time, 10)
-
-            self.assertTrue(np.all(elevation(window_times) > vis_threshold))
-
-        vis_threshold, lighting, tol = 10, 1, 1e-5
-        windows = generate_vtw(self.satellite, location, vis_threshold, lighting, tol)
-
-        for window in windows:
-            rise_time, set_time = window.rise_time, window.set_time
-            window_times = np.linspace(rise_time, set_time, 10)
-            sun_elev = sun_elevation(window_times)
-
-            self.assertTrue(np.all(elevation(window_times) > vis_threshold))
-            self.assertTrue(np.all((0 < sun_elev)))
-
-        vis_threshold, lighting, tol = 10, -1, 1e-5
-        windows = generate_vtw(self.satellite, location, vis_threshold, lighting, tol)
-
-        for window in windows:
-            rise_time, set_time = window.rise_time, window.set_time
-            window_times = np.linspace(rise_time, set_time, 10)
-            sun_elev = sun_elevation(window_times)
-
-            self.assertTrue(np.all(elevation(window_times) > vis_threshold))
-            self.assertTrue(np.all((sun_elev <= 0)))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_file_save.py b/tests/test_file_save.py
index 8c07932..817fd34 100644
--- a/tests/test_file_save.py
+++ b/tests/test_file_save.py
@@ -1,27 +1,58 @@
 
 
-from celest.file_save import get_parameter_string, _save_data_as_txt
+from celest.file_save import TextFileWriter, _process_file_name
 from unittest import TestCase
+import numpy as np
 import os
 
 
-class TestSaveDataAsTxt(TestCase):
+class TestTextFileWriter(TestCase):
 
     def setUp(self):
-        self.file_name = "test_file"
+        self.file_name = "file_name"
+        self.processed_file_name = _process_file_name(self.file_name)
+        self.header = "header"
+        self.writer = TextFileWriter(self.file_name, self.header)
+
+        self.subheading = "subheading"
+        self.parameters = [
+            ["label 1", "value 1"],
+            ["label 2", "value 2"]
+        ]
+        self.data = [
+            ["data 1", np.random.rand(5)],
+            ["data 2", np.random.rand(5)]
+        ]
+        self.data_format = "%.5f"
 
     def tearDown(self):
-        if os.path.exists(self.file_name + ".txt"):
-            os.remove(self.file_name + ".txt")
-
-    def test_save_data_creates_file(self):
-        _save_data_as_txt(self.file_name)
-        self.assertTrue(os.path.exists(self.file_name + ".txt"))
-
-    def test_get_parameter_string_for_non_final_parameter(self):
-        parameter_string = get_parameter_string(0, 2, "label", "value")
-        self.assertEqual("label: value\n", parameter_string)
-
-    def test_get_parameter_string_for_final_parameter(self):
-        parameter_string = get_parameter_string(0, 1, "label", "value")
-        self.assertEqual("label: value", parameter_string)
+        if os.path.exists(self.processed_file_name):
+            os.remove(self.processed_file_name)
+
+    def test_initialize(self):
+        self.assertIsInstance(self.writer, TextFileWriter)
+
+    def test_process_file_name_ends_in_txt_when_not_in_file_name(self):
+        file_name = "test_file"
+        self.assertEqual(file_name + ".txt", _process_file_name(file_name))
+
+    def test_process_file_name_ends_in_txt_when_in_file_name(self):
+        file_name = "test_file.txt"
+        self.assertEqual(file_name, _process_file_name(file_name))
+
+    def test_adding_layer_adds_layer(self):
+        self.assertEqual(0, len(self.writer._layers))
+        self._add_test_layer()
+        self.assertEqual(1, len(self.writer._layers))
+        self._add_test_layer()
+        self.assertEqual(2, len(self.writer._layers))
+
+    def _add_test_layer(self):
+        self.writer.add_layer(self.subheading, self.parameters, self.data,
+                              self.data_format)
+
+    def test_save_creates_file(self):
+        self._add_test_layer()
+        self._add_test_layer()
+        self.writer.save()
+        self.assertTrue(os.path.exists(self.processed_file_name))
diff --git a/tests/test_satellite.py b/tests/test_satellite.py
index 2c36672..f985e8d 100644
--- a/tests/test_satellite.py
+++ b/tests/test_satellite.py
@@ -1,7 +1,10 @@
 
 
-from celest.coordinates import WGS84, ITRS, GroundLocation
+from celest.coordinates.frames.attitude import Attitude
 from celest.coordinates.frames.gcrs import GCRS
+from celest.coordinates.frames.itrs import ITRS
+from celest.coordinates.frames.wgs84 import WGS84
+from celest.coordinates.ground_location import GroundLocation
 from celest.satellite import Satellite
 from celest import units as u
 from unittest import TestCase
@@ -12,25 +15,34 @@
 class TestSatellite(TestCase):
 
     def setUp(self):
-        fname = "tests/test_data/coordinate_validation_long.txt"
-        cols = (0, 5, 6, 7, 11, 12, 13, 14, 15, 16, 17, 18, 19)
-        skiprows = 1
-        max_rows = 5000
-        data = np.loadtxt(fname=fname, usecols=cols, skiprows=skiprows,
-                          max_rows=max_rows)
-
-        self.times = data[:, 0] + 2430000
+        data = np.loadtxt(
+            fname="tests/test_data/coordinate_validation_long.txt",
+            usecols=(0, 5, 6, 7, 11, 12, 13, 14, 15, 16, 17, 18, 19),
+            skiprows=1,
+            max_rows=5000
+        )
+
+        self.julian = data[:, 0] + 2430000
         self.wgs84 = data[:, 1:3]
         self.altitude = data[:, 3]
         self.gcrs = data[:, 4:7]
         self.gcrs_velocity = data[:, 7:10]
         self.itrs = data[:, 10:]
 
-        self.gcrs_position = GCRS(self.times, self.gcrs[:, 0], self.gcrs[:, 1],
-                                  self.gcrs[:, 2], u.km)
-        self.gcrs_velocity = GCRS(self.times, self.gcrs_velocity[:, 0],
-                                  self.gcrs_velocity[:, 1],
-                                  self.gcrs_velocity[:, 2], u.m / u.s)
+        self.gcrs_position = GCRS(
+            self.julian,
+            self.gcrs[:, 0],
+            self.gcrs[:, 1],
+            self.gcrs[:, 2],
+            u.km
+        )
+        self.gcrs_velocity = GCRS(
+            self.julian,
+            self.gcrs_velocity[:, 0],
+            self.gcrs_velocity[:, 1],
+            self.gcrs_velocity[:, 2],
+            u.m / u.s
+        )
 
         self.satellite = Satellite(self.gcrs_position, self.gcrs_velocity)
 
@@ -40,13 +52,19 @@ def test_initialization(self):
         self.assertIsInstance(self.satellite, Satellite)
 
     def test_value_error_raised_from_improper_data_frames(self):
-        wgs84_position = WGS84(self.times, self.wgs84[:, 0], self.wgs84[:, 1],
-                               self.altitude, u.deg, u.km)
+        wgs84_position = WGS84(
+            self.julian,
+            self.wgs84[:, 0],
+            self.wgs84[:, 1],
+            self.altitude,
+            u.deg,
+            u.km
+        )
         self.assertRaises(ValueError, Satellite, wgs84_position)
 
     def test_attitude(self):
         self.satellite.attitude(self.location)
-        self.assertIsNotNone(self.satellite.attitude(self.location))
+        self.assertIsInstance(self.satellite.attitude(self.location), Attitude)
 
     def test_look_angle(self):
         """
@@ -59,8 +77,13 @@ def test_look_angle(self):
         generate the test cases.
         """
 
-        gcrs = GCRS(self.times[210:220], self.gcrs[210:220, 0],
-                    self.gcrs[210:220, 1], self.gcrs[210:220, 2], u.km)
+        gcrs = GCRS(
+            self.julian[210:220],
+            self.gcrs[210:220, 0],
+            self.gcrs[210:220, 1],
+            self.gcrs[210:220, 2],
+            u.km
+        )
         satellite = Satellite(gcrs)
 
         expected_look_angle = np.array([66.88, 65.09, 63.90, 63.22, 62.46,
diff --git a/tests/test_satellite_old/__init__.py b/tests/test_satellite_old/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/tests/test_satellite_old/test_angle_representations.py b/tests/test_satellite_old/test_angle_representations.py
deleted file mode 100644
index 96a5c4f..0000000
--- a/tests/test_satellite_old/test_angle_representations.py
+++ /dev/null
@@ -1,35 +0,0 @@
-
-
-from celest.satellite._angle_representations import (
-    sexagesimal,
-    _ISO6709_representation
-)
-from unittest import TestCase
-import numpy as np
-
-
-class TestAngleRepresentations(TestCase):
-
-    def test_sexagesimal(self):
-
-        angles = np.array([43.6532, -79.3832, -33.2833, 149.1000])
-
-        true_sexagesimal = ["+43\u00B039\u203211.52\u2033",
-                            "-79\u00B022\u203259.52\u2033",
-                            "-33\u00B016\u203259.88\u2033",
-                            "+149\u00B006\u203200.00\u2033"]
-        calc_sexagesimal = sexagesimal(angles)
-
-        self.assertTrue(np.array_equal(true_sexagesimal, calc_sexagesimal))
-
-    def test_ISO6709_representation(self):
-
-        latitude = np.array([43.6532, -33.2833])
-        longitude = np.array([-79.3832, 149.1000])
-        altitude = np.array([430.23, 532.98])
-
-        true_iso_format = ["43\u00B039\u203211.52\u2033N 79\u00B022\u203259.52\u2033W 430.23km",
-                           "33\u00B016\u203259.88\u2033S 149\u00B006\u203200.00\u2033E 532.98km"]
-        calc_iso_format = _ISO6709_representation(latitude, longitude, altitude)
-
-        self.assertTrue(np.array_equal(true_iso_format, calc_iso_format))
diff --git a/tests/test_satellite_old/test_astronomical_quantities.py b/tests/test_satellite_old/test_astronomical_quantities.py
deleted file mode 100644
index 0b84871..0000000
--- a/tests/test_satellite_old/test_astronomical_quantities.py
+++ /dev/null
@@ -1,225 +0,0 @@
-
-
-from celest.satellite._astronomical_quantities import *
-from unittest import TestCase
-import julian as jd
-import numpy as np
-import unittest
-
-
-class TestAstronomicalQuantities(TestCase):
-
-    def test_nutation_angles(self):
-        """Test `AstronomicalQuantities.nutation_angles`.
-
-        Notes
-        -----
-        Test case is taken from Astronomical Algorithms.[1]_
-
-        References
-        ----------
-        .. [1] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 148. isbn: 9780943396613.
-        """
-
-        julian = np.array([2446895.5])
-        D, M, N, F, O = nutation_angles(julian)
-
-        self.assertAlmostEqual(D[0], 136.9623, places=4)
-        self.assertAlmostEqual(M[0], 94.9792, places=4)
-        self.assertAlmostEqual(N[0], 229.2784, places=4)
-        self.assertAlmostEqual(F[0], 143.4079, places=4)
-        self.assertAlmostEqual(O[0], 11.2531, places=4)
-
-    def test_nutation_components(self):
-        """Test `AstronomicalQuantities.nutation_components`.
-
-        Notes
-        -----
-        Test cases are taken from the PHP Science Labs.[1]_
-
-        References
-        ----------
-        .. [1] Jay Tanner. NeoProgrammics - Science Computations. 2021.
-           url: http://www.neoprogrammics.com/nutations/.
-        """
-
-        julian = np.array([2449634.50000, 2453420.56250, 2477418.21181])
-
-        true_lon_nutation = np.array([11.694, -5.993, 2.937])
-        true_obliquity_nutation = np.array([-5.946, 8.431, -8.871])
-        test_lon_nutation, test_obliquity_nutation = nutation_components(julian)
-
-        self.assertTrue(np.allclose(test_lon_nutation, true_lon_nutation,
-                        atol=0.5))
-        self.assertTrue(np.allclose(test_obliquity_nutation,
-                        true_obliquity_nutation, atol=0.1))
-
-    def test_mean_obliquity(self):
-        """Test `AstronomicalQuantities.mean_obliquity`.
-
-        Notes
-        -----
-        Test cases are taken from PHP Science Labs.[1]_
-
-        References
-        ----------
-        .. [1] Jay Tanner. Obliquity of the Ecliptic - PHP Science Labs. 2021.
-           url: https://www.neoprogrammics.com/obliquity_of_the_ecliptic/Obliq
-           uity_Of_The_Ecliptic_Calculator.php
-        """
-
-        julian = np.array([2459437.81600, 2477404.57292, 2422327.21875])
-
-        true_mean_obliquity = np.array([23.43648, 23.43008, 23.44969])
-        test_mean_obliquity = mean_obliquity(julian)
-
-        self.assertTrue(np.allclose(test_mean_obliquity, true_mean_obliquity,
-                        atol=0.0001))
-
-    def test_apparent_obliquity(self):
-        """Test `AstronomicalQuantities.apparent_obliquity`.
-
-        Notes
-        -----
-        Test cases are taken from PHP Science Labs.[1]_
-
-        References
-        ----------
-        .. [1] Jay Tanner. Obliquity of the Ecliptic - PHP Science Labs. 2021.
-           url: https://www.neoprogrammics.com/obliquity_of_the_ecliptic/Obliq
-           uity_Of_The_Ecliptic_Calculator.php
-        """
-
-        julian = np.array([2459437.81597, 2477404.57292, 2422327.21875])
-
-        true_apparent_obliquity = np.array([23.43763, 23.42763, 23.44798])
-        calc_apparent_obliquity = apparent_obliquity(julian)
-
-        self.assertTrue(np.allclose(calc_apparent_obliquity,
-                        true_apparent_obliquity, atol=0.0001))
-
-    def test_from_julian(self):
-        """Test `AstronomicalQuantities.from_julian`.
-
-        Notes
-        -----
-        Test cases are generated from the `Julian` Python library.
-        """
-
-        julian = np.array([2436116.31000, 2445246.65000, 2456124.09000])
-        test_year, test_month, test_day = from_julian(julian)
-
-        true_year, true_month, true_day = [], [], []
-
-        for current_julian in julian:
-            dt = jd.from_jd(current_julian)
-
-            seconds = dt.second + dt.microsecond / 1e6
-            minutes = dt.minute + seconds / 60
-            hours = dt.hour + minutes / 60
-            day = dt.day + hours / 24
-            month = dt.month
-            year = dt.year
-
-            true_year.append(year)
-            true_month.append(month)
-            true_day.append(day)
-
-        self.assertTrue(np.array_equal(test_year, true_year))
-        self.assertTrue(np.array_equal(test_month, true_month))
-        self.assertTrue(np.allclose(test_day, true_day, atol=0.0001))
-
-    def test_day_of_year(self):
-        """Test `AstronomicalQuantities.day_of_year`.
-
-        Notes
-        -----
-        Test cases are taken from "Astronomical Algorithms" by Jean Meeus.[1]_
-
-        References
-        ----------
-        .. [1] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 65. isbn: 9780943396613.
-        """
-
-        julian = np.array([2443826.5, 2447273.5, 2447273.8, 2447274.4])
-
-        true_day = np.array([318, 113, 113, 113])
-        test_day = day_of_year(julian)
-
-        self.assertTrue(np.array_equal(test_day, true_day))
-
-    def test_equation_of_time(self):
-        """Test `AstronomicalQuantities.equation_of_time`.
-
-        Notes
-        -----
-        Test cases are taken from "Astronomical Algorithms" by Jean Meeus,
-        PLANETCALC, and the Global Monitoring Laboratory.[1]_[2]_[3]_
-
-        References
-        ----------
-        .. [1] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 185. isbn: 9780943396613.
-        .. [2] Anton. Online calculator: Equation of time. url:
-           https://planetcalc.com/9235/.
-        .. [3] NOAA US Department of Commerce. ESRL Global Monitoring
-           Laboratory -Global Radiation and Aerosols. url:
-           https://gml.noaa.gov/grad/solcalc/.
-        """
-
-        julian = np.array([2455368.75, 2448908.50, 2459448.50])
-
-        true_equation_of_time = np.array([-0.42658, 3.42735, -0.71054])
-        test_equation_of_time = equation_of_time(julian)
-
-        self.assertTrue(np.allclose(test_equation_of_time,
-                        true_equation_of_time, atol=0.04))
-
-    def test_equation_of_equinoxes(self):
-        """Test `AstronomicalQuantities.equation_of_equinoxes`.
-
-        Notes
-        -----
-        Test cases taken from "Astronomical Algorithms" by Jean Meeus.[1]_
-
-        References
-        ----------
-        .. [1] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 88. isbn: 9780943396613.
-        """
-
-        julian = np.array([2446895.5])
-
-        true_equation_of_equinoxes = np.array([-0.23170])
-        test_equation_of_equinoxes = equation_of_equinoxes(julian)
-
-        self.assertTrue(np.allclose(test_equation_of_equinoxes,
-                        true_equation_of_equinoxes, atol=0.005))
-
-    def test_sun_right_ascension(self):
-        """Test `AstronomicalQuantities.sun_right_ascension`.
-
-        Notes
-        -----
-        The test case is taken from "Astronomical Algorithms" by Jean
-        Meeus.[1]_
-
-        References
-        ----------
-        .. [1] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-           1998, pp. 185. isbn: 9780943396613.
-        """
-
-        julian = np.array([2448908.5])
-
-        true_sun_right_ascension = np.array([198.38083])
-        test_sun_right_ascension = sun_right_ascension(julian)
-
-        self.assertTrue(np.allclose(true_sun_right_ascension,
-                        test_sun_right_ascension, atol=0.001))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_satellite_old/test_coordinate.py b/tests/test_satellite_old/test_coordinate.py
deleted file mode 100644
index f517664..0000000
--- a/tests/test_satellite_old/test_coordinate.py
+++ /dev/null
@@ -1,255 +0,0 @@
-
-
-from celest.encounter.groundposition import GroundPosition
-from celest.satellite.coordinate import Coordinate
-from polare import Stroke
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-WGS84_MAJOR_AXIS = 6378.137
-WGS84_MINOR_AXIS = 6356.752314245
-
-
-class TestCoordinate(TestCase):
-
-    def setUp(self):
-
-        fname = "tests/test_data/coordinate_validation_long.txt"
-        cols = (0, 5, 6, 7, 11, 12, 13, 14, 15, 16, 17, 18, 19)
-        skiprows = 1
-        max_rows = 5000
-        data = np.loadtxt(fname=fname, usecols=cols, skiprows=skiprows,
-                          max_rows=max_rows)
-
-        self.times = data[:, 0]
-        self.geographical = data[:, 1:3]
-        self.altitude = data[:, 3]
-        self.GCRS = data[:, 4:7]
-        self.GCRS_vel = data[:, 7:10]
-        self.ITRS = data[:, 10:]
-        self._length = data.shape[0]
-
-        self.tri_geo = np.concatenate((self.geographical,
-                                       self.altitude.reshape((-1, 1))), axis=1)
-
-        self.offset = 2430000
-        self.coor_gcrs = Coordinate(self.GCRS, self.GCRS_vel, "gcrs",
-                                    self.times, self.offset)
-        self.coor_itrs = Coordinate(self.ITRS, self.GCRS_vel, "itrs",
-                                    self.times, self.offset)
-
-    def test_set_base_data(self):
-
-        self.assertIsNotNone(self.coor_gcrs._julian)
-        self.assertIsNotNone(self.coor_gcrs._GCRS)
-        self.assertIsNotNone(self.coor_gcrs._ITRS)
-        self.assertEqual(self._length, self.coor_gcrs._length)
-
-        self.assertIsNotNone(self.coor_itrs._julian)
-        self.assertIsNotNone(self.coor_itrs._GCRS)
-        self.assertIsNotNone(self.coor_itrs._ITRS)
-        self.assertEqual(self._length, self.coor_itrs._length)
-
-    def test_geo_to_itrs(self):
-
-        test_itrs = Coordinate._geo_to_itrs(self, self.tri_geo)
-
-        self.assertTrue(np.allclose(test_itrs[:, 0], self.ITRS[:, 0],
-                                    atol=0.001))
-        self.assertTrue(np.allclose(test_itrs[:, 1], self.ITRS[:, 1],
-                                    atol=0.001))
-        self.assertTrue(np.allclose(test_itrs[:, 2], self.ITRS[:, 2],
-                                    atol=0.001))
-
-    def test_itrs_to_geo(self):
-
-        itrs = np.array([Stroke(self.times, self.ITRS[:, i]) for i in range(3)])
-        test_lat, test_lon = Coordinate._itrs_to_geo(self, itrs)
-
-        self.assertTrue(np.allclose(test_lat(self.times),
-                        self.geographical[:, 0], atol=0.18))
-        self.assertTrue(np.allclose(test_lon(self.times),
-                        self.geographical[:, 1], atol=0.00001))
-
-    def test_geo(self):
-
-        lat1, lon1, alt1 = self.coor_gcrs.geo()
-        lat2, lon2, alt2 = self.coor_itrs.geo()
-
-        for i in range(self._length):
-
-            eps = 1
-            calc_lon, test_lon = lon1[i], self.geographical[i, 1]
-            cond_one = abs(180 - calc_lon) < eps and abs(180 + test_lon) < eps
-            cond_two = abs(180 + calc_lon) < eps and abs(180 - test_lon) < eps
-
-            if cond_one:
-                test_lon = (test_lon + 360) % 360
-            elif cond_two:
-                calc_lon = (calc_lon + 360) % 360
-
-            self.assertAlmostEqual(lat1[i], self.geographical[i, 0], delta=0.5)
-            self.assertAlmostEqual(calc_lon, test_lon, delta=0.9)
-            self.assertAlmostEqual(alt1[i], self.altitude[i], delta=0.1)
-
-        self.assertTrue(np.allclose(lat2, self.geographical[:, 0], atol=0.18))
-        self.assertTrue(np.allclose(lon2, self.geographical[:, 1], atol=0.00001))
-        self.assertTrue(np.allclose(alt2, self.altitude, atol=0.06))
-
-    def test_era(self):
-
-        julian = [2454545]
-        position = [[6343.82, -2640.87, -11.26]]
-        velocity = [[0, 0, 0]]
-        coor = Coordinate(position, velocity, "itrs", julian)
-
-        true_era = np.degrees(np.array([6.2360075]))
-        test_era = coor.era()
-
-        self.assertAlmostEqual(true_era[0], test_era[0], delta=0.01)
-
-    def test_gcrs_and_itrs(self):
-
-        calc_itrs = self.coor_gcrs._gcrs_and_itrs(self.GCRS, frame="gcrs")
-        calc_gcrs = self.coor_itrs._gcrs_and_itrs(self.ITRS, frame="itrs")
-
-        self.assertTrue(np.allclose(calc_itrs[:, 0], self.ITRS[:, 0], atol=0.35))
-        self.assertTrue(np.allclose(calc_itrs[:, 1], self.ITRS[:, 1], atol=0.35))
-        self.assertTrue(np.allclose(calc_itrs[:, 2], self.ITRS[:, 2], atol=0.35))
-
-        self.assertTrue(np.allclose(calc_gcrs[:, 0], self.GCRS[:, 0], atol=0.35))
-        self.assertTrue(np.allclose(calc_gcrs[:, 1], self.GCRS[:, 1], atol=0.35))
-        self.assertTrue(np.allclose(calc_gcrs[:, 2], self.GCRS[:, 2], atol=0.35))
-
-    def test_gcrs(self):
-
-        calc_gcrs_x, calc_gcrs_y, calc_gcrs_z, _, _, _ = self.coor_itrs.gcrs()
-
-        self.assertTrue(np.allclose(calc_gcrs_x, self.GCRS[:, 0], atol=0.35))
-        self.assertTrue(np.allclose(calc_gcrs_y, self.GCRS[:, 1], atol=0.35))
-        self.assertTrue(np.allclose(calc_gcrs_z, self.GCRS[:, 2], atol=0.35))
-
-    def test_itrs(self):
-
-        calc_itrs_x, calc_itrs_y, calc_itrs_z, _, _, _ = self.coor_gcrs.itrs()
-
-        self.assertTrue(np.allclose(calc_itrs_x, self.ITRS[:, 0], atol=0.35))
-        self.assertTrue(np.allclose(calc_itrs_y, self.ITRS[:, 1], atol=0.35))
-        self.assertTrue(np.allclose(calc_itrs_z, self.ITRS[:, 2], atol=0.35))
-
-    def test_lvlh(self):
-
-        lvlh_data = self.coor_gcrs.lvlh()
-        self.assertIsNotNone(lvlh_data)
-
-    def test_get_ang(self):
-
-        vector_one = np.array([[56, 92, 76], [9238, 8479, 9387], [2, 98, 23]])
-        vector_two = np.array([[36, 29, 38], [2703, 947, 8739], [9827, 921, 1]])
-
-        true_angle = np.array([16.28, 37, 83.65])
-        test_angle = Coordinate._get_ang(self, vector_one, vector_two)
-
-        self.assertTrue(np.allclose(test_angle, true_angle, atol=0.01))
-
-    def test_horizontal(self):
-
-        from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz
-        from astropy import units as u
-        from astropy import time
-
-        latitude, longitude, height = 52.1579, -106.6702, 0.482
-        location = EarthLocation.from_geodetic(longitude * u.deg,
-                                               latitude * u.deg,
-                                               height * u.km)
-
-        times = time.Time(self.times + self.offset, format="jd")
-        x, y, z = self.GCRS.T
-
-        gcrs = SkyCoord(x=x, y=y, z=z, unit='km', frame=GCRS(obstime=times),
-                        representation_type='cartesian')
-        el_az = gcrs.transform_to(AltAz(obstime=times, location=location))
-
-        true_elevation = el_az.alt.degree
-        true_azimuth = el_az.az.degree
-
-        c = Coordinate(self.GCRS, self.GCRS_vel, "gcrs", self.times, self.offset)
-        location = GroundPosition(latitude, longitude, height)
-        test_elevation, test_azimuth = c.horizontal(location)
-
-        self.assertTrue(np.allclose(true_elevation, test_elevation, atol=0.32))
-        self.assertTrue(np.allclose(true_azimuth, test_azimuth, atol=7.3))
-
-    def test_off_nadir(self):
-        """
-        Notes
-        -----
-        Test cases generated using a geometric model designed by Mingde Yin.
-
-        The tolerance required to pass all test cases is artificially inflated
-        due to the low precision in the model's output and parameters used to
-        generate the test cases.
-        """
-
-        location = GroundPosition(52.1579, -106.6702, 0.482)
-        coor = Coordinate(self.GCRS[210:220], self.GCRS_vel[210:220], "gcrs",
-                          self.times[210:220], self.offset)
-
-        true_off_nadir = np.array([66.88, 65.09, 63.90, 63.22, 62.46, 61.67,
-                                   58.42, 38.27, 23.73, 56.29])
-        test_off_nadir = coor.off_nadir(location)
-
-        self.assertTrue(np.allclose(test_off_nadir, true_off_nadir, atol=0.25))
-
-    def test_WGS84_radius(self):
-        """
-        Notes
-        -----
-        Validation against online calculator methodology.[1]_
-
-        References
-        ----------
-        .. [1] Timur. Earth Radius by Latitude (WGS 84). 2018. url:
-           https://planetcalc.com/7721/.
-        """
-
-        latitude = np.radians(self.geographical[:, 0])
-        numerator = (WGS84_MAJOR_AXIS ** 2 * np.cos(latitude)) ** 2 + \
-            (WGS84_MINOR_AXIS ** 2 * np.sin(latitude)) ** 2
-        denominator = (WGS84_MAJOR_AXIS * np.cos(latitude)) ** 2 + \
-            (WGS84_MINOR_AXIS * np.sin(latitude)) ** 2
-        true_radius = np.sqrt(numerator / denominator)
-
-        latitude = Stroke(self.times, self.geographical[:, 0], "cubic")
-        test_radius = Coordinate._WGS84_radius(self, latitude)(self.times)
-
-        self.assertTrue(np.allclose(true_radius, test_radius, atol=0.01))
-
-    def test_altitude(self):
-
-        test_altitude = self.coor_itrs.altitude()
-        self.assertTrue(np.allclose(test_altitude, self.altitude, atol=0.06))
-
-    def test_distance(self):
-        """Test `Coordinate.distance`."""
-
-        julian = [30462.5, 30462.50069, 30462.50172, 30462.50279, 30462.50386]
-        position = [[6343.81620, -2640.87223, -11.25542802],
-                    [6295.64584, -2718.09271, 443.08233],
-                    [6173.04658, -2808.91831, 1108.58544],
-                    [5980.91111, -2872.96258, 1792.17964],
-                    [5724.30203, -2904.09810, 2460.25799]]
-        velocity = np.zeros((5, 3))
-        location = GroundPosition(52.1579, -106.6702, 0.482)
-        coor = Coordinate(position, velocity, "itrs", julian, self.offset)
-
-        true_distance = [9070.49269, 8776.71795, 8330.99543, 7851.70083, 7359.09190]
-        test_distance = coor.distance(location)
-
-        self.assertTrue(np.allclose(true_distance, test_distance, atol=0.35))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_satellite_old/test_satellite.py b/tests/test_satellite_old/test_satellite.py
deleted file mode 100644
index c820096..0000000
--- a/tests/test_satellite_old/test_satellite.py
+++ /dev/null
@@ -1,88 +0,0 @@
-
-
-from celest.encounter.groundposition import GroundPosition
-from celest.satellite.satellite import Satellite
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-class TestSatellite(TestCase):
-
-    def setUp(self):
-
-        fname = "tests/test_data/coordinate_validation_long.txt"
-        cols = (0, 11, 12, 13, 14, 15, 16)
-        skiprows = 1
-        max_rows = 5000
-        data = np.loadtxt(fname=fname, usecols=cols, skiprows=skiprows,
-                          max_rows=max_rows)
-
-        self.times = data[:, 0]
-        self.GCRS = data[:, 1:4]
-        self.GCRS_vel = data[:, 4:7]
-
-        self.offset = 2430000
-        self.satellite = Satellite(self.GCRS, self.GCRS_vel, "gcrs",
-                                   self.times, self.offset)
-
-        self.julian = self.satellite.julian()
-        self.ut1 = self.satellite.ut1()
-        self.gmst = self.satellite.gmst()
-        self.gast = self.satellite.gast()
-        self.geo = self.satellite.geo()
-        self.gcrs = self.satellite.gcrs()
-        self.itrs = self.satellite.itrs()
-
-    def test_attitude(self):
-
-        location = GroundPosition(52.1579, -106.6702, 0.482)
-        roll, pitch, yaw = self.satellite.attitude(location)
-
-        self.assertIsNotNone(roll)
-        self.assertIsNotNone(pitch)
-        self.assertIsNotNone(yaw)
-
-    def test_save(self):
-
-        times = ("julian", "ut1", "gmst", "gast")
-        positions = ("geo", "gcrs", "itrs")
-        self.satellite.save(times, positions, path="test_data.csv", sep=",")
-
-        data = np.loadtxt("test_data.csv", delimiter=",", skiprows=1)
-        load_julian = data[:, 1]
-        load_ut1 = data[:, 2]
-        load_gmst = data[:, 3]
-        load_gast = data[:, 4]
-        load_lat = data[:, 5]
-        load_lon = data[:, 6]
-        load_alt = data[:, 7]
-        load_gcrs = data[:, 8:14]
-        load_itrs = data[:, 14:20]
-
-        self.assertTrue(np.array_equal(load_julian, self.julian))
-        self.assertTrue(np.array_equal(load_ut1, self.ut1))
-        self.assertTrue(np.array_equal(load_gmst, self.gmst))
-        self.assertTrue(np.array_equal(load_gast, self.gast))
-        self.assertTrue(np.array_equal(load_lat, self.geo[0]))
-        self.assertTrue(np.array_equal(load_lon, self.geo[1]))
-        self.assertTrue(np.array_equal(load_alt, self.geo[2]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 0], self.gcrs[0]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 1], self.gcrs[1]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 2], self.gcrs[2]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 3], self.gcrs[3]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 4], self.gcrs[4]))
-        self.assertTrue(np.array_equal(load_gcrs[:, 5], self.gcrs[5]))
-        self.assertTrue(np.array_equal(load_itrs[:, 0], self.itrs[0]))
-        self.assertTrue(np.array_equal(load_itrs[:, 1], self.itrs[1]))
-        self.assertTrue(np.array_equal(load_itrs[:, 2], self.itrs[2]))
-        self.assertTrue(np.array_equal(load_itrs[:, 3], self.itrs[3]))
-        self.assertTrue(np.array_equal(load_itrs[:, 4], self.itrs[4]))
-        self.assertTrue(np.array_equal(load_itrs[:, 5], self.itrs[5]))
-
-        import os
-        os.remove("test_data.csv")
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_satellite_old/test_time.py b/tests/test_satellite_old/test_time.py
deleted file mode 100644
index a947d0f..0000000
--- a/tests/test_satellite_old/test_time.py
+++ /dev/null
@@ -1,167 +0,0 @@
-
-
-from astropy import coordinates, time
-from celest.satellite.time import Time
-from unittest import TestCase
-import julian as jd
-import unittest
-import numpy as np
-
-
-class TestTime(TestCase):
-
-    def setUp(self):
-
-        self.julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        self.astropy_time = time.Time(self.julian, format="jd")
-
-    def test_true_solar_time(self):
-        """Test `Time.true_solar_time`.
-
-        Notes
-        -----
-        Test cases are generated from the Global Monitoring Labratory and the
-        Government of Australia.[1]_[2]_
-
-        References
-        ----------
-        .. [1] NOAA US Department of Commerce. ESRL Global Monitoring
-           Laboratory -Global Radiation and Aerosols. url:
-           https://gml.noaa.gov/grad/solcalc/.
-        .. [2] Time Conventions. url:
-           http://www.bom.gov.au/climate/data-services/solar/content/data-time.html.
-        """
-
-        julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        longitude = [-105, -118.24, 147.46]
-
-        true_true_solar_time = np.array([23.0715, 0.378059, 10.76556])
-        test_true_solar_time = Time(julian).true_solar_time(longitude)
-
-        self.assertTrue(np.allclose(true_true_solar_time, test_true_solar_time,
-                        atol=0.11))
-
-    def test_mean_solar_time(self):
-        """Test `Time.mean_solar_time`.
-
-        Notes
-        -----
-        Test cases are generated from the Government of Australia.[1]_
-
-        References
-        ----------
-        .. [1] Time Conventions. url:
-           http://www.bom.gov.au/climate/data-services/solar/content/data-time.html.
-        """
-
-        julian = 2456293.5416666665
-        longitude = 147.46
-
-        true_mean_solar_time = np.array([10.83056])
-        test_mean_solar_time = Time(julian).mean_solar_time(longitude)
-
-        self.assertTrue(np.allclose(true_mean_solar_time, test_mean_solar_time,
-                        atol=0.001))
-
-    def test_true_hour_angle(self):
-        """Test `Time.true_hour_angle`.
-
-        Notes
-        -----
-        Test cases are taken from Global Monitoring Labratory.[1]_
-
-        References
-        ----------
-        .. [1] NOAA US Department of Commerce. ESRL Global Monitoring
-           Laboratory -Global Radiation and Aerosols. url:
-           https://gml.noaa.gov/grad/solcalc/.
-        """
-
-        julian = [2455368.75, 2459450.85, 2456293.5416666665]
-        longitude = [-105, -118.24, 147.46]
-
-        true_true_hour_angle = np.array([166.0734, 185.671, 341.53]) / 15 % 24
-        test_true_hour_angle = Time(julian).true_hour_angle(longitude)
-
-        self.assertTrue(np.allclose(true_true_hour_angle, test_true_hour_angle,
-                        atol=1.51))
-
-    def test_mean_hour_angle(self):
-
-        julian = 2456293.5416666665
-        longitude = 147.46
-
-        true_mean_hour_angle = np.array([342.4584]) / 15 % 24
-        test_mean_hour_angle = Time(julian).mean_hour_angle(longitude)
-
-        self.assertTrue(np.allclose(true_mean_hour_angle, test_mean_hour_angle,
-                        atol=0.01))
-
-    def test_ut1(self):
-
-        ut1_utc_diff = self.astropy_time.get_delta_ut1_utc().value / 3600
-
-        true_ut1 = self.astropy_time.to_value("decimalyear") % 1
-        true_ut1 = (true_ut1 * 365 * 24) % 24 + ut1_utc_diff
-        test_ut1 = Time(self.julian).ut1()
-
-        self.assertTrue(np.allclose(true_ut1, test_ut1, atol=0.001))
-
-    def test_julian(self):
-
-        test_julian = Time(julian=self.julian).julian()
-        self.assertTrue(np.array_equal(self.julian, test_julian))
-
-    def test_datetime(self):
-
-        test_datetime = Time(julian=self.julian).datetime()
-
-        for julian, test_datetime_instace in zip(self.julian, test_datetime):
-            true_datetime_instance = jd.from_jd(julian)
-
-            self.assertEqual(test_datetime_instace.year,
-                             true_datetime_instance.year)
-            self.assertEqual(test_datetime_instace.month,
-                             true_datetime_instance.month)
-            self.assertEqual(test_datetime_instace.day,
-                             true_datetime_instance.day)
-            self.assertEqual(test_datetime_instace.second,
-                             true_datetime_instance.second)
-            self.assertAlmostEqual(test_datetime_instace.microsecond,
-                                   true_datetime_instance.microsecond, delta=1)
-
-    def test_gmst(self):
-
-        true_gmst = self.astropy_time.sidereal_time("mean", "greenwich")
-        true_gmst = coordinates.Angle(true_gmst).hour
-        test_gmst = Time(self.julian).gmst()
-
-        self.assertTrue(np.allclose(true_gmst, test_gmst, atol=0.0001))
-
-    def test_lmst(self):
-
-        true_lmst = self.astropy_time.sidereal_time("mean", longitude="150")
-        true_lmst = coordinates.Angle(true_lmst).hour
-        test_lmst = Time(self.julian).lmst(150)
-
-        self.assertTrue(np.allclose(true_lmst, test_lmst, atol=0.1))
-
-    def test_gast(self):
-
-        true_gast = self.astropy_time.sidereal_time("apparent", "greenwich")
-        true_gast = coordinates.Angle(true_gast).hour
-        test_gast = Time(self.julian).gast()
-
-        self.assertTrue(np.allclose(true_gast, test_gast, atol=0.0001))
-
-    def test_last(self):
-
-        true_last = self.astropy_time.sidereal_time("apparent", longitude="150")
-        true_last = coordinates.Angle(true_last).hour
-        test_last = Time(self.julian).last(150)
-        
-        self.assertTrue(np.allclose(true_last, test_last, atol=0.1))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_schedule/request_handler_testing_utils.py b/tests/test_schedule/request_handler_testing_utils.py
index 8d0602f..3d6cfe2 100644
--- a/tests/test_schedule/request_handler_testing_utils.py
+++ b/tests/test_schedule/request_handler_testing_utils.py
@@ -1,20 +1,23 @@
 
 
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._window_handling import VTW
+from celest.coordinates.frames.attitude import Attitude
+from celest.coordinates.ground_location import GroundLocation
+from celest.encounter.window_handling import VisibleTimeWindow
+from celest.units.quantity import Quantity
+from celest import units as u
 import numpy as np
 
 
-toronto = GroundPosition(latitude=43.65, longitude=-79.38, height=0.076)
-north_bay = GroundPosition(latitude=46.31, longitude=-79.46, height=0.193)
-sudbury = GroundPosition(latitude=46.49, longitude=-80.99, height=0.348)
-ottawa = GroundPosition(latitude=45.42, longitude=-75.70, height=0.070)
-kingston = GroundPosition(latitude=44.23, longitude=-76.49, height=0.093)
-niagara_falls = GroundPosition(latitude=43.09, longitude=-79.08, height=0.099)
-london = GroundPosition(latitude=42.98, longitude=-81.24, height=0.251)
-mississauga = GroundPosition(latitude=43.59, longitude=-79.64, height=0.156)
-timmins = GroundPosition(latitude=48.48, longitude=-81.33, height=0.295)
-tobermory = GroundPosition(latitude=45.25, longitude=-81.66, height=0.271)
+toronto = GroundLocation(43.65, -79.38, 0.076, u.deg, u.km)
+north_bay = GroundLocation(46.31, -79.46, 0.193, u.deg, u.km)
+sudbury = GroundLocation(46.49, -80.99, 0.348, u.deg, u.km)
+ottawa = GroundLocation(45.42, -75.70, 0.070, u.deg, u.km)
+kingston = GroundLocation(44.23, -76.49, 0.093, u.deg, u.km)
+niagara_falls = GroundLocation(43.09, -79.08, 0.099, u.deg, u.km)
+london = GroundLocation(42.98, -81.24, 0.251, u.deg, u.km)
+mississauga = GroundLocation(43.59, -79.64, 0.156, u.deg, u.km)
+timmins = GroundLocation(48.48, -81.33, 0.295, u.deg, u.km)
+tobermory = GroundLocation(45.25, -81.66, 0.271, u.deg, u.km)
 
 
 VTW_FILE_PATH = "tests/test_data/test_vtw_data/"
@@ -36,39 +39,38 @@
 
 
 def initialize_request_list():
-
     return [construct_request_from_files(vtw_data_list) for vtw_data_list in VTW_DATA]
 
 
 def construct_request_from_files(vtw_data_list):
-
     name, location, deadline, duration, priority, quality, look_ang = vtw_data_list
-    vtws = get_vtw_list_from_file(name)
+    vtws = get_vtw_list_from_file(name, location)
 
     return [
         location,
-        deadline,
-        duration,
+        Quantity(deadline, u.jd2000),
+        Quantity(duration, u.s),
         priority,
         quality,
-        look_ang,
+        Quantity(look_ang, u.deg) if look_ang is not None else None,
         vtws
     ]
 
 
-def get_vtw_list_from_file(location_string):
-
-    roll, pitch, yaw = get_roll_pitch_yaw_from_file(location_string)
+def get_vtw_list_from_file(location_string, location):
+    attitude = get_roll_pitch_yaw_from_file(location_string, location)
     vtws = load_file(VTW_FILE_PATH + location_string + VTW_FILE_FINALE)
-    return [VTW(vtw[0], vtw[1], roll, pitch, yaw) for vtw in vtws]
+    return [VisibleTimeWindow(vtw[0], vtw[1], attitude) for vtw in vtws]
 
 
-def get_roll_pitch_yaw_from_file(location_string):
-
+def get_roll_pitch_yaw_from_file(location_string, location):
     attitude = load_file(ATTITUDE_FILE_PATH + location_string + ATTITUDE_FILE_FINALE)
-    return attitude[:, 0], attitude[:, 1], attitude[:, 2]
+    julian = np.linspace(0, len(attitude) - 1, len(attitude))
+    roll = attitude[:, 0].reshape((-1,))
+    pitch = attitude[:, 1].reshape((-1,))
+    yaw = attitude[:, 2].reshape((-1,))
+    return Attitude(julian, roll, pitch, yaw, u.deg, location)
 
 
 def load_file(file_name):
-
     return np.loadtxt(file_name, delimiter=",")
diff --git a/tests/test_schedule/test_insertion_operators.py b/tests/test_schedule/test_insertion_operators.py
index e671566..bf1e43a 100644
--- a/tests/test_schedule/test_insertion_operators.py
+++ b/tests/test_schedule/test_insertion_operators.py
@@ -10,7 +10,8 @@
     minimum_opportunity_insertion
 )
 from celest.schedule.request_handler import RequestHandler
-from polare import Stroke
+from celest.units.quantity import Quantity
+from celest import units as u
 from tests.test_schedule.request_handler_testing_utils import initialize_request_list
 from unittest import TestCase
 import numpy as np
@@ -19,7 +20,6 @@
 class TestInsertionOperators(TestCase):
 
     def setUp(self):
-
         self.request_handler = RequestHandler()
 
         request_list = initialize_request_list()
@@ -27,14 +27,12 @@ def setUp(self):
             self.request_handler.add_request(*request)
 
     def test_greedy_insert_one_item(self):
-
         greedy_insertion(self.request_handler, 1)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 1)
         self.assertTrue(self.request_handler[0][0])
         self.assertEqual(self.request_handler[0][7], 7)
 
     def test_greedy_insert_two_items(self):
-
         greedy_insertion(self.request_handler, 2)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 2)
         self.assertTrue(self.request_handler[0][0])
@@ -43,59 +41,58 @@ def test_greedy_insert_two_items(self):
         self.assertEqual(self.request_handler[1][7], 7)
 
     def test_greedy_insert_all_items(self):
-
         greedy_insertion(self.request_handler, 10)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 7)
 
     def test_look_angle_time_with_root(self):
+        look_angle = Quantity(np.linspace(0, 1, 100), u.deg)
+        desired_look_angle = Quantity(0.5, u.deg)
+        julian = Quantity(np.linspace(0, 1, 100), u.jd2000)
+        start = Quantity(0, u.jd2000)
+        end = Quantity(1, u.jd2000)
+        look_angle_time = _look_angle_time(look_angle, desired_look_angle, julian, start, end)
 
-        time = np.linspace(0, 1, 100)
-        look_angle_stroke = Stroke(time, time)
-        look_angle_time = _look_angle_time(look_angle_stroke, 0.5, 0, 1, )
-
-        self.assertAlmostEqual(look_angle_time, 0.5, delta=1e-6)
+        self.assertAlmostEqual(look_angle_time.data, 0.5, delta=1e-2)
 
     def test_look_angle_time_with_no_root(self):
+        look_angle = Quantity(np.linspace(0, 1, 100), u.deg)
+        desired_look_angle = Quantity(2, u.deg)
+        julian = Quantity(np.linspace(0, 1, 100), u.jd2000)
+        start = Quantity(0, u.jd2000)
+        end = Quantity(1, u.jd2000)
+        look_angle_time = _look_angle_time(look_angle, desired_look_angle,
+                                           julian, start, end)
 
-        time = np.linspace(0, 1, 100)
-        look_angle_stroke = Stroke(time, time)
-
-        self.assertIsNone(_look_angle_time(look_angle_stroke, 2, 0, 1))
+        self.assertIsNone(look_angle_time)
 
     def test_image_quality_with_nadir_time(self):
-
         vtw = self.request_handler[0][10][0]
         nadir_time = (vtw.rise_time + vtw.set_time) / 2
         self.assertEqual(image_quality(vtw, nadir_time), 10)
 
     def test_image_quality_with_start_time(self):
-
         vtw = self.request_handler[0][10][0]
         start_time = vtw.rise_time
         self.assertEqual(image_quality(vtw, start_time), 1)
 
     def test_image_quality_with_set_time(self):
-
         vtw = self.request_handler[0][10][0]
         start_time = vtw.set_time
         self.assertEqual(image_quality(vtw, start_time), 1)
 
     def test_image_quality_is_met_with_nadir_time(self):
-
         vtw = self.request_handler[0][10][0]
         nadir_time = (vtw.rise_time + vtw.set_time) / 2
         self.assertTrue(image_quality_is_met(vtw, nadir_time, 9))
         self.assertTrue(image_quality_is_met(vtw, nadir_time, 10))
 
     def test_image_quality_is_met_with_start_time(self):
-
         vtw = self.request_handler[0][10][0]
         start_time = vtw.rise_time
         self.assertTrue(image_quality_is_met(vtw, start_time, 1))
         self.assertFalse(image_quality_is_met(vtw, start_time, 2))
 
     def test_conflict_exists_with_no_conflicts(self):
-
         self._schedule_first_request()
 
         rise_time = self.request_handler[0][10][0].rise_time
@@ -107,14 +104,13 @@ def test_conflict_exists_with_no_conflicts(self):
         self.assertFalse(does_conflict_exists(self.request_handler, start, duration_in_days))
 
     def _schedule_first_request(self):
-
         rise_time = self.request_handler[0][10][0].rise_time
         set_time = self.request_handler[0][10][0].set_time
+        duration = Quantity(86400 * (set_time.to(u.jd2000).data - rise_time.to(u.jd2000).data), u.s)
 
-        self.request_handler.schedule_request(0, 0, rise_time, 86400 * (set_time - rise_time))
+        self.request_handler.schedule_request(0, 0, rise_time, duration)
 
     def test_conflict_exists_with_window_around_scheduled_window(self):
-
         self._schedule_first_request()
 
         rise_time = self.request_handler[0][10][0].rise_time
@@ -126,7 +122,6 @@ def test_conflict_exists_with_window_around_scheduled_window(self):
         self.assertTrue(does_conflict_exists(self.request_handler, start, duration_in_days))
 
     def test_conflict_exists_with_start_time_in_scheduled_window(self):
-
         self._schedule_first_request()
 
         rise_time = self.request_handler[0][10][0].rise_time
@@ -138,7 +133,6 @@ def test_conflict_exists_with_start_time_in_scheduled_window(self):
         self.assertTrue(does_conflict_exists(self.request_handler, start, duration_in_days))
 
     def test_conflict_exists_with_end_time_in_scheduled_window(self):
-
         self._schedule_first_request()
 
         rise_time = self.request_handler[0][10][0].rise_time
@@ -150,7 +144,6 @@ def test_conflict_exists_with_end_time_in_scheduled_window(self):
         self.assertTrue(does_conflict_exists(self.request_handler, start, duration_in_days))
 
     def test_conflict_exists_with_window_in_scheduled_window(self):
-
         self._schedule_first_request()
 
         rise_time = self.request_handler[0][10][0].rise_time
@@ -162,14 +155,12 @@ def test_conflict_exists_with_window_in_scheduled_window(self):
         self.assertTrue(does_conflict_exists(self.request_handler, start, duration_in_days))
 
     def test_minimum_opportunity_insert_one_item(self):
-
         minimum_opportunity_insertion(self.request_handler, 1)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 1)
         self.assertTrue(self.request_handler[0][0])
         self.assertEqual(len(self.request_handler[0][10]), 20)
 
     def test_minimum_opportunity_insert_two_items(self):
-
         minimum_opportunity_insertion(self.request_handler, 2)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 2)
         self.assertTrue(self.request_handler[0][0])
@@ -178,24 +169,20 @@ def test_minimum_opportunity_insert_two_items(self):
         self.assertEqual(len(self.request_handler[1][10]), 19)
 
     def test_minimum_opportunity_insert_all_items(self):
-
         minimum_opportunity_insertion(self.request_handler, 10)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 8)
 
     def test_minimum_conflict_insert_one_item(self):
-
         minimum_conflict_insertion(self.request_handler, 1)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 1)
         self.assertTrue(self.request_handler[0][0])
 
     def test_minimum_conflict_insert_two_items(self):
-
         minimum_conflict_insertion(self.request_handler, 2)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 2)
         self.assertTrue(self.request_handler[0][0])
         self.assertTrue(self.request_handler[1][0])
 
     def test_minimum_conflict_insert_all_items(self):
-
         minimum_conflict_insertion(self.request_handler, 10)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 4)
diff --git a/tests/test_schedule/test_request_handler.py b/tests/test_schedule/test_request_handler.py
index fcbbb06..e90c551 100644
--- a/tests/test_schedule/test_request_handler.py
+++ b/tests/test_schedule/test_request_handler.py
@@ -8,22 +8,19 @@
 class TestRequestHandler(TestCase):
 
     def setUp(self) -> None:
-
         self.request_handler = RequestHandler()
         for request in initialize_request_list():
             self.request_handler.add_request(*request)
 
     def test_getitem(self):
-
-        self.assertListEqual(self.request_handler[0], self.request_handler.requests[0])
+        self.assertListEqual(self.request_handler[0],
+                             self.request_handler.requests[0])
 
     def test_request_handler_iteration(self):
-
         for request in self.request_handler:
             self.assertIsInstance(request, list)
 
     def test_add_request(self):
-
         request_list = initialize_request_list()
 
         self.assertListEqual(self.request_handler[0][4:-1], request_list[0][:-1])
@@ -40,7 +37,6 @@ def test_add_request(self):
         self.assertEqual(self.request_handler.number_of_requests, len(request_list))
 
     def test_schedule_request(self):
-
         self.request_handler.schedule_request(0, 1, 2, 3)
 
         self.assertTrue(self.request_handler[0][RequestIndices.is_scheduled])
@@ -49,7 +45,6 @@ def test_schedule_request(self):
         self.assertEqual(self.request_handler[0][RequestIndices.scheduled_duration], 3)
 
     def test_unschedule_all_requests(self):
-
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.request_handler.schedule_request(1, 1, 2, 3)
         self.request_handler.unschedule_all_requests()
@@ -66,7 +61,6 @@ def test_unschedule_all_requests(self):
         self.assertFalse(self.request_handler[9][RequestIndices.is_scheduled])
 
     def test_unschedule_request(self):
-
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.request_handler.unschedule_request(0)
 
@@ -76,79 +70,66 @@ def test_unschedule_request(self):
         self.assertIsNone(self.request_handler[0][RequestIndices.scheduled_duration])
 
     def test_is_request_scheduled(self):
-
         self.assertFalse(self.request_handler.is_request_scheduled(0))
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.assertTrue(self.request_handler.is_request_scheduled(0))
 
     def test_is_request_scheduled(self):
-
         self.assertFalse(self.request_handler.is_request_scheduled(0))
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.assertTrue(self.request_handler.is_request_scheduled(0))
 
     def test_vtw_index(self):
-
         self.assertIsNone(self.request_handler.vtw_index(0))
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.assertEqual(self.request_handler.vtw_index(0), 1)
 
     def test_start_time(self):
-
         self.assertIsNone(self.request_handler.start_time(0))
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.assertEqual(self.request_handler.start_time(0), 2)
 
     def test_scheduled_duration(self):
-
         self.assertIsNone(self.request_handler.scheduled_duration(0))
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.assertEqual(self.request_handler.scheduled_duration(0), 3)
 
     def test_deadline(self):
-
         test_deadline = self.request_handler.deadline(0)
         actual_deadline = self.request_handler[0][RequestIndices.deadline]
         self.assertEqual(test_deadline, actual_deadline)
 
     def test_duration(self):
-
         test_duration = self.request_handler.duration(0)
         actual_duration = self.request_handler[0][RequestIndices.duration]
         self.assertEqual(test_duration, actual_duration)
 
     def test_priority(self):
-
         test_priority = self.request_handler.priority(0)
         actual_priority = self.request_handler[0][RequestIndices.priority]
         self.assertEqual(test_priority, actual_priority)
 
     def test_image_quality(self):
-
         test_image_quality = self.request_handler.image_quality(0)
         actual_image_quality = self.request_handler[0][RequestIndices.quality]
         self.assertEqual(test_image_quality, actual_image_quality)
 
     def test_look_angle(self):
-
         test_look_angle = self.request_handler.look_angle(0)
         actual_look_angle = self.request_handler[0][RequestIndices.look_angle]
         self.assertEqual(test_look_angle, actual_look_angle)
 
     def test_vtw_list(self):
-
         test_vtw_list = self.request_handler.vtw_list(0)
         actual_vtw_list = self.request_handler[0][RequestIndices.vtw_list]
         self.assertEqual(test_vtw_list, actual_vtw_list)
 
     def test_number_of_scheduled_items(self):
-
         self.request_handler.schedule_request(0, 1, 2, 3)
         self.request_handler.schedule_request(1, 1, 2, 3)
         self.assertEqual(self.request_handler.number_of_scheduled_requests, 2)
 
     def test_sort_by_decreasing_priority(self):
-
         self.request_handler.sort_by_decreasing_priority()
 
         self.assertEqual(self.request_handler[0][RequestIndices.priority], 7)
@@ -163,7 +144,6 @@ def test_sort_by_decreasing_priority(self):
         self.assertEqual(self.request_handler[9][RequestIndices.priority], 1)
 
     def test_sort_by_decreasing_opportunity(self):
-
         self.request_handler.sort_by_decreasing_opportunity()
 
         self.assertEqual(len(self.request_handler[0][RequestIndices.vtw_list]), 20)
@@ -178,20 +158,20 @@ def test_sort_by_decreasing_opportunity(self):
         self.assertEqual(len(self.request_handler[9][RequestIndices.vtw_list]), 18)
 
     def test_sort_vtws_by_increasing_conflict_degree(self):
-
         self.request_handler.sort_vtws_by_increasing_conflict_degree()
         for request in self.request_handler:
-            previous_conflict_degree = 0
+            previous_conflict_degree = None
             for vtw in request[RequestIndices.vtw_list]:
                 current_conflict_degree = self.request_handler._conflict_degree(vtw)
-                self.assertLessEqual(previous_conflict_degree, current_conflict_degree)
+                if previous_conflict_degree is not None:
+                    self.assertLessEqual(previous_conflict_degree, current_conflict_degree)
                 previous_conflict_degree = current_conflict_degree
 
     def test_sort_vtws_by_increasing_rise_time(self):
-
         self.request_handler.sort_vtws_by_increasing_rise_time()
         for request in self.request_handler:
-            previous_rise_time = 0
+            previous_rise_time = None
             for vtw in request[RequestIndices.vtw_list]:
-                self.assertLessEqual(previous_rise_time, vtw.rise_time)
+                if previous_rise_time is not None:
+                    self.assertLessEqual(previous_rise_time, vtw.rise_time)
                 previous_rise_time = vtw.rise_time
diff --git a/tests/test_schedule/test_schedule.py b/tests/test_schedule/test_schedule.py
index 07d831f..60a60bf 100644
--- a/tests/test_schedule/test_schedule.py
+++ b/tests/test_schedule/test_schedule.py
@@ -1,13 +1,20 @@
 
 
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._window_handling import OWHandler
-from celest.satellite.satellite import Satellite
+from celest.coordinates.frames.gcrs import GCRS
+from celest.coordinates.ground_location import GroundLocation
+from celest.encounter.window_handling import WindowHandler
+from celest.satellite import Satellite
 from celest.schedule.insertion_operators import _INSERTION_FUNCTIONS
-from celest.schedule.request_handler import RequestIndices
 from celest.schedule.removal_operators import _REMOVAL_FUNCTIONS
+from celest.schedule.request_handler import RequestIndices
 from celest.schedule.schedule import Schedule
-from celest.schedule.scheduling_utils import cost, is_complete, initialize_solution
+from celest.schedule.scheduling_utils import (
+    cost,
+    initialize_solution,
+    is_complete
+)
+from celest.units.quantity import Quantity
+from celest import units as u
 from unittest import TestCase
 import numpy as np
 
@@ -18,10 +25,13 @@ def setUp(self):
 
         fname = 'Tests/test_data/coordinate_validation_long.txt'
         data = np.loadtxt(fname, usecols=[0, 11, 12, 13, 14, 15, 16], skiprows=1, max_rows=5000)
-        julian, gcrs_pos, gcrs_vel = data[:, 0], data[:, 1:4], data[:, 4:]
+        julian, gcrs_pos, gcrs_vel = data[:, 0] + 2430000, data[:, 1:4], data[:, 4:]
+
+        gcrs_position = GCRS(julian, gcrs_pos[:, 0], gcrs_pos[:, 1], gcrs_pos[:, 2], u.km)
+        gcrs_velocity = GCRS(julian, gcrs_vel[:, 0], gcrs_vel[:, 1], gcrs_vel[:, 2], u.m / u.s)
 
-        self.satellite = Satellite(gcrs_pos, gcrs_vel, "gcrs", julian, 2430000)
-        self.toronto = GroundPosition(43.6532, -79.3832, 0.076)
+        self.satellite = Satellite(gcrs_position, gcrs_velocity)
+        self.toronto = GroundLocation(43.6532, -79.3832, 0.076, u.deg, u.km)
 
     def test_schedule_raises_exception_if_null_satellite_is_passed(self):
 
@@ -54,37 +64,40 @@ def test_generate_raises_exception_when_no_requests_are_added(self):
         schedule = Schedule(self.satellite, 10)
         self.assertRaises(Exception, schedule.generate, 100, 0.8, 0.5)
 
-    def test_generate_returns_OWHandler_when_requests_are_added(self):
+    def test_generate_returns_WindowHandler_when_requests_are_added(self):
 
         schedule = Schedule(self.satellite, 10)
-        schedule.add_request(self.toronto, 2460467, 30, 1, 1, None)
-        self.assertIsInstance(schedule.generate(100, 0.8, 0.5), OWHandler)
+        schedule.add_request(self.toronto, Quantity(2460467, u.jd2000), Quantity(30, u.s), 1, 1, None)
+        self.assertIsInstance(schedule.generate(100, 0.8, 0.5), WindowHandler)
 
     def test_generate_returns_improved_solutions(self):
 
-        toronto = GroundPosition(43.65, -79.38, 0.076)
-        north_bay = GroundPosition(46.31, -79.46, 0.193)
-        sudbury = GroundPosition(46.49, -80.99, 0.348)
-        ottawa = GroundPosition(45.42, -75.70, 0.070)
-        kingston = GroundPosition(44.23, -76.49, 0.093)
-        niagara_falls = GroundPosition(43.09, -79.08, 0.099)
-        london = GroundPosition(42.98, -81.24, 0.251)
-        mississauga = GroundPosition(43.59, -79.64, 0.156)
-        timmins = GroundPosition(48.48, -81.33, 0.295)
-        tobermory = GroundPosition(45.25, -81.66, 0.271)
+        toronto = GroundLocation(43.65, -79.38, 0.076, u.deg, u.km)
+        north_bay = GroundLocation(46.31, -79.46, 0.193, u.deg, u.km)
+        sudbury = GroundLocation(46.49, -80.99, 0.348, u.deg, u.km)
+        ottawa = GroundLocation(45.42, -75.70, 0.070, u.deg, u.km)
+        kingston = GroundLocation(44.23, -76.49, 0.093, u.deg, u.km)
+        niagara_falls = GroundLocation(43.09, -79.08, 0.099, u.deg, u.km)
+        london = GroundLocation(42.98, -81.24, 0.251, u.deg, u.km)
+        mississauga = GroundLocation(43.59, -79.64, 0.156, u.deg, u.km)
+        timmins = GroundLocation(48.48, -81.33, 0.295, u.deg, u.km)
+        tobermory = GroundLocation(45.25, -81.66, 0.271, u.deg, u.km)
 
         schedule = Schedule(self.satellite, 10)
 
-        schedule.add_request(toronto, 2460467, 30, 1, 1, None)
-        schedule.add_request(north_bay, 2460467, 30, 1, 1, None)
-        schedule.add_request(sudbury, 2460467, 30, 4, 1, None)
-        schedule.add_request(ottawa, 2460467, 30, 2, 1, None)
-        schedule.add_request(kingston, 2460467, 30, 7, 1, None)
-        schedule.add_request(niagara_falls, 2460467, 30, 3, 1, None)
-        schedule.add_request(london, 2460467, 30, 4, 1, None)
-        schedule.add_request(mississauga, 2460467, 30, 5, 1, None)
-        schedule.add_request(timmins, 2460467, 30, 1, 1, None)
-        schedule.add_request(tobermory, 2460467, 30, 7, 1, None)
+        deadline = Quantity(2460467, u.jd2000)
+        duration = Quantity(30, u.s)
+
+        schedule.add_request(toronto, deadline, duration, 1, 1, None)
+        schedule.add_request(north_bay, deadline, duration, 1, 1, None)
+        schedule.add_request(sudbury, deadline, duration, 4, 1, None)
+        schedule.add_request(ottawa, deadline, duration, 2, 1, None)
+        schedule.add_request(kingston, deadline, duration, 7, 1, None)
+        schedule.add_request(niagara_falls, deadline, duration, 3, 1, None)
+        schedule.add_request(london, deadline, duration, 4, 1, None)
+        schedule.add_request(mississauga, deadline, duration, 5, 1, None)
+        schedule.add_request(timmins, deadline, duration, 1, 1, None)
+        schedule.add_request(tobermory, deadline, duration, 7, 1, None)
 
         initialize_solution(schedule.request_handler)
         initial_solution_cost = cost(schedule.request_handler)
diff --git a/tests/test_time.py b/tests/test_time.py
index 287a245..1aaed6a 100644
--- a/tests/test_time.py
+++ b/tests/test_time.py
@@ -24,7 +24,7 @@ def test_true_solar_time_for_validation(self):
 
         Notes
         -----
-        Test cases are generated from the Global Monitoring Labratory and the
+        Test cases are generated from the Global Monitoring Laboratory and the
         Government of Australia.[1]_[2]_
 
         References
@@ -72,7 +72,7 @@ def test_true_hour_angle_for_validation(self):
 
         Notes
         -----
-        Test cases are taken from Global Monitoring Labratory.[1]_
+        Test cases are taken from Global Monitoring Laboratory.[1]_
 
         References
         ----------
diff --git a/tests/test_units/test_quantity.py b/tests/test_units/test_quantity.py
index afdd889..290ab0c 100644
--- a/tests/test_units/test_quantity.py
+++ b/tests/test_units/test_quantity.py
@@ -156,6 +156,72 @@ def test_rtruediv_quantity_of_different_dimensions(self):
         self.assertEqual(repr(self.compound_unit / self.simple_unit),
                          repr(result.unit))
 
+    def test_eq_when_equal_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity == self.simple_quantity)
+
+    def test_eq_when_not_equal_with_same_dimensions(self):
+        self.assertFalse(self.simple_quantity == (self.simple_quantity + 1))
+
+    def test_eq_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__eq__,
+                          self.compound_quantity)
+
+    def test_neq_when_not_equal_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity != (self.simple_quantity + 1))
+
+    def test_neq_when_equal_with_same_dimensions(self):
+        self.assertFalse(self.simple_quantity != self.simple_quantity)
+
+    def test_neq_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__ne__,
+                          self.compound_quantity)
+
+    def test_lt_when_less_than_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity < (self.simple_quantity + 1))
+
+    def test_lt_when_greater_than_with_same_dimensions(self):
+        self.assertFalse((self.simple_quantity + 1) < self.simple_quantity)
+
+    def test_lt_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__lt__,
+                          self.compound_quantity)
+
+    def test_le_when_less_than_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity <= (self.simple_quantity + 1))
+
+    def test_le_when_equal_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity <= self.simple_quantity)
+
+    def test_le_when_greater_than_with_same_dimensions(self):
+        self.assertFalse((self.simple_quantity + 1) <= self.simple_quantity)
+
+    def test_le_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__le__,
+                          self.compound_quantity)
+
+    def test_gt_when_greater_with_same_dimensions(self):
+        self.assertTrue((self.simple_quantity + 1) > self.simple_quantity)
+
+    def test_gt_when_less_with_same_dimensions(self):
+        self.assertFalse(self.simple_quantity > (self.simple_quantity + 1))
+
+    def test_gt_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__gt__,
+                          self.compound_quantity)
+
+    def test_ge_when_greater_than_with_same_dimensions(self):
+        self.assertTrue((self.simple_quantity + 1) >= self.simple_quantity)
+
+    def test_ge_when_equal_with_same_dimensions(self):
+        self.assertTrue(self.simple_quantity >= self.simple_quantity)
+
+    def test_ge_when_less_than_with_same_dimensions(self):
+        self.assertFalse(self.simple_quantity >= (self.simple_quantity + 1))
+
+    def test_ge_with_different_dimensions_raises_arithmetic_error(self):
+        self.assertRaises(ArithmeticError, self.simple_quantity.__ge__,
+                          self.compound_quantity)
+
     def test_to_with_same_dimension(self):
         simple_quantity_in_km = self.simple_quantity.to(u.km)
         self.assertEqual(0.005, simple_quantity_in_km.data)