Repository: JaiWillems/celest (language: python)

Issue:
MAINT: Add Polare integration on 0.3.x
The integration of the Polare library changes the backend operation of Celest. Primarily, it allows conversions to occur in an on-time-compilation fashion that takes advantage of data pre-interpolation. Having conversions return Polare's `Stroke` objects on command were used in window generation to get window encounter times to arbitrary precision. The change in the base framework will also provide significant benefits when determining satellite positions during window encounters for work in target tracking.


Discussion / hints:

Fix patch:
diff --git a/CITATION.cff b/CITATION.cff
new file mode 100644
index 0000000..f81291d
--- /dev/null
+++ b/CITATION.cff
@@ -0,0 +1,11 @@
+cff-version: 1.2.0
+message: Citation for the Celest library.
+title: Celest
+abstract: Satellite dynamics and mission planning software.
+authors:
+  - family-names: Willems
+    given-names: Jai
+version: 0.2.0
+date-released: "2021-11-13"
+license: BSD-3-Clause
+repository-code: "https://github.com/JaiWillems/celest"
\ No newline at end of file
diff --git a/README.md b/README.md
index d6d85e8..0fa0ad2 100644
--- a/README.md
+++ b/README.md
@@ -1,17 +1,17 @@
-# Celest 0.2.0
+# <a href="https://celest.readthedocs.io/en/latest/"><img alt="Celest" src="/branding/logo/transparent_background_large.png" height="150"></a>
 
-[![PyPI version](https://badge.fury.io/py/celest.svg)](https://badge.fury.io/py/celest) [![license](https://img.shields.io/pypi/l/celest)](https://img.shields.io/pypi/l/celest) [![downloads](https://img.shields.io/pypi/dm/celest)](https://img.shields.io/pypi/dm/celest) [![Documentation Status](https://readthedocs.org/projects/celest/badge/?version=latest)](https://celest.readthedocs.io/en/latest/?badge=latest) [![python](https://img.shields.io/pypi/pyversions/celest)](https://img.shields.io/pypi/pyversions/celest) [![format](https://img.shields.io/pypi/wheel/celest)](https://img.shields.io/pypi/wheel/celest)
+[![PyPI version](https://badge.fury.io/py/celest.svg)](https://badge.fury.io/py/celest) [![license](https://img.shields.io/pypi/l/celest)](https://img.shields.io/pypi/l/celest) [![Downloads](https://static.pepy.tech/personalized-badge/celest?period=total&units=international_system&left_color=grey&right_color=brightgreen&left_text=downloads)](https://pepy.tech/project/celest) [![Documentation Status](https://readthedocs.org/projects/celest/badge/?version=latest)](https://celest.readthedocs.io/en/latest/?badge=latest) [![python](https://img.shields.io/pypi/pyversions/celest)](https://img.shields.io/pypi/pyversions/celest) [![format](https://img.shields.io/pypi/wheel/celest)](https://img.shields.io/pypi/wheel/celest)
 
-Celest is a satellite dynamics and mission planning library designed for the University of Toronto Aerospace Team which is applying the program to a hyperspectral imaging CubeSat mission. The library aims to provide the necessary tools to plan satellite-ground interactions such as ground target tracking, and calculating imaging and transmission opportunities.
+Celest is a satellite dynamics and mission planning library designed for the University of Toronto Aerospace Team which is applying the program to a hyperspectral imaging CubeSat mission. The library aims to provide tools to plan satellite-ground interactions such as ground target tracking, and calculating imaging and transmission opportunities.
 
 * **Documentation:** https://celest.readthedocs.io/en/latest/
-* **Source Code:** https://github.com/JaiWillems/Celest
+* **Source Code:** https://github.com/JaiWillems/celest
 * **PyPI:** https://pypi.org/project/celest/
 * **Bug Report or Feature Request:** https://github.com/JaiWillems/Celest/issues
 
 ## Installation
 
-Celest can be installed from PyPI using the following command:
+Celest can be installed from PyPI using the command line:
 
 ```terminal
 pip install celest
@@ -19,12 +19,12 @@ pip install celest
 
 ## What can the program do?
 
-The program only requires an input data set for a satellite that contains the positions in either a gcrs, itrs, or geographical frame and the assosciated Julian times. The program can then be used to calculate various position and time representations, satellite to ground encounters, and transformations necessary for ground location tracking.
+The program requires an input data set containing satellite positions in the gcrs, itrs, or geographical frame with corresponding Julian times. The program can then calculate various position and time representations, satellite to ground encounters, and transformations necessary for ground location tracking.
 
 ## What can be expected of the project?
 
 In future releases, a high order c-implemented orbit propagator will be incorporated into the library to allow for a complete mission planning workflow from inputting orbital parameters to calculating satellite-to-ground encounter opportunities.
 
-We also intend to incorporate a window scheduling system that will be used to develop an optimal satellite itinerary for satellite mission planning.
+We also intend to incorporate a window scheduling algorithm to create an optimal satellite itinerary for satellite mission planning.
 
-Other features that are being considered include a terminal interface, higher fidelity conversion calculations, window/data visuals, data loaders, and more robust ground tracking by incorporating forward motion compensation.
+Other features under consideration include a terminal interface, higher fidelity conversions, window/data visuals, data loaders, and more robust ground tracking by incorporating forward motion compensation.
diff --git a/branding/logo/transparent_background_large.png b/branding/logo/transparent_background_large.png
new file mode 100644
index 0000000..ed819b6
Binary files /dev/null and b/branding/logo/transparent_background_large.png differ
diff --git a/branding/logo/transparent_background_small.png b/branding/logo/transparent_background_small.png
new file mode 100644
index 0000000..dd4ddac
Binary files /dev/null and b/branding/logo/transparent_background_small.png differ
diff --git a/branding/logo/white_background.png b/branding/logo/white_background.png
new file mode 100644
index 0000000..416b2f8
Binary files /dev/null and b/branding/logo/white_background.png differ
diff --git a/celest/core/__init__.py b/celest/core/__init__.py
deleted file mode 100644
index 2f2527b..0000000
--- a/celest/core/__init__.py
+++ /dev/null
@@ -1,3 +0,0 @@
-
-from celest.core.decorators import set_module
-from celest.core.interpolation import _interpolate
diff --git a/celest/core/decorators.py b/celest/core/decorators.py
deleted file mode 100644
index f37d8c8..0000000
--- a/celest/core/decorators.py
+++ /dev/null
@@ -1,14 +0,0 @@
-"""Define decorators used throughout Celest."""
-
-
-def set_module(module):
-    """Override the module of a class or function."""
-
-    def decorator(func):
-
-        if module is not None:
-            func.__module__ = module
-
-        return func
-
-    return decorator
diff --git a/celest/core/interpolation.py b/celest/core/interpolation.py
deleted file mode 100644
index 84ffc22..0000000
--- a/celest/core/interpolation.py
+++ /dev/null
@@ -1,104 +0,0 @@
-"""General interpolation functionality.
-
-The module contains the functions to allow for general interpolation of column
-data without an independent data array.
-"""
-
-
-from celest.core.decorators import set_module
-from scipy.interpolate import interp1d
-import numpy as np
-
-
-def _interpolate(data: np.ndarray, factor: int=1, dt: int=0, indices:
-                np.ndarray=None) -> np.ndarray:
-    """Interpolate `data` per the given parameters.
-
-    This method takes a multicolumn data array and interpolates it along
-    each column, in the same manner, using the interpolation parameters.
-
-    Parameters
-    ----------
-    data : np.ndarray
-        Array of shape (i, j) containing columns of similar information. Note
-        that the data must be two dimensional.
-    factor : int, optional
-        The factor increase in the number of sample points of interpolated
-        regions.
-    dt : int, optional
-        A positive integer defining the number of points laying beyond the
-        `indices` range to interpolate within.
-    indices : np.ndarray, optional
-        Array of non-empty arrays containing the interpolation indices in a
-        monotonically increasing fashion.
-
-    Returns
-    -------
-    np.ndarray
-        The `data` parameter interpolated per the given parameters.
-
-    Notes
-    -----
-    The `_interp` method was designed to allow for greater precision data
-    in regions where it is important. As a result, it allows for the
-    interpolation of full data sets or only where may be necessary.
-
-    The `dt` parameter is only relavent if an `indices` parameter is given.
-
-    Examples
-    --------
-    >>> data_1 = np.linspace(0, 9, 10).reshape((-1, 1))
-    >>> data_2 = np.linspace(5, 14, 10).reshape((-1, 1)))
-    >>> data = np.concatenate((data_1, data_2, axis=1)
-    >>> ind = np.array([[2, 4]]
-    >>> Interpolation()._interp(data=data, factor=2, dt=1, indices=ind))
-    np.array([[ 0.          5.        ]
-              [ 1.          6.        ]
-              [ 1.57142857  6.57142857]
-              [ 2.14285714  7.14285714]
-              [ 2.71428571  7.71428571]
-              [ 3.28571429  8.28571429]
-              [ 3.85714286  8.85714286]
-              [ 4.42857143  9.42857143]
-              [ 5.         10.        ]
-              [ 6.         11.        ]
-              [ 7.         12.        ]
-              [ 8.         13.        ]
-              [ 9.         14.        ]])
-    """
-
-    j = data.shape[0]
-
-    indep = np.linspace(0, j - 1, j, dtype=int)
-    interp = interp1d(x=indep, y=data, axis=0, kind="cubic")
-
-    if (indices is None) or (indices.size == 0):
-
-        indep_new = np.linspace(0, j - 1, j * factor)
-
-    else:
-
-        indep_new = np.copy(indep)
-
-        for reg in np.flipud(indices):
-
-            min_I, max_I = int(reg[0] - dt), int(reg[-1] + dt)
-
-            if min_I < 0:
-                min_I = 0
-
-            if max_I > indep[-1]:
-                max_I = indep[-1]
-
-            alpha = 1 if factor == 1 else 0
-
-            left_reg = indep_new[:min_I]
-
-            new_reg = np.linspace(min_I, max_I, factor * (max_I - min_I + alpha))
-            right_reg = indep_new[max_I + 1:]
-
-            indep_new = np.concatenate((left_reg, new_reg, right_reg))
-
-    data_new = interp(indep_new)
-
-    return data_new
diff --git a/celest/encounter/__init__.py b/celest/encounter/__init__.py
index 86d34a4..e31b725 100644
--- a/celest/encounter/__init__.py
+++ b/celest/encounter/__init__.py
@@ -1,4 +1,3 @@
 
 from celest.encounter.groundposition import GroundPosition
-from celest.encounter.attitudes import attitudes
 from celest.encounter import windows
diff --git a/celest/encounter/_encounter_math_utils.py b/celest/encounter/_encounter_math_utils.py
deleted file mode 100644
index 08efaa3..0000000
--- a/celest/encounter/_encounter_math_utils.py
+++ /dev/null
@@ -1,193 +0,0 @@
-"""Encounter indices utilities.
-
-This module contains the geometric utilities necessary for analytical pass
-analysis implementation.
-"""
-
-
-from typing import Literal
-import numpy as np
-
-
-def _cone_constraint(theta: np.array, U: np.array, X: np.array) -> np.array:
-    """Return an array of indices for `X` points that fall within a cone.
-
-    This method defines a double-sided cone with an apex located at `U` and an
-    aperture angle `theta` and returns the indices for `X` points that fall
-    within this geometry.
-
-    Parameters
-    ----------
-    theta : np.array
-        Array of size (n,) containing cone aperture angles in degrees.
-    U : np.array
-        Array of shape (n, 3) containing the cartesian coordinates of the apex.
-    X : np.array
-        Array of shape (n, 3) containing the cartesian coordinate points to
-        check.
-
-    Returns
-    -------
-    np.array
-        Array of indices corresponding to points that fall within the cone.
-    """
-
-    U_norm = np.repeat(np.linalg.norm(U, axis=1).reshape((-1, 1)), 3, axis=1)
-    D = np.divide(U, U_norm)
-    Y = X - U
-
-    d1, d2, d3 = D[:, 0], D[:, 1], D[:, 2]
-    y1, y2, y3 = Y[:, 0], Y[:, 1], Y[:, 2]
-
-    ct2 = np.cos(np.radians(theta)) ** 2
-
-    r1 = (d1 ** 2 - ct2) * y1 + d1 * d2 * y2 + d1 * d3 * y3
-    r2 = d2 * d1 * y1 + (d2 ** 2 - ct2) * y2 + d2 * d3 * y3
-    r3 = d3 * d1 * y1 + d3 * d2 * y2 + (d3 ** 2 - ct2) * y3
-
-    point = y1 * r1 + y2 * r2 + y3 * r3
-
-    ind = np.where(point >= 0)[0]
-
-    return ind
-
-
-def _plane_constraint(U: np.array, X: np.array) -> np.array:
-    """Return an array of indices containing points that fall above a plane.
-
-    Parameters
-    ----------
-    U, X : np.array
-        Arrays of shape (n, 3) containing cartesian vectors defining a plane.
-
-    Returns
-    -------
-    np.array
-        Array of points that fall above the plane.
-
-    Notes
-    -----
-    A point :math:`\textbf{X}=\langle x, y, z\rangle^T` falls falls above the
-    plane defined with :math:`\textbf{U}=\langle x_0, y_0, z_0\rangle^T` if it
-    satisfies the following inequality:
-
-    .. math:: 0\leq\left(\textbf{X}-\textbf{U}\right)^T\cdot\textbf{U}
-    """
-
-    Y = X - U
-
-    u1, u2, u3 = U[:, 0], U[:, 1], U[:, 2]
-    y1, y2, y3 = Y[:, 0], Y[:, 1], Y[:, 2]
-
-    point = y1 * u1 + y2 * u2 + y3 * u3
-
-    ind = np.where(point >= 0)[0]
-
-    return ind
-
-
-def _aperature_theta(ang: float, form: Literal[0, 1], n: int=None, U:
-                    np.array=None, X: np.array=None) -> np.array:
-    """Calculate cone aperature angles from constraint angles.
-
-    This method calculates the aperture angles for the cone that defines a
-    quasi-valid encounter region based on the input angle type.
-
-    Parameters
-    ----------
-    ang : float
-        The constraint angle in degrees.
-    form : {0, 1}
-        Defines the angle as altitude if `angType=0` and off-nadir if
-        `angType=1`.
-    n : int, optional
-        Length of returned theta array.
-    U : np.array, optional
-        Array of shape (n, 3) containing the cartesian coordinates of the apex.
-    X : np.array, optional
-        Array of shape (n, 3) containing the cartesian coordinate points to
-        check.
-
-    Returns
-    -------
-    np.array
-        Array of shape (n,) containing theta angles in degrees.
-
-    Notes
-    -----
-    The `n` parameter is necessary when `angType=0`. The `U` and `X` parameters
-    are necessary when `angtype=1`.
-
-    If the constraint angle type is an altitude angle, then :math:`\theta` can
-    be calculated from :math:`\theta=90\degree-a\degree` where :math:`a` is the
-    altitude constraint angle.
-
-    If the constraint angle type is an off-nadir angle, then :math:`theta` can
-    be calculated from:
-
-    .. math::
-
-        \theta = \gamma+\sin^{-1}\left(\frac{\left|\textbf{X}-\textbf{U}
-        \right|}{\left|\textbf{U}\right|}\sin{\gamma}\right)
-
-    where :math:`\gamma` is the off-nadir constraint angle.
-    """
-
-    if form:
-        ang = np.radians(ang)
-
-        num = np.linalg.norm(X - U, axis=1)
-        denom = np.linalg.norm(U, axis=1)
-
-        theta = ang + np.arcsin(np.sin(ang) * num / denom)
-        theta = np.degrees(theta)
-
-    else:
-        theta = np.full((n,), 90 - ang)
-
-    return theta
-
-
-def _analytical_encounter_ind(sat_position: np.array, gnd_position: np.array,
-                             ang: float, form: Literal[0, 1]) -> np.array:
-    """Return encounter indices.
-
-    This method returns the encounter indices corresponding only to valid
-    satellite positions. It does not consider the sun constraint angle or
-    lighting conditions.
-
-    Parameters
-    ----------
-    sat_position : np.array
-        Array of shape (n, 3) containing satellite ECEF positions.
-    gnd_position : np.array
-        Array of shape (n, 3) containing ground location ECEF positions.
-    ang : float
-        Constraint angle in degrees.
-    angType : {0, 1}
-        Defines the angle as altitude if `angType=0` and off-nadir if
-        `angType=1`.
-
-    Returns
-    -------
-    np.array
-        Array of indices defining valid encounter positions.
-
-    Notes
-    -----
-    This method finds all indices that fall within a single-sided cone
-    extending above the Earth's surface with the apex at a ground location of
-    interest and its aperture angle defined by the constraint angle.
-    """
-
-    n = sat_position.shape[0]
-    U = gnd_position
-    X = sat_position
-    theta = _aperature_theta(ang, form, n, U, X)
-
-    ind_1 = _cone_constraint(theta, U, X)
-    ind_2 = _plane_constraint(U, X)
-
-    ind_final = np.intersect1d(ind_1, ind_2)
-
-    return ind_final
diff --git a/celest/encounter/_window_handling.py b/celest/encounter/_window_handling.py
index 6fbaec6..8539dfb 100644
--- a/celest/encounter/_window_handling.py
+++ b/celest/encounter/_window_handling.py
@@ -1,75 +1,50 @@
-"""Window information localization and handling.
 
-This module contains the `Window` class to isolate window information in
-addition to the `Windows` class to handle `Window` objects in a linked
-list structure.
-"""
 
-
-from typing import Any, Literal
 from celest.encounter.groundposition import GroundPosition
-import pandas as pd
+from typing import Any
 import numpy as np
+import pandas as pd
 
 
-class Window(object):
-    """Isolate window information.
+class Window:
+    """Window(satellite, location, start, end, enc, ang, lighting)
+
+    Encounter information.
 
     Parameters
     ----------
     satellite : Satellite
-        Satellite associated with the encounter window.
+        Satellite for ground location encounter.
     location : GroundPosition
-        Ground location for imaging or ground station communcations.
-    start : float
-        Window start time in Julian.
-    end : float
-        Window end time in Julian.
+        Ground location for satellite encounter.
+    start, end: float
+        Window start and end times in J2000 Julian.
     enc : {"data link", "image"}
-        Type of encounter.
+        Encounter type.
     ang : float
         Encounter constraint angle.
     lighting : {-1, 0, 1}
-        Lighting constraint to specify night only, any time, or day only
-        encounters.
-    sca : float
-        Solar constraint angle.
+        Lighting constraint for night only, all time, or day only encounters.
 
     Attributes
     ----------
     satellite : Satellite
-        Satellite associated with the encounter window.
+        Satellite for ground location encounter.
     coor : tuple
-        Ground location coordinate tuple in the (lat, lon) format.
-    start : float
-        Window start time in Julian.
-    end : float
-        Window end time in Julian.
+        Ground location (lat, lon) coordinates.
+    start, end : float
+        Window start and end times in J2000 Julian.
     dt : float
         Window duration in seconds.
     type : {"data link", "image"}
-        Type of encounter.
-    angle_type : {"A", "N"}
-        Employed constraint angle as either being the off-nadir or altitude
-        angle.
+        Encounter type.
     ang : float
-        Encounter constraint angle.
+        Encounter constraint angle in degrees.
     lighting : {-1, 0, 1}
-        Lighting constraint to specify night only, any time, or day only
-        encounters.
-    sca : float
-        Solar constraint angle.
-
-    Methods
-    -------
-    copy()
-        Copy local window data to a new `Window` object.
+        Lighting constraint for night only, all time, or day only encounters.
     """
 
-    def __init__(self, satellite: Any, location: Any, start: float, end: float,
-                 enc: Literal["data link", "image"], ang: float, lighting:
-                 Literal[-1, 0, 1], sca: float) -> None:
-        """Initialize attributes."""
+    def __init__(self, satellite, location, start, end, enc, ang, lighting) -> None:
 
         self.satellite = satellite
         self.coor = (location.lat, location.lon)
@@ -79,23 +54,19 @@ def __init__(self, satellite: Any, location: Any, start: float, end: float,
         self.duration = 86400 * (end - start)
 
         self.type = enc
-        self.angle_type = "N" if enc == "image" else "A"
         self.ang = ang
         self.lighting = lighting
-        self.sca = sca
 
     def __str__(self) -> str:
-        """Define `Window` information string."""
 
         str_1 = str(self.coor) + " "
         str_2 = self.type + " "
-        str_3 = "ang:" + self.angle_type + str(self.ang) + " "
+        str_3 = "ang:" + str(self.ang) + " "
         str_4 = "lighting:" + str(self.lighting) + " "
-        str_5 = "sca:" + str(self.sca) + " "
-        str_6 = "start:" + str(self.start) + " "
-        str_7 = "end:" + str(self.end)
+        str_5 = "start:" + str(self.start) + " "
+        str_6 = "end:" + str(self.end)
 
-        str_final = str_1 + str_2 + str_3 + str_4 + str_5 + str_6 + str_7
+        str_final = str_1 + str_2 + str_3 + str_4 + str_5 + str_6
 
         return str_final
 
@@ -105,7 +76,6 @@ def copy(self) -> Any:
         Returns
         -------
         Window
-            Return `Window` with the same attributes as `self`.
 
         Examples
         --------
@@ -121,23 +91,26 @@ def copy(self) -> Any:
         enc = self.type
         ang = self.ang
         lighting = self.lighting
-        sca = self.sca
 
-        return_window = Window(sat, gnd, start, end, enc, ang, lighting, sca)
+        return_window = Window(sat, gnd, start, end, enc, ang, lighting)
 
         return return_window
 
 
 class Windows(object):
+    """Data structure to hold `Window` objects.
+
+    This class is designed to hold encounter information and facilitate
+    user-data interactions.
+    """
 
     def __init__(self) -> None:
-        """Initialize attributes."""
 
         self.windows = pd.Series(dtype="object")
         self.index_head = 0
 
     def __iter__(self) -> Any:
-        
+
         return self
 
     def __next__(self) -> Any:
@@ -148,49 +121,82 @@ def __next__(self) -> Any:
         else:
             index = self.windows.index[self.index_head]
             self.index_head += 1
+
             return self.windows[index]
 
+    def __getitem__(self, key) -> Window:
+        """Return window closest to key time.
+
+        This method allows for window indexing. If a scalar index is provided,
+        the window with the closest start time will be returned. If the index
+        is a range, all windows that have start times within the range will be
+        returned. If the index is a tuple, then all the windows with start
+        times closes to each tuple entry will be returned.
+        """
+
+        inds = self.windows.index.to_numpy()
+
+        if isinstance(key, (int, float)):
+
+            ind = np.argmin(np.abs(inds - key))
+            rtn = self.windows[inds[ind]]
+
+        elif isinstance(key, slice):
+
+            rtn = self.windows[key].to_numpy()
+
+        elif isinstance(key, tuple):
+
+            inds_new = []
+            for k in key:
+                ind = np.argmin(np.abs(inds - k))
+                inds_new.append(inds[ind])
+
+            inds_new = list(set(inds_new))
+            rtn = self.windows[inds_new].to_numpy()
+
+        return rtn
+
     def __len__(self) -> int:
-        """Return number of windows."""
 
         return len(self.windows)
 
-    def _add_window(self, window: Window) -> None:
+    def _add_window(self, window) -> None:
         """Add new window to data base.
-        
+
         Parameters
         ----------
         window : Window
-            Window to add to data base.
 
         Notes
         -----
         This method adds a new window to the data base by appending it to the
         internal Pandas Series object indexed by the window's start time.
         """
-        
+
         temp_series = pd.Series(window, index=[window.start])
-        self.windows = self.windows.append(temp_series)
+        self.windows = pd.concat([self.windows, temp_series])
         self.windows.sort_index(inplace=True)
 
     def stats(self) -> pd.DataFrame:
         """Return enocunter statistics.
 
         This method returns statistics on the stored encounters. The statistics
-        include the 5th, 50th, and 95th percentile of encounter durations as
-        well as the enocounter frequency.
-        
+        include 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 `Windows` object returned from the windows
         generation function.
 
-        >>> toronto_img.stats()
+        >>> stats = toronto_img.stats()
+        >>> print(stats)
                                        (45, -73), image
         Q5 Duration (s)                      108.468446
         Q50 Duration (s)                     167.687443
@@ -200,23 +206,24 @@ def stats(self) -> pd.DataFrame:
 
         if len(self.windows) == 0:
             return None
-        
+
         enc_dict = {}
         start_times = np.empty((0,))
         end_times = np.empty((0,))
-        
+
         for window in self.windows:
 
             if (window.coor, window.type) not in enc_dict.keys():
-                enc_dict[window.coor, window.type] = np.array([window.duration])
+                enc_dict[window.coor, window.type] = np.array(
+                    [window.duration])
             else:
                 arr = enc_dict[window.coor, window.type]
                 new_arr = np.append(arr, window.duration)
                 enc_dict[window.coor, window.type] = new_arr
-            
+
             start_times = np.append(start_times, window.start)
             end_times = np.append(end_times, window.end)
-        
+
         out_data = {}
         num_days = np.max(end_times) - np.min(start_times)
 
@@ -228,23 +235,23 @@ def stats(self) -> pd.DataFrame:
             freq = len(dt_data) / num_days
 
             out_data[enc_str] = pd.Series([qs[0], qs[1], qs[2], freq])
-        
+
         df = pd.DataFrame.from_dict(out_data)
         df.index = ["Q5 Duration (s)",
                     "Q50 Duration (s)",
                     "Q75 Duration (s)",
                     "Encounter Frequency (per day)"]
-        
+
         return df
 
-    def to_numpy(self) -> np.array:
+    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
@@ -260,9 +267,9 @@ def to_numpy(self) -> np.array:
 
         return self.windows.to_numpy()
 
-    def save(self, fname: str, delimiter: Literal[",", "\\t"]) -> None:
+    def save(self, fname, delimiter) -> None:
         """Save encounter information.
-        
+
         Parameters
         ----------
         fname : str
@@ -270,7 +277,7 @@ def save(self, fname: str, delimiter: Literal[",", "\\t"]) -> None:
         delimiter : {",", "\\t"}
             String or character separating columns.
         """
-        
+
         np_data = self.to_numpy()
 
         out_data = np.empty(shape=(np_data.shape[0] + 1, 6), dtype="<U25")
@@ -284,7 +291,7 @@ def save(self, fname: str, delimiter: Literal[",", "\\t"]) -> None:
         ]
 
         for i, window in enumerate(np_data):
-            
+
             out_data[i + 1, 0] = window.coor[0]
             out_data[i + 1, 1] = window.coor[1]
             out_data[i + 1, 2] = window.start
@@ -292,4 +299,5 @@ def save(self, fname: str, delimiter: Literal[",", "\\t"]) -> None:
             out_data[i + 1, 4] = window.duration
             out_data[i + 1, 5] = window.type
 
-        np.savetxt(fname, out_data, delimiter=delimiter, fmt="%s,%s,%s,%s,%s,%s")
+        np.savetxt(fname, out_data, delimiter=delimiter,
+                   fmt="%s,%s,%s,%s,%s,%s")
diff --git a/celest/encounter/_window_utils.py b/celest/encounter/_window_utils.py
deleted file mode 100644
index f315a35..0000000
--- a/celest/encounter/_window_utils.py
+++ /dev/null
@@ -1,122 +0,0 @@
-"""Window generation dependencies.
-
-This module contains implementation specific functions used in identifying
-viable encounters.
-"""
-
-from celest.encounter._encounter_math_utils import _analytical_encounter_ind
-from celest.satellite.coordinate import Coordinate
-from celest.satellite.time import Time
-from typing import Any, Literal
-from jplephem.spk import SPK
-import numpy as np
-import pkg_resources
-
-
-def _sun_itrs(julian: np.ndarray) -> np.ndarray:
-    """Return Sun itrs position.
-
-    Parameters
-    ----------
-    julian : np.ndarray
-        Array of shape (n,) containing Julian time data.
-
-    Returns
-    -------
-    np.ndarray
-        Array of shape (n, 3) containing Sun ECEF positions.
-    """
-
-    ephem = pkg_resources.resource_filename(__name__, '../data/de421.bsp')
-    kernal = SPK.open(ephem)
-
-    ssb2sunb = kernal[0, 10].compute(julian)
-    ssb2eb = kernal[0, 3].compute(julian)
-    eb2e = kernal[3, 399].compute(julian)
-    e2sun = (ssb2sunb - ssb2eb - eb2e).T
-
-    SPK.close(kernal)
-
-    sun_itrs = Coordinate(e2sun, "gcrs", Time(julian)).itrs()
-
-    return sun_itrs
-
-def _get_ang(u: np.ndarray, v: np.ndarray) -> np.ndarray:
-    """Calculate degree angle bewteen two vectors.
-
-    Parameters
-    ----------
-    u, v : np.ndarray
-        Arrays of shape (n,3) with rows of XYZ cartesian data.
-
-    Returns
-    -------
-    np.ndarray
-        Array of shape (n,) containing the degree angle between the two
-        arrays.
-    """
-
-    num = np.sum(u * v, axis=1)
-    denom = np.linalg.norm(u, axis=1) * np.linalg.norm(v, axis=1)
-    ang = np.degrees(np.arccos(num / denom))
-
-    return ang
-
-def _window_encounter_ind(satellite: Any, location: Any, ang: float, form:
-                          Literal[0, 1], sca: float=0, lighting:
-                          Literal[-1, 0, 1]=0) -> np.ndarray:
-    """Return indices for possible encounters.
-
-    Parameters
-    ----------
-    satellite : Satellite
-        Satellite involved with a ground location encounter.
-    location : GroundPosition
-        Ground location involved with a satellite encounter.
-    ang : float
-        Encounter constraint angle.
-    form : {0, 1}
-        Angle type used where 0 is for the altitude angle and 1 is for the
-        off-nadir angle.
-    sca : float, optional
-        Solar constraint angle, by default 0
-    lighting : {-1, 0, 1}, optional
-        lighting conditions for encounters where -1 is for night time
-        encounters, 0 is for all-time encoutners, and 1 is for day time
-        encounters, by default 0
-
-    Returns
-    -------
-    np.ndarray
-        Array containing indices defining viable satellite encounter positions.
-    """
-
-    sat_itrs = satellite.position.itrs()
-    time_data = satellite.time.julian()
-
-    ground_GEO = [location.lat, location.lon]
-    ground_GEO = np.repeat(np.array([ground_GEO]), time_data.size, axis=0)
-    gnd_itrs = Coordinate(ground_GEO, "geo", Time(time_data)).itrs()
-
-    enc_params = [sat_itrs, gnd_itrs, ang, form]
-    analytical_ind = _analytical_encounter_ind(*enc_params)
-
-    sun_itrs = _sun_itrs(time_data)
-    sca_angs = _get_ang(sat_itrs - gnd_itrs, sun_itrs - gnd_itrs)
-    sca_ind = np.where(sca_angs > sca)[0]
-
-    enc_ind = np.intersect1d(analytical_ind, sca_ind)
-
-    if lighting != 0:
-        sun_angs = _get_ang(gnd_itrs, sun_itrs)
-
-        if lighting == 1:
-            # Isolate indices occuring at day.
-            day_ind = np.where(sun_angs < 90)
-            enc_ind = np.intersect1d(enc_ind, day_ind)
-        else:
-            # Isolate indices occuring at night.
-            night_ind = np.where(sun_angs >= 90)
-            enc_ind = np.intersect1d(enc_ind, night_ind)
-
-    return enc_ind
diff --git a/celest/encounter/attitudes.py b/celest/encounter/attitudes.py
deleted file mode 100644
index 9680db1..0000000
--- a/celest/encounter/attitudes.py
+++ /dev/null
@@ -1,197 +0,0 @@
-"""Generate pointing profiles for a satellite.
-
-This module contains the functionality to determine satellite orientations for
-ground tracking.
-"""
-
-
-from celest.satellite.coordinate import Coordinate
-from scipy.spatial.transform import Rotation, Slerp
-from typing import Any
-import numpy as np
-
-
-
-def skew(matrix: np.ndarray) -> np.ndarray:
-    """Return the skew of a 3x3 input matrix.
-
-    Parameters
-    ----------
-    matrix : np.ndarray
-        The 3x3 matrix to skew.
-    
-    Returns
-    -------
-    np.ndarray
-        The skewed version of `matrix`.
-    
-    Notes
-    -----
-    The skew of a column vector given by
-    :math:`\vec{V}=\langle a_x, a_y, a_z\rangle^T` can be found as follows:
-
-    .. math::
-    
-        \vec{v}_{skew} = \begin{bmatrix}
-                            0 & -a_z & a_y \\
-                            a_z & 0 & -a_x \\
-                            -a_y & a_x & 0 \\
-                         \end{bmatrix}
-    """
-
-    n = matrix.shape[0]
-    zeros = np.zeros(n)
-    a_x, a_y, a_z = matrix[:, 0], matrix[:, 1], matrix[:, 2]
-
-    matrix_skewed = np.array([[zeros, -a_z, a_y],
-                              [a_z, zeros, -a_x],
-                              [-a_y, a_x, zeros]])
-
-    return np.transpose(matrix_skewed, axes=(2, 0, 1))
-
-
-def sat_rotation(sat_z: np.ndarray) -> Rotation:
-    """Return satellite rotations for encounter alignment.
-
-    This function defines the satellite orientation as a rotation to be applied
-    to the +Z gcrs axis to align with the Satellites +Z axis.
-
-    Parameters
-    ----------
-    z_sat: np.ndarray
-        Array of shape (n, 3) containing the Satellite's +Z vectors in the gcrs
-        frame.
-    
-    Returns
-    -------
-    Rotation
-        The rotations defined to transform the gcrs z-axis to the Satellite's
-        z-axis.
-    
-    Notes
-    -----
-    Let the gcrs z-axis be :math:`\vec{a}` and the satellite's z-axis be
-    :math:`\vec{b}`, then the rotation matrix, :math:`R`, is defined as follows
-
-    .. math:: R = I + [\vec{v}]_{\times} + [\vec{v}]_{\times}^2\frac{1 - c}{s^2}
-
-    where :math:`I` is the :math:`3x3` identity matrix,
-    :math:`[\vec{v}]_{\times}` is the skew of
-    :math:`\vec{v}=\vec{a}\times\vec{b},
-    :math:`s = \left|\left|\vec{v}\right|\right|`, and
-    :math:`c=\vec{a}\cdot\vec{b}`.[1]
-
-    References
-    ----------
-    .. [1] “Calculate rotation matrix to align vector A to vector B in 3d?,”
-       Stackexchange.com. [Online]. Available:
-       https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d.
-       [Accessed: 29-Jul-2021].
-    """
-
-    # Reference upward direction.
-    z_ECI = np.array([0, 0, 1])
-
-    # Normalize the vectors.
-    normalized_sat_z = sat_z / np.linalg.norm(sat_z, axis=1)[:, np.newaxis]
-
-    # Cross of difference between the unit vectors.
-    v = np.cross(z_ECI[np.newaxis, :], normalized_sat_z, axis=1)
-
-    sin = np.linalg.norm(v, axis=1)  # Sin of angle.
-    cos = np.sum(z_ECI[np.newaxis, :] * normalized_sat_z, axis=1)  # Cos of angle.
-
-    identity = np.eye(3, 3)[np.newaxis, :, :]
-    v_skew = skew(v)
-    v_skew_squared = np.matmul(v_skew, v_skew)
-    ratio = ((1 - cos) / sin**2)[:, np.newaxis, np.newaxis]
-
-    matrix = identity + v_skew + v_skew_squared * ratio
-
-    return Rotation.from_matrix(matrix)
-
-
-def attitudes(satellite: Any, location: Any, encInd: np.ndarray, maneuverTime: int) -> Rotation:
-    """Generate satellite rotations for ground tracking.
-
-    This function takes a ground location along and indices associated with
-    delested windows at which the spacecraft makes imaging passes over the
-    location. It then generates the satellite rotations to track the ground
-    location.
-
-    Parameters
-    ----------
-    location : GroundPosition
-        Ground location and encounter information.
-    encInd: np.ndarray
-        Array of indices that correspond to times and positions where the
-        spacecraft is in an imaging encounter window with the given ground
-        location.
-    maneuverTime: float
-        Number of array indices to pad on either side of
-        an encounter window to use for maneuvering time.
-
-    Notes
-    -----
-    The strategy for pointing profile generation is as follows: (1) The default
-    orientation is to have the spacecraft camera pointing towards its zenith.
-    (2) When the spacecraft is imaging, orient the satellite such that the
-    camera is facing the target. (3) Generate an initial set of pointing
-    profiles assuming the above. (4) Interpolate the rotations between the
-    normal and target-acquired states to smooth out transitions.
-    """
-
-    ground_GEO = [location.lat, location.lon, location.radius]
-    ground_GEO = np.repeat(np.array([ground_GEO]), satellite.position.length, axis=0)
-    target_site = Coordinate(ground_GEO, "geo", satellite.time)
-
-    # Stage 1: preliminary rotations.
-    # Get difference vector between spacecraft and target site.
-    SC_to_site = target_site.gcrs() - satellite.position.gcrs()
-
-    # Point toward the zenith except when over the imaging site.
-    pointing_directions = satellite.position.gcrs()
-    pointing_directions[encInd, :] = SC_to_site[encInd, :]
-
-    # Preliminary rotation set.
-    # Temporarily represent as quaternion for interpolation.
-    rotations = sat_rotation(pointing_directions).as_quat()
-
-    # Set flight modes. By default, point normal.
-    flight_ind = 0 * np.ones(SC_to_site.shape[0])
-
-    # Point to target during encounters.
-    flight_ind[encInd] = 2
-
-    # Stage 2: interpolation.
-    # Generate sets of encounters which are clustered together.
-    # This takes individual encInd into clusters which we can use later to
-    # figure out when to start interpolation.
-    split_ind = np.where(np.diff(encInd) > 1)[0]+1
-    print(split_ind)
-    encounter_segments = np.split(encInd, split_ind)
-    print(encounter_segments)
-
-    for seg in encounter_segments:
-        print(seg)
-
-        start_step = seg[0] - maneuverTime
-        end_step = seg[-1] + maneuverTime
-
-        # Get starting and ending quaternions.
-        start_rotation = rotations[start_step]
-        end_rotation = rotations[end_step]
-
-        slerp_1 = Slerp([start_step, seg[0]], Rotation.from_quat([start_rotation, rotations[seg[0]]]))
-        interp_rotations_1 = slerp_1(np.arange(start_step, seg[0]))
-
-        rotations[start_step:seg[0]] = interp_rotations_1.as_quat()
-        flight_ind[start_step:seg[0]] = 1
-
-        slerp_2 = Slerp([seg[-1], end_step], Rotation.from_quat([rotations[seg[-1]], end_rotation]))
-        interp_rotations_2 = slerp_2(np.arange(seg[-1], end_step))
-
-        rotations[seg[-1]:end_step] = interp_rotations_2.as_quat()
-        flight_ind[seg[-1]:end_step] = 3
-
-    return Rotation.from_quat(rotations)
diff --git a/celest/encounter/groundposition.py b/celest/encounter/groundposition.py
index 6e43f1c..08ef5a5 100644
--- a/celest/encounter/groundposition.py
+++ b/celest/encounter/groundposition.py
@@ -1,78 +1,55 @@
-"""Localize ground position information.
 
-This module contains the `GoundPosition` class which stores ground location
-data and is an instrumental part of encounter planning by managing the ground
-location and satellite encounter relationship.
-"""
 
-
-from celest.core.decorators import set_module
 import numpy as np
 
 
-@set_module('celest.encounter')
 class GroundPosition(object):
-    """Localize Earth bound location information.
+    """GroundPosition(latitude, longitude)
 
-    The `GoundPosition` class stores ground location information and is an
-    instrumental part of encounter planning functionality.
+    Localize Earth surface location information.
 
     Parameters
     ----------
     latitude : float
-        Locations geodetic lattitude in degrees and decimals.
+        Lattitude of the location in decimal degrees.
     longitude : float
-        Locations geodetic longitude in degrees and decimals.
+        Longitude of the location in decimal degrees.
 
     Attributes
     ----------
-    lat : float
-        Geodetic latitude for the location given in degrees and decimals.
-    lon : float
-        Geodetic longitude for the location given in degrees and decimals.
+    lat, lon : float
+        Latitude and longitude of the location in decimal degrees.
     radius : float
-        Earth's radius at the given coordinates.
+        Earth radius at (`latitude`, `longitude`).
     """
 
     def __init__(self, latitude: float, longitude: float) -> None:
-        """Initialize attributes."""
 
         self.lat = latitude
         self.lon = longitude
         self.radius = self._radius(latitude)
 
     def __str__(self) -> str:
-        """Define information string."""
-
-        output = [
-            f"coor={(self.lat, self.lon)}",
-            f"radius={round(self.radius, 5)}"
-        ]
 
-        return ", ".join(output)
+        return f"coor={(self.lat, self.lon)}, radius={round(self.radius, 5)}"
 
     def _radius(self, latitude: float) -> float:
-        """Calculate the Earth's geocentric radius using WGS84.
+        """Calculate geocentric radius using WGS84.
 
         Parameters
         ----------
         latitude : float
-            Latitude of the ground location in degrees and decimals.
+            Location's latitude in decimal degrees.
 
         Returns
         -------
         float
-            Earth's geocentric radius in kilometres and decimals.
+            Earth's geocentric radius in decimal kilometers.
 
         Notes
         -----
-        By using an Earth ellipsoid with the WGS84 parameters of
-        :math:`a=6378.137` and :math:`b=6356.7523142`, the geocentric radius
-        can be calculated using the following formulation:
-
-        .. math:: r = \sqrt{\frac{(a^2\cos(\beta))^2 + (b^2\sin(\beta))^2}{(a\cos(\beta))^2 + (b\sin(\beta))^2}}
-
-        where :math:`\beta` is the observer's latitude.[1]_
+        The WGS84 Earth ellipsoid model is used as discussed in "Earth Radius
+        by Latitude (WGS 84)" by Timur. [1]_
 
         References
         ----------
diff --git a/celest/encounter/windows.py b/celest/encounter/windows.py
index 43d41ad..0b24d03 100644
--- a/celest/encounter/windows.py
+++ b/celest/encounter/windows.py
@@ -1,126 +1,211 @@
-"""Generate window sets for specified encounters.
 
-This module contains functions to generate window sets for various satellite
-and ground-location encounters.
-"""
 
-
-from celest.core.decorators import set_module
-from celest.core.interpolation import _interpolate
 from celest.encounter._window_handling import Window, Windows
-from celest.encounter._window_utils import _window_encounter_ind
-from typing import Any, Literal
+from celest.satellite.coordinate import Coordinate
+from jplephem.spk import SPK
 import numpy as np
+import pkg_resources
 
 
-_image_encounter = {
-    "type": "I",
-    "angType": 1,  # Off-nadir angle.
-    "lighting": 1,
-    "sca": 0
-}
+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.
 
-_data_link_encounter = {
-    "type": "T",
-    "angType": 0,  # Altitude angle.
-    "lighting": 0,
-    "sca": 30
-}
+    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.
+    """
 
-@set_module('celest.encounter.windows')
-def generate(satellite: Any, location: Any, enc: Literal["image", "data link"],
-             ang: float, factor: int=5) -> Windows:
-    """Return encounter windows.
+    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 _sun_coor(julian):
+    """Return the sun's coordinates.
 
     Parameters
     ----------
-    satellite : Satellite
-        Satellite taking part in ground interactions.
-    location : GroundPosition
-        Ground location of imaging site or ground station.
-    enc : {"image", "data link"}
-        Type of encounter as being an imaging or data link encounter.
-    ang : float
-        Encounter contraint angle in degrees.
-    factor : int, optional
-        Data interpolation factor for more precise encounter information, by
-        default 5.
+    julian : np.ndarray
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
-    Windows
-        The encounter opportunities between the satellite and ground location
-        of the type specified.
-    
-    Notes
-    -----
-    Imaging encounters are use the off-nadir angle, measured in increasing
-    degrees from the satellite's nadir to the ground location. When the
-    off-nadir angle is used, the input `ang` provides a maximum constraint.
-    Data link encounters use the altitude angle, measured in increasing
-    degrees from the ground location's horizon to the satellite. When the
-    altitude angle is used, the input `ang` provides a minimum constraint.
-
-    Examples
-    --------
-    >>> toronto = GroundPosition(latitude=43.65, longitude=-79.38)
-    >>> toronto_dl = windows.generate(satellite, toronto, "data link", 30)
+    Coordinate
+        Coordinate object containing the sun's position.
     """
 
-    windows = Windows()
+    ephem = pkg_resources.resource_filename(__name__, '../data/de421.bsp')
+    kernal = SPK.open(ephem)
 
-    if enc == "image":
-        ang_type = _image_encounter["angType"]
-        lighting = _image_encounter["lighting"]
-        sca = _image_encounter["sca"]
-    else:
-        ang_type = _data_link_encounter["angType"]
-        lighting = _data_link_encounter["lighting"]
-        sca = _data_link_encounter["sca"]
+    ssb2sunb = kernal[0, 10].compute(julian)
+    ssb2eb = kernal[0, 3].compute(julian)
+    eb2e = kernal[3, 399].compute(julian)
+    e2sun = (ssb2sunb - ssb2eb - eb2e).T
+
+    SPK.close(kernal)
 
-    if factor > 1:
+    sun_coor = Coordinate(e2sun, "gcrs", julian)
 
-        if ang_type:
-            off_nadir = satellite.position.off_nadir(location)
-            ind = np.where(off_nadir < ang)[0]
+    return sun_coor
 
-        elif not ang_type:
-            alt, _ = satellite.position.horizontal(location)
-            ind = np.where(alt > ang)[0]
 
-        ind = np.split(ind, np.where(np.diff(ind) != 1)[0] + 1)
-        ind = np.array(ind, dtype=object)
+def _root_find(f, tl, tr, tol):
+    """Return root of f(t) between lg and rg.
+
+    Parameters
+    ----------
+    f : Stroke
+        Stroke object representing the function to be evaluated.
+    tl : float
+        Lower bound of the search interval.
+    tr : float
+        Upper bound of the search interval.
+    tol : float
+        Search tolerance.
 
-        julian_interp = _interpolate(satellite.time.julian(), factor, 2, ind)
-        eci_interp = _interpolate(satellite.position.gcrs(), factor, 2, ind)
+    Returns
+    -------
+    float
+        Root of `f(t)` between `lg` and `rg`.
+    """
 
-        satellite.position._GCRS = eci_interp
-        satellite.position._ITRS = None
-        satellite.position._GEO = None
-        satellite.position.length = eci_interp.shape[0]
-        satellite.position.time._julian = julian_interp
-        satellite.position.time._length = julian_interp.size
-        satellite.time = satellite.position.time
+    l, r = tl, tr
+    fl = int(f(l))
 
-    enc_ind = _window_encounter_ind(satellite, location, ang, ang_type, sca, lighting)
+    while abs(r - l) > tol:
 
-    window_ind = np.split(enc_ind, np.where(np.diff(enc_ind) != 1)[0] + 1)
-    window_ind = np.array(window_ind, dtype=object)
+        c = (l + r) / 2
+        fc = int(f(c))
+
+        if abs(fl - fc) == 1:
+            r = c
+        else:
+            l = c
+            fl = fc
+
+    return (l if fl == 1 else r)
+
+
+def _get_ang(u: np.ndarray, v: np.ndarray) -> np.ndarray:
+    """Calculate degree angle bewteen two vectors.
+
+    Parameters
+    ----------
+    u, v : np.ndarray
+        2-D arrays containing row-vectors.
+
+    Returns
+    -------
+    np.ndarray
+        1-D array containing the degree angle between rows of `u` and `v`.
+    """
+
+    ua = None if u.ndim == 1 else 1
+    va = None if v.ndim == 1 else 1
+
+    n = np.sum(u * v, axis=(ua or va))
+    d = np.linalg.norm(u, axis=ua) * np.linalg.norm(v, axis=va)
+    ang = np.degrees(np.arccos(n / d))
+
+    return ang
+
+
+def generate(satellite, location, enc, ang, lighting=0, tol=1e-5):
+    """Return windows for a given encounter.
+
+    Parameters
+    ----------
+    satellite : Satellite
+        Satellite object.
+    location : Location
+        Location of ground station.
+    enc : {"image", "data_link"}
+        Type of encounter.
+    ang : float
+        Constraint angle 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
+    -------
+    Windows
+        The satellite window details.
+    """
+
+    # Determine windows based on the type of encounter.
+    if enc == "image":
+        dv = satellite._altitude(location)
+        raw_windows_1 = _constraint_highlighter(dv, "ge", 0)
+
+        dv = satellite.off_nadir(location, stroke=True)
+        raw_windows_2 = _constraint_highlighter(dv, "lt", ang)
+
+        raw_windows = raw_windows_1 * raw_windows_2
+    elif enc == "data_link":
+        dv = satellite._altitude(location)
+        raw_windows = _constraint_highlighter(dv, "gt", ang)
+    else:
+        raise ValueError("Invalid encounter type.")
+
+    if lighting != 0:
+        sun_coor = _sun_coor(satellite._julian)
+        sun_alt, _ = sun_coor.horizontal(location, stroke=True)
+        comp = "gt" if lighting == 1 else "le"
+        raw_windows = raw_windows * _constraint_highlighter(sun_alt, comp, 0)
+
+    julian = satellite._julian
+    ind = np.arange(0, julian.size, 1) * raw_windows(julian).astype(int)
+    ind = ind[ind != 0]
+    ind = np.split(ind, np.where(np.diff(ind) != 1)[0] + 1)
+    ind = np.array(ind, dtype=object)
+
+    windows = Windows()
 
-    times = satellite.time.julian()
+    # Populate Windows object.
+    if ind.size == 0:
+        return windows
 
-    if window_ind.size != 0:
+    for i in ind:
 
-        n = len(window_ind)
+        st1, st2 = julian[i[0] - 1], julian[i[0]]
+        start = _root_find(raw_windows, st1, st2, tol)
 
-        for j in range(n):
+        et1, et2 = julian[i[-1]], julian[i[-1] + 1]
+        end = _root_find(raw_windows, et1, et2, tol)
 
-            start = times[window_ind[j][0]]
-            end = times[window_ind[j][-1]]
-            window = Window(satellite, location, start, end, enc, ang, lighting, sca)
+        window = Window(satellite, location, start, end, enc, ang, lighting)
+        windows._add_window(window)
 
-            windows._add_window(window)
-            
     return windows
diff --git a/celest/satellite/_angle_representations.py b/celest/satellite/_angle_representations.py
index a349458..bae9db8 100644
--- a/celest/satellite/_angle_representations.py
+++ b/celest/satellite/_angle_representations.py
@@ -1,23 +1,20 @@
-"""Define various angle representations."""
 
 
-from celest.core.decorators import set_module
 import numpy as np
 
 
-@set_module('celest.satellite')
 def sexagesimal(angles: np.ndarray) -> np.ndarray:
     """Convert decimal angles into sexagesimal angles.
 
     Parameters
     ----------
     angles : np.ndarray
-        Array of shape (n,) containing angles in decimal degrees.
+        1-D array containing angles in decimal degrees.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing sexagesimal angles as strings.
+        1-D array containing sexagesimal angles as strings.
 
     Examples
     --------
@@ -30,14 +27,12 @@ def sexagesimal(angles: np.ndarray) -> np.ndarray:
               '+149°06′00.00″'])
     """
 
-    deg = u"\u00B0"
-    min = u"\u2032"
-    sec = u"\u2033"
+    deg, min, sec = u"\u00B0", u"\u2032", u"\u2033"
 
-    length = angles.shape[0]
-    out_arr = np.empty((length,), dtype="<U32")
+    n = angles.shape[0]
+    out_arr = np.empty((n,), dtype="<U32")
 
-    for i in range(length):
+    for i in range(n):
 
         ang = angles[i]
 
@@ -56,27 +51,20 @@ def sexagesimal(angles: np.ndarray) -> np.ndarray:
 
 
 def _ISO6709_representation(position: np.ndarray) -> np.ndarray:
-    """Format geographical data in an international standard.
-
-    This method formats the input geographical point location data in the
-    ISO6709 standard format.
+    """Format geographical data in the ISO6709 standard.
 
     Parameters
     ----------
     position : np.ndarray
-        Array of shape (n, 3) with columns of geodetic latitude,
-        terrestrial longitude, and geodetic altitude data.
+        2-D array containing columns of latitude, longitude, and altitude data.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing position strings in accordance to
-        ISO6709 standard.
+        1-D array containing ISO6709 standardized position strings.
     """
 
-    deg = "\u00B0"
-    min = "\u2032"
-    sec = "\u2033"
+    deg, min, sec = "\u00B0", "\u2032", "\u2033"
 
     out_arr = np.empty((position.shape[0],), dtype="<U37")
 
diff --git a/celest/satellite/_astronomical_quantities.py b/celest/satellite/_astronomical_quantities.py
index 7ae284d..b1623c7 100644
--- a/celest/satellite/_astronomical_quantities.py
+++ b/celest/satellite/_astronomical_quantities.py
@@ -1,8 +1,3 @@
-"""Astronomical quantities for astronomical calculations.
-
-This module contains the functionality to cary out advanced astronomical
-calculations.
-"""
 
 
 from typing import Tuple
@@ -10,10 +5,10 @@
 
 
 def nutation_angles(julian: np.ndarray) -> Tuple:
-    """Return the five Earth nutation angles.
+    """Return five Earth nutation angles.
 
-    This method takes a Julian day array and calculates the five nutation
-    angles at each time. The calculated angles include the mean elongation of
+    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
@@ -22,14 +17,13 @@ def nutation_angles(julian: np.ndarray) -> Tuple:
     Parameters
     ----------
     jullian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     tuple
-        Returns a tuple, `(D, M, N, F, O)`, where the items are Numpy arrays
-        of shape (n,) containing the Earth nutation angles in degrees and
-        decimals.
+        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
     -----
@@ -38,7 +32,7 @@ def nutation_angles(julian: np.ndarray) -> Tuple:
 
     .. math:: T = \\frac{JD - 2451545}{36525}
 
-    The nutation angles can then be calculated in degrees and decimals. [Mee98b]_
+    The nutation angles can then be calculated in decimal degrees. [Mee98b]_
 
     .. math:: D = 297.85036 + 445267.111480T - 0.0019142T^2 + T^3 / 189474
 
@@ -53,7 +47,7 @@ def nutation_angles(julian: np.ndarray) -> Tuple:
     References
     ----------
     .. [Mee98b] Jean Meeus. Astronomical algorithms. 2nd ed. Willmann-Bell,
-       1998, pp. 143–144. isbn: 9780943396613.
+       1998, pp. 143-144. isbn: 9780943396613.
 
     Examples
     --------
@@ -72,52 +66,31 @@ def nutation_angles(julian: np.ndarray) -> Tuple:
 
     return D, M, N, F, O
 
+
 def nutation_components(julian: np.ndarray) -> Tuple:
     """Return the nutations in longitude and obliquity.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     tuple
-        Returns a tuple, `(longitude, obliquity)`, where the items are NumPy
-        arrays of shape (n,) containing the nutation components of longitude
-        and obliquity in seconds and decimals.
+        Tuple of the form, `(longitude, obliquity)`, containing 1-D NumPy
+        arrays with the nutation components of longitude and obliquity in
+        decimal seconds.
 
     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 longitude of the ascending node of the Moon's orbit on the ecliptic
-    measured from the mean equinox can be calculated in decimal degrees
-    using the following:
-
-    .. math:: \Omega = 125.04452 - 1934.136261T + 0.0020708T^2 + T^3 / 450000
-
-    The mean longitudes of the Sun, :math:`L`, and Moon, :math:`L'`, can be
-    calculated in decimal degrees using the following:
-
-    .. math:: L = 280.4665 + 36000.7698 * T
-
-    .. math:: L' = 218.3165 + 481267.8813 * T
-
-    The nutation in longitude, :math:`\Delta\psi`, and nutation in obliquity,
-    :math:`\Delta\epsilon`, can then be calculated in decimal minutes. [Mee98c]_
-
-    .. math:: \Delta\Psi = -17.20\sin\Omega - 1.32\sin 2L - 0.23\sin 2L' + 0.21\sin 2\Omega
-
-    .. math:: \Delta\epsilon = 9.20\cos\Omega + 0.57\cos 2L + 0.10\cos 2L' - 0.09\cos 2\Omega
+    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.
+       1998, pp. 143-144. isbn: 9780943396613.
 
     Examples
     --------
@@ -149,42 +122,73 @@ def nutation_components(julian: np.ndarray) -> Tuple:
 
     return delta_psi, delta_epsilon
 
+
+def precession_angles(julian: np.ndarray) -> Tuple:
+    """Return 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.
+    """
+
+    # Get time constants.
+    t = (julian - 2451545.0) / 36525
+    t2 = t * t
+    t3 = t2 * t
+    t4 = t3 * t
+    t5 = t4 * t
+
+    # Get precession angles.
+    zeta = 2.59796176 + 2306.0809506 * t + 0.3019015 * t2 + 0.0179663 * t3 - \
+        0.0000327 * t4 - 0.0000002 * t5
+    theta = 2004.1917476 * t - 0.4269353 * t2 - 0.0418251 * t3 - \
+        0.0000601 * t4 - 0.0000001 * t5
+    z = - 2.5976176 + 2306.0803226 * t + 1.0947790 * t2 + 0.0182273 * t3 + \
+        0.0000470 * t4 - 0.0000003 * t5
+
+    return zeta, theta, z
+
+
 def mean_obliquity(julian: np.ndarray) -> np.ndarray:
-    """Return the mean obliquity of the ecliptic in the J2000 epoch.
+    """Return the ecliptic's mean obliquity in the J2000 epoch.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Returns an array of shape (n,) containing the mean obliquity of the
-        ecliptic in degrees and decimals.
+        1-D array containing the ecliptic's mean obliquity in degrees and
+        decimals.
 
     See Also
     --------
     apparent_obliquity :
-        Return the apparent obliquity of the ecliptic in the J2000 epoch.
+        Return the ecliptic's apparent obliquity in the J2000 epoch.
 
     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}
-
-    If we define :math:`U = T / 100`, then the mean obliquity of the ecliptic.
-    :math:`\epsilon_o` can be found by the following:
-
-    .. math:: \epsilon_0 = 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
-
-    The algorithm used is only valid for 10000 years on either side of
-    J2000. [Mee98d]_
+    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
     ----------
@@ -216,34 +220,31 @@ def mean_obliquity(julian: np.ndarray) -> np.ndarray:
 
     return epsilon_0
 
+
 def apparent_obliquity(julian: np.ndarray) -> np.ndarray:
-    """Return the apparent obliquity of the ecliptic in the J2000 epoch.
+    """Return the ecliptic's apparent obliquity in the J2000 epoch.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Returns an array of shape (n,) containing the apparent obliquity of
-        the ecliptic in degrees and decimals.
+        1-D array containing the ecliptic's apparent obliquity in degrees and
+        decimals.
 
     See Also
     --------
     mean_obliquity :
-        Return the mean obliquity of the ecliptic in the J2000 epoch.
+        Return the ecliptic's mean obliquity in the J2000 epoch.
 
     Notes
     -----
-    The apparent obliquity of the ecliptic, :math:`\epsilon` can be calculated
-    as :math:`\epsilon = \epsilon_0 + \Delta\epsilon` where :math:`\epsilon_0`
-    is the mean obliquity of the ecliptic and :math:`\Delta\epsilon` is the
-    nutation of obliquity. [Mee98e]_
-
-    The due to the limitations of the algorithm used in the calculation of
-    :math:`\epsilon_0` is only valid for 10000 years on either side of J2000.
+    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
     ----------
@@ -262,20 +263,20 @@ def apparent_obliquity(julian: np.ndarray) -> np.ndarray:
 
     return epsilon
 
+
 def from_julian(julian: np.ndarray) -> Tuple:
-    """Return year, month, day from Julian day data.
+    """Return year, month, day from Julian times.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     Tuple
-        Returns a tuple, `(year, month, day)`, where the items are NumPy
-        arrays of shape (n,) containing the year, month, and decimal days of
-        the input Julian data.
+        Returns a tuple, `(year, month, day)`, containing 1-D NumPy arrays
+        representing the year, month, and decimal days of the input times.
 
     Notes
     -----
@@ -309,40 +310,37 @@ def from_julian(julian: np.ndarray) -> Tuple:
     day = B - D - (30.6001 * E).astype(int) + F
 
     month = E
-    ind_1 = np.where(month < 14)[0]
-    ind_2 = np.where((month == 14) | (month == 15))[0]
-    month[ind_1] = month[ind_1] - 1
-    month[ind_2] = month[ind_2] - 13
+    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
-    ind_1 = np.where(month > 2)[0]
-    ind_2 = np.where((month == 1) | (month == 2))[0]
-    year[ind_1] = year[ind_1] - 4716
-    year[ind_2] = year[ind_2] - 4715
+    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
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing the day of the year.
+        1-D array containing the day of the year.
 
     Notes
     -----
-    The day of the year, :math:`N` can be calculated from the following:
-
-    .. math:: N = INT\left(\frac{275M}{9}\right) - K\times\left(\frac{M+9}{12}\right) + D - 30
-
-    where :math:`M` is the month number, :math:`D` is the day of the month,
-    and :math:`K=1` if the year is a leap year otherwise :math:`K=2`. [Mee98g]_
+    The day of year is calculated using the methods in "Astronomical
+    Algorithms" by Jean Meeus. [Mee98g]_
 
     References
     ----------
@@ -366,29 +364,24 @@ def day_of_year(julian: np.ndarray) -> np.ndarray:
 
     return N.astype(int)
 
+
 def equation_of_time(julian: np.ndarray) -> np.ndarray:
-    """Return the Equation of Time in degrees.
+    """Return the equation of time in decimal degrees.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing the Equation of Time in decimal degrees.
+        1-D array containing the Equation of Time in decimal degrees.
 
     Notes
     -----
-    The equation of time in radians can be calculated from the following:
-
-    .. math:: E = y\sin 2L_0 - 2e\sin M + 4ey\sin M\cos 2L_0 - \frac{1}{2}y^2\sin 4L_0 - \frac{5}{4}e^2\sin 2M
-
-    where :math:`y=\tan^2\frac{\epsilon}{2}`, :math:`\epsilon` is the apparent
-    obliquity of the ecliptic, :math:`L_0` is the Sun's mean longitude,
-    :math:`e` is the eccentricity of the Earth's orbit, and :math:`M` is the
-    Sun's mean anomaly. [Mee98h]_
+    The equation of time is calculated using the methods described in
+    "Astronomical Algorithms" by Jean Meeus. [Mee98h]_
 
     References
     ----------
@@ -420,19 +413,19 @@ def equation_of_time(julian: np.ndarray) -> np.ndarray:
 
     return EOT
 
+
 def equation_of_equinoxes(julian: np.ndarray) -> np.ndarray:
-    """Return the equation of the equinoxes in arcseconds.
+    """Return the equation of the equinoxes in decimal arcseconds.
 
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing the Equation of Equinoxes in decimal
-        arcseconds.
+        1-D array containing the equation of equinoxes in decimal arcseconds.
 
     Notes
     -----
@@ -469,34 +462,18 @@ def sun_right_ascension(julian: np.ndarray) -> np.ndarray:
     Parameters
     ----------
     julian : np.ndarray
-        Array of shape (n,) containing Julian times in decimal days.
+        1-D array containing Julian times in the J2000 epoch.
 
     Returns
     -------
     np.ndarray
-        Array of shape (n,) containing the right ascension of the mean sun in
-        degrees and decimals.
+        1-D array containing containing the right ascension of the mean sun in
+        decimal degrees.
 
     Notes
     -----
-    To calculate the right ascension of the Sun, we must first calculate the
-    following dependencies:
-
-    .. math:: L0 = 280^\circ.46646 + 36000^\circ.76983t_U + 0^\circ.0003032t_U^2
-    .. math:: M = 357^\circ.52911 + 35999^\circ.05029t_U - 0^\circ.0001537t_U^2
-
-    .. math:: C = (1^\circ.914602 - 0^\circ.004817t_U - 0^\circ.000014t_U^2)\sin(m)
-                + (0^\circ.019993 - 0^\circ.000101t_U)\sin(2m)
-                + 0^\circ.000289\sin(3m)
-
-    \odot = L0 + C
-
-    where :math:`t_U = \frac{JD-2451545.0}{36525}` and `JD` is the Julian day.
-    We can then calculate the right ascension, :math:`\alpha`, as
-
-    .. math:: \tan\alpha = \frac{\cos\epsilon\sin\odot}{\cos\odot}
-
-    where :math:`\epsilon` is the mean obliquity of the ecliptic. [Mee98j]_
+    The right ascension of the Sun is calculated using the methods discussed in
+    "Astronomical Algorithms" by Jean Meeus. [Mee98j]_
 
     References
     ----------
diff --git a/celest/satellite/coordinate.py b/celest/satellite/coordinate.py
index 8bc016e..be4eb0d 100644
--- a/celest/satellite/coordinate.py
+++ b/celest/satellite/coordinate.py
@@ -1,293 +1,303 @@
-"""Localize position information representations.
 
-This module provides the `Coordinate` class to allow for a simple interface
-for converting a position type into various representations. The class is also
-used for position inputs and outputs for other `Celest` functionality.
-"""
 
-
-from celest.core.decorators import set_module
 from celest.satellite._angle_representations import _ISO6709_representation
-from typing import Any, Literal, Tuple
+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
 import numpy as np
 
 
-@set_module('celest.satellite')
-class Coordinate(object):
-    """Localize position information representations.
+class Coordinate(Time):
+    """Coordinate(position, frame, julian, offset=0)
+
+    Satellite coordinate transformations.
+
+    The input position and time can be converted into different coordinate
+    representations useful for satellite applications. Here, `julian + offset`
+    is the Julian time in J2000 epoch and is used in various transformations.
 
-    The `Coordinate` class provides a simple user interface for converting a
-    position type into various representations. The class is also used for
-    position inputs and outputs for other `Celest` functionality.
+    If `frame=="geo"`, the `position` input can have two columns, (latitude,
+    longitude) or three columns (latitude, longitude, altitude). If no
+    geodetic altitude is provided, it is assumed to be zero. Other frames
+    require three columns.
 
     Parameters
     ----------
-    position : np.ndarray
-        Base position to initialize the `Coodinate` class.
+    position : array_like
+        2-D array containing position coordinates. Axis 0 is the time series
+        axis.
+
+        Supported coordinate frames include the Geocentric Celestial Reference
+        System (gcrs), International Terrestrial Reference System (itrs), and
+        geographical (geo) system.
     frame : {"gcrs", "geo", "itrs"}
-        Specifies the input position frame.
-    time : Time
-        Times associated with the position data in the J2000 epoch. The length
-        of the `time` parameter must match the length of the `position`
-        parameter.
+        Frame specifier for the input position.
+    julian : array_like
+        1-D array containing time data in Julian days.
 
-    Attributes
-    ----------
-    time : Time
-        Times corresponding to the position data.
-    length : int
-        Length of the input position and time arrays.
-
-    Methods
-    -------
-    geo(iso=False)
-        Return geographical position data.
-    era()
-        Return the Earth rotation angles in radians and decimals.
-    gcrs()
-        Return cartesian gcrs position data.
-    itrs()
-        Return cartesian itrs position data.
-    horizontal(location)
-        Return horizontal position data in degrees and decimals.
-    off_nadir(location)
-        Return the off-nadir angle to a ground location.
-    altitude()
-        Return the altitude above Earth's surface in kilometres.
-    distance(location)
-        Return the distance to a ground location.
+        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 positions:
+
+    >>> julian = [30462.50, 30462.50]
+    >>> position = [[-4681.50824149 -5030.09119386     0.        ]
+    ...             [-4714.35351825 -4978.74325953   454.41492765]]
+    >>> c = Coordinate(position=position, frame="gcrs", julian=julian, offset=0)
+
+    Get horizontal coordinates for the ground location at (43.65, -79.38):
+
+    >>> location = GroundPosition(latitude=43.65, longitude=-79.38)
+    >>> altitude, azimuth = c.horizontal(location=location)
     """
 
-    def __init__(self, position: np.ndarray, frame: Literal["gcrs", "geo", "itrs"], time: Any) -> None:
-        """Initialize attributes."""
-        
+    def __init__(self, position, frame, julian, offset=0) -> None:
+
+        super().__init__(julian, offset)
+
+        time, position = np.array(julian), np.array(position)
+
         if frame not in ["gcrs", "geo", "itrs"]:
             raise ValueError(f"{frame} is not a valid frame.")
-        
-        if position.shape[0] != len(time):
-            raise ValueError(f"position and time data lengths are mismatched being {position.shape[0]} and {len(time)}")
 
-        self.time = time
+        if position.ndim != 2:
+            raise ValueError("position input must be 2-D.")
+
+        if time.ndim != 1:
+            raise ValueError("time input must be 1-D.")
+
+        if position.shape[0] != time.size:
+            raise ValueError(f"position and time data lengths are mismatched: "
+                             f"lengths are {position.shape[0]} and {len(time)}")
+
+        kind_dict = {1: "zero", 2: "linear", 3: "quadratic"}
+        if position.shape[0] > 3:
+            self._kind = "cubic"
+        else:
+            self._kind = kind_dict[position.shape[0]]
 
-        self._GEO = None
         self._GCRS = None
         self._ITRS = None
+        self._length = None
 
-        self.length = None
         self._set_base_position(position, frame)
-    
+
     def __len__(self) -> int:
-        """Return length of position and time data."""
 
-        return self.length
+        return self._length
 
-    def _set_base_position(self, position: np.ndarray, frame:
-                           Literal["gcrs", "geo", "itrs"]) -> None:
+    def _set_base_position(self, position, frame) -> None:
         """Initialize base position.
 
-        This method takes an input position to initialize the object's base
-        position.
-
         Parameters
         ----------
-        position : np.ndarray
-            Array of shape (n, 2) or (n, 3) containing the inputted position
-            data.
+        position : array_like
+            2-D array containing position coordinates.
+
+            Supported coordinate frames include the Geocentric Celestial
+            Reference System (gcrs), International Terrestrial Reference
+            System (itrs), and geographical (geo) system.
         frame : {"gcrs", "geo", "itrs"}
-            String defining the type of input position data as the Geocentric
-            Celestial Reference System (gcrs), International Terrestrial
-            Reference System (itrs), or geographical (geo) data.
+            Frame specifier for the input position.
 
         Notes
         -----
-        The input data must be of shape (n, 3) if `type="itrs"` or `type="gcrs"`
-        where the columns are XYZ cartesian data. The data can be of the shape
-        (n, 2) or (n, 3) if `type="geo"` where the columns are geodetic
-        latitude, terrestrial longitude, and geodetic altitude. When
-        geographical data is entered of shape (n, 2), the height data is
-        assumed to be zero.
+        GCRS and ITRS positions shall have 3 columns representing XYZ
+        cartesian points. GEO positions can have either 2 or 3 columns where
+        the first two represent geodetic latitude and terrestrial longitude;
+        the last column represent geodetic altitude and is assumed zero if not
+        passed in.
         """
-        
-        basePos = position
-        self.length = basePos.shape[0]
-
-        if frame == "geo":
-            if basePos.shape[1] == 2:
-                height = np.zeros((self.length, 1))
-                basePos = np.concatenate((basePos, height), axis=1)
 
-            self._GEO = basePos
+        self._length = position.shape[0]
 
-        elif frame == "gcrs":
-            self._GCRS = basePos
+        if position.shape[1] == 2:
+            height = np.zeros((self._length, 1))
+            position = np.concatenate((position, height), axis=1)
 
+        if frame == "gcrs":
+            x, y, z = [Stroke(self._julian, position[:, i], self._kind) for i in range(3)]
+            self._GCRS = np.array([x, y, z])
         elif frame == "itrs":
-            self._ITRS = basePos
+            x, y, z = [Stroke(self._julian, position[:, i], self._kind) for i in range(3)]
+            self._ITRS = np.array([x, y, z])
+        else:
+            itrs = self._geo_to_itrs(position)
+            x, y, z = [Stroke(self._julian, itrs[:, i], self._kind) for i in range(3)]
+            self._ITRS = np.array([x, y, z])
 
-    def _geo_to_itrs(self, position: np.ndarray) -> np.ndarray:
-        """Convert geographical to itrs coordinates.
+    def _geo_to_itrs(self, position) -> np.ndarray:
+        """Geographical to itrs transformation.
 
         Parameters
         ----------
         position : np.ndarray
-            Array of shape (n, 3) containing geographical coordinates of a
-            position with columns of geodetic latitude, terrestrial longitude,
+            2-D array with columns of geodetic latitude, terrestrial longitude,
             and geodetic altitude given in decimal degrees and kilometres.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n, 3) containing the XYZ itrs position data.
+            2-D array with columns of x, y, z cartesian coordinates in ITRS.
 
         See Also
         --------
-        _itrs_to_geo : Convert itrs to geographical coordinates.
+        _itrs_to_geo : Itrs to geographical transformation.
 
         Notes
         -----
-        This method uses an ellipsoid based model of the Earth to convert a
-        geographical position to itrs cartesian coordinates using the methods
-        described in "Coordinate Systems in Geodesy" by E. J. Krakiwsky and
-        D.E. Wells as presented by Christopher Lum. [KW98a]_ [Lum20]_
+        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.
+        .. [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.
+           Longitude. June 2020.url:https://www.youtube.com/watch?v=4BJ-GpYbZlU.
         """
 
-        a = 6378.137
-        b = 6356.752314245
+        a = 6378.137  # WGS84 major axis in km.
+        b = 6356.752314245  # WGS84 minor axis in km.
 
-        lat, lon = np.radians(position[:, 0]), np.radians(position[:, 1])
+        lat, lon = np.copy(np.radians(position[:, 0:2])).T
+        clat, slat = np.cos(lat), np.sin(lat)
+        clon, slon = np.cos(lon), np.sin(lon)
 
         e = np.sqrt(1 - b ** 2 / a ** 2)
-        N = a / np.sqrt(1 - e ** 2 * np.sin(lat) ** 2)
+        n = a / np.sqrt(1 - e ** 2 * slat ** 2)
         h = position[:, 2]
 
-        x = ((N + h) * np.cos(lat) * np.cos(lon)).reshape((-1, 1))
-        y = ((N + h) * np.cos(lat) * np.sin(lon)).reshape((-1, 1))
-        z = ((N * (1 - e ** 2) + h) * np.sin(lat)).reshape((-1, 1))
-        itrs = np.concatenate((x, y, z), axis=1)
+        x = ((n + h) * clat * clon).reshape((-1, 1))
+        y = ((n + h) * clat * slon).reshape((-1, 1))
+        z = ((n * (1 - e ** 2) + h) * slat).reshape((-1, 1))
 
-        return itrs
+        return np.concatenate((x, y, z), axis=1)
 
-    def _itrs_to_geo(self, position: np.ndarray) -> np.ndarray:
-        """Convert itrs to geographical coordinates.
+    def _itrs_to_geo(self, position) -> np.ndarray:
+        """Itrs to geographical transformation.
 
         Parameters
         ----------
         position : np.ndarray
-            Array of shape (n, 3) with rows containing XYZ itrs position data.
+            2-D array with columns of x, y, z cartesian coordinates in ITRS.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n, 3) with columns of geodetic latitude,
-            terrestrial longitude, and geodetic altitude data in degrees and
-            kilometres.
+            2-D array with columns of latitude, longitude, and altitude given
+            in decimal degrees and kilometres.
 
         See Also
         --------
-        _GEO_to_itrs : Convert geographical to itrs coordinates.
-
-        Notes
-        -----
-        Let :math:`x`, :math:`y`, and :math:`z` be the itrs vector components.
-        We can then calculate the latitude, :math:`\phi`, using the following
-        equation:
-
-        .. math:: \phi = 90^\circ - \cos^{-1}\left(\frac{z}{R_\bigoplus}\right)
-
-        where :math:`R_\bigoplus` is the geocentric radius of the Earth. We can
-        also calculate the longitude, :math:`\lambda` using the following:
-
-        .. math::\lambda = \tan^{-1}_2\left(y, x\right)
+        _geo_to_itrs : Geographical to itrs transformation.
         """
 
-        # Cartesian coordinates.
-        x = position[:, 0]
-        y = position[:, 1]
-        z = position[:, 2]
-
-        # Convert to spherical coordinates.
-        radius = np.sqrt(x ** 2 + y ** 2 + z ** 2)
-        theta = np.degrees(np.arccos(z / radius))
-        phi = np.degrees(np.arctan2(y, x))
-
-        # Get Geodetic altitude.
-        alt = self.altitude()
+        # Calculate latitude and longitude.
+        x, y, z = position
+        lat = np.degrees(np.arctan(z / np.sqrt(x ** 2 + y ** 2)))
+        lon = np.degrees(np.arctan2(y, x))
 
-        # Fomulate output array.
-        lat = (90 - theta).reshape(-1, 1)
-        lon = phi.reshape(-1, 1)
-        alt = alt.reshape(-1, 1)
-        geo = np.concatenate((lat, lon, alt), axis=1)
+        return lat, lon
 
-        return geo
+    def geo(self, iso=False, stroke=False) -> np.ndarray:
+        """Return geographical coordinates.
 
-    def geo(self, iso: bool=False) -> np.ndarray:
-        """Return geographical position data.
+        The method can not return the data as both a stroke and iso
+        formatted strings. If both `iso=True` and `stroke=True`, the method
+        will return a tuple of stroke objects.
 
         Parameters
         ----------
         iso : bool, optional
-            Formats the output as ISO6709 sexagesimal position strings if true.
+            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
         -------
-        np.ndarray
-            Array of shape (n, 3) with columns of geodetic latitude,
-            terrestrial longitude, and geodetic altitude data. If
-            `iso=True`, the output is an array of shape (n,) containing
-            standard representation position strings.
+        tuple
+            Tuple of latitude, longitude, and altitude given in decimal
+            degrees.
+
+            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 position data.
-        gcrs : Return gcrs position data.
+        itrs : Return itrs coordinates.
+        gcrs : Return gcrs coordinates.
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
-        >>> coor = Coordinate(position=position, frame="itrs", time=time)
-        >>> coor.geo()
+        Initialize object and get geographical data:
+
+        >>> julian=[2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
+        >>> c = Coordinate(position=position, frame="itrs", julian=julian)
+        >>> c.geo()
         np.array([[-9.38870528e-02, -2.26014826e+01, 5.04126976e+02]])
 
-        We can generate user-friendly position strings by setting `iso=True`.
+        Return formatted output strings by setting `iso=True`:
 
-        >>> coor.geo(iso=True)
+        >>> c.geo(iso=True)
         np.array(['00°05′37.99″S 22°36′05.34″W 504.12697'])
         """
 
-        if self._GEO is None:
-            if self._ITRS is not None:
-                self._GEO = self._itrs_to_geo(position=self._ITRS)
-            else:
-                self._ITRS = self._gcrs_and_itrs(position=self._GCRS, frame="gcrs")
-                self._GEO = self._itrs_to_geo(position=self._ITRS)
+        if self._ITRS is None:
+            x, y, z = self._GCRS
+            x, y, z = [i(self._julian).reshape((-1, 1)) for i in [x, y, z]]
+            gcrs = np.concatenate((x, y, z), axis=1)
+
+            itrs = self._gcrs_and_itrs(gcrs, "gcrs")
+            x, y, z = [Stroke(self._julian, itrs[:, i], self._kind) for i in range(3)]
+            self._ITRS = np.array([x, y, z])
+
+        lat, lon = self._itrs_to_geo(self._ITRS)
+        alt = self.altitude(stroke=True)
+
+        if stroke:
+            return lat, lon, alt
 
-        geo = _ISO6709_representation(position=self._GEO) if iso else self._GEO
+        lat, lon, alt = lat(self._julian), lon(self._julian), alt(self._julian)
 
-        return geo
+        if iso:
+            return _ISO6709_representation(lat, lon, alt)
+        else:
+            return lat, lon, alt
 
-    def era(self) -> np.ndarray:
-        """Return the Earth rotation angles in degrees and decimals.
+    def era(self, stroke=False) -> np.ndarray:
+        """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
         -------
-        np.ndarray
-            Array of shape (n,) containing Earth rotation angles in degrees and
-            decimals.
+        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 denotes the diurnal rotation component of the
-        coordinate misalignment between the ECI and ECEF frames. The angle can
-        be calculated from the following formulation:
+        The Earth rotation angle is calculated using the following:
 
         .. math:: \gamma^\circ = 360.9856123035484\Delta T + 280.46
 
@@ -302,282 +312,365 @@ def era(self) -> np.ndarray:
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
-        >>> coor = Coordinate(position=position, type="itrs", time=time)
-        >>> coor.era()
+        >>> julian = [2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
+        >>> c = Coordinate(position=position, type="itrs", julian=julian)
+        >>> c.era()
         np.array([6.2360075])
         """
 
-        ang = np.zeros((self.length,))
-        jul_data = self.time.julian()
-
         # Multiply time elapsed since J2000 by Earth rotation rate and add
         # J2000 orientation.
-        dJulian = jul_data - 2451545
+        dJulian = self._julian - 2451545
         ang = (360.9856123035484 * dJulian + 280.46) % 360
 
-        return ang
+        if stroke:
+            return Stroke(self._julian, ang, self._kind)
+        else:
+            return ang
 
-    def _gcrs_and_itrs(self, position: np.ndarray, frame: Literal["itrs", "gcrs"]) -> np.ndarray:
-        """Convert between gcrs and itrs positions.
+    def _gcrs_and_itrs(self, position, frame) -> np.ndarray:
+        """Transform between gcrs and itrs coordinates.
 
         Parameters
         ----------
         position : np.ndarray
-            Array of shape (n, 3) representing the input data as XYZ cartesian
-            data.
+            2-D array containing x, y, z cartesian coordinates.
         frame : {"itrs", "gcrs"}
-            The type of input data.
+            Input data frame.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n, 3) representing the output data as XYZ cartesian
-            data.
+            2-D array containing the transformed x, y, z cartesian data.
 
         Notes
         -----
-        This method applies a simplification to converting between gcrs and
-        itrs coordinate by disregarding precession, nutation, and polar motion
-        effects (which will be incorporated in future releases). By
-        disregarding such factors, we can assume the alignment of the gcrs and
-        itrs axes and simplify the conversion to a rotation about the z-axis.
+        The conversion between the gcrs and itrs frames accounts for Earth
+        rotation, precession, and nutation. Polar motion effects are ignored
+        due to their poor predictive nature.
+        """
 
-        We can convert between the gcrs and itrs frames using the following
-        equations:
+        pos_data = np.copy(position)
 
-        .. math:: \vec{V}_{gcrs} = E_3^\gamma\;\vec{V}_{itrs}
+        # Get dependencies.
+        era = self.era()
+        bias_mtrx = bias_matrix()
+        precession_mtrx = precession_matrix(self._julian)
+        nutation_mtrx = nutation_matrix(self._julian)
 
-        .. math:: \vec{V}_{itrs} = E_3^{-\gamma}\;\vec{V}_{gcrs}
+        # Apply nutation-precession model before ERA rotation.
+        if frame == "gcrs":
+            pos_data = np.einsum('ij, kj -> ki', bias_mtrx, pos_data)
+            pos_data = np.einsum('ijk, ik -> ij', precession_mtrx, pos_data)
+            pos_data = np.einsum('ijk, ik -> ij', nutation_mtrx, pos_data)
 
-        where :math:`\gamma` is the Earth rotation angle, :math:`E_3^\theta` is
-        the rotation matrix that rotates a vector around the z-axis by an angle
-        :math:`\theta`, and :math:`\vec{V}` is a position vector in either the
-        ECI or ECEF frames. [Kok17b]_
+        # Pre-process ERA angle and sign.
+        theta = np.radians(era)
+        theta = -theta if frame == "gcrs" else theta
 
-        References
-        ----------
-        .. [Kok17b] Don Koks. Changing Coordinates in the Context of Orbital
-           Mechanics. Cyber and Electronic Warfare Division, Defence Science,
-           and Technology Group, Jan.2017, p. 21.
-        """
+        # Rotate position by the ERA rotational matrix.
+        c, s = np.cos(theta), np.sin(theta)
+        c1, c2, c3 = np.copy(pos_data).T
 
-        theta = -self.era() if frame == "gcrs" else self.era()
-        theta = np.radians(theta)
+        pos_data[:, 0] = c * c1 - s * c2
+        pos_data[:, 1] = s * c1 + c * c2
+        pos_data[:, 2] = c3
 
-        # Construct rotational matrix.
-        A11 = np.cos(theta)
-        A12 = -np.sin(theta)
-        A21 = np.sin(theta)
-        A22 = np.cos(theta)
+        # Apply nutation-precession model after ERA rotation if frame is ITRS.
+        if frame == "itrs":
+            pos_data = np.einsum('ijk, ik -> ij',
+                                 np.linalg.inv(nutation_mtrx), pos_data)
+            pos_data = np.einsum('ijk, ik -> ij',
+                                 np.linalg.inv(precession_mtrx), pos_data)
+            pos_data = np.einsum('ij, kj -> ki',
+                                 np.linalg.inv(bias_mtrx), pos_data)
 
-        # Rotate position data around z-axis by ERA.
-        output = np.zeros((self.length, 3))
-        output[:, 0] = A11 * position[:, 0] + A12 * position[:, 1]
-        output[:, 1] = A21 * position[:, 0] + A22 * position[:, 1]
-        output[:, 2] = position[:, 2]
+        return pos_data
 
-        return output
+    def gcrs(self, stroke=False) -> np.ndarray:
+        """Return gcrs coordinates.
 
-    def gcrs(self) -> np.ndarray:
-        """Return cartesian gcrs position data.
+        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
         -------
-        np.ndarray
-            Array of shape (n,3) with columns of XYZ gcrs position data.
+        tuple
+            Tuple containing x, y, z cartesian coordinates in the gcrs frame.
+
+            The coordinates are stroke objects if `stroke=True`. Otherwise,
+            1-D arrays are returned.
 
         See Also
         --------
-        itrs : Return cartesian itrs position data.
-        geo : Return geographical position data.
+        itrs : Return itrs coordinates.
+        geo : Return geographical coordinates.
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
-        >>> coor = Coordinate(position=position, type="itrs", time=time)
-        >>> coor.gcrs()
+        Calculate gcrs coordinates from itrs coordinates:
+
+        >>> juluian = [2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
+        >>> c = Coordinate(position=position, type="itrs", julian=julian)
+        >>> c.gcrs()
         np.array([[6212.21719598, -2937.10811161, -11.26]])
         """
 
         if self._GCRS is None:
-            if self._ITRS is None:
-                self._ITRS = self._geo_to_itrs(self._GEO)
-            self._GCRS = self._gcrs_and_itrs(self._ITRS, frame="itrs")
+            itrs = np.array([self._ITRS[i](self._julian) for i in range(3)]).T
+            gcrs = self._gcrs_and_itrs(itrs, "itrs")
+            x, y, z = [Stroke(self._julian, gcrs[:, i], self._kind) for i in range(3)]
+            self._GCRS = np.array([x, y, z])
 
-        return self._GCRS
+        x, y, z = self._GCRS[0], self._GCRS[1], self._GCRS[2]
+        if stroke:
+            return x, y, z
+        else:
+            return x(self._julian), y(self._julian), z(self._julian)
 
-    def itrs(self) -> np.ndarray:
-        """Return cartesian itrs position data.
+    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
         -------
-        np.ndarray
-            Array of shape (n,3) with columns of XYZ itrs position data.
+        tuple
+            Tuple containing x, y, z cartesian data in the itrs frame.
+
+            The coordinates are stroke objects if `stroke=True`. Otherwise,
+            1-D arrays are returned.
 
         See Also
         --------
-        gcrs : Return cartesian gcrs position data.
-        geo : Return geographical position data.
+        gcrs : Return gcrs coordinates.
+        geo : Return geographical coordinates.
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
-        >>> coor = Coordinate(position=position, type="gcrs", time=time)
-        >>> coor.itrs()
+        Calculate itrs coordinates from gcrs coordinates:
+
+        >>> julian = [2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
+        >>> c = Coordinate(position=position, type="gcrs", julian=julian)
+        >>> c.itrs()
         np.array([[6461.30569276, -2338.75507354, -11.26]])
         """
 
         if self._ITRS is None:
-            if self._GCRS is not None:
-                self._ITRS = self._gcrs_and_itrs(position=self._GCRS, frame="gcrs")
-            else:
-                self._ITRS = self._geo_to_itrs(position=self._GEO)
+            gcrs = np.array([self._GCRS[i](self._julian) for i in range(3)]).T
+            itrs = self._gcrs_and_itrs(gcrs, "gcrs")
+            x, y, z = [Stroke(self._julian, itrs[:, i], self._kind) for i in range(3)]
+            self._ITRS = np.array([x, y, z])
 
-        return self._ITRS
+        x, y, z = self._ITRS[0], self._ITRS[1], self._ITRS[2]
+        if stroke:
+            return x, y, z
+        else:
+            return x(self._julian), y(self._julian), z(self._julian)
 
-    def _get_ang(self, u: np.ndarray, v: np.ndarray) -> np.ndarray:
-        """Calculate degree angle bewteen two vectors.
+    def _get_ang(self, u, v) -> np.ndarray:
+        """Return the degree angle bewteen two vectors.
 
         Parameters
         ----------
         u, v : np.ndarray
-            Arrays of shape (n,3) with rows of XYZ cartesian data.
+            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
-            Array of shape (n,) containing the degree angle between the two
-            arrays.
+            1-D array containing the degree angle between to vector arrays.
         """
 
-        num = np.sum(u * v, axis=1)
-        denom = np.linalg.norm(u, axis=1) * np.linalg.norm(v, axis=1)
-        ang = np.degrees(np.arccos(num / denom))
+        ua = None if u.ndim == 1 else 1
+        va = None if v.ndim == 1 else 1
+
+        n = np.sum(u * v, axis=(ua or va))
+        d = np.linalg.norm(u, axis=ua) * np.linalg.norm(v, axis=va)
+        ang = np.degrees(np.arccos(n / d))
 
         return ang
 
-    def horizontal(self, location: Any) -> Tuple:
-        """Return horizontal position data in degrees and decimals.
+    def _altitude(self, location) -> Stroke:
+        """Return the altitude angle of a satellite.
 
-        In this method, the azimuth is calculated as increasing degrees
-        clockwise from North and ranges from 0 to 360 degrees. The altitude is
-        calculated as increasing degrees above the observer's local horizon and
-        ranges from 0 to 90 degrees.
+        Parameters
+        ----------
+        location : GroundPosition
+            Eath bound origin of the horizontal system.
+
+        Returns
+        -------
+        Stroke
+            Stroke containing the altitude of the satellite.
+        """
+
+        # Get origin of horizontal system.
+        lat, lon, radius = location.lat, location.lon, location.radius
+        gnd_itrs = self._geo_to_itrs(np.array([[lat, lon, 0]])).reshape((3,))
+
+        return 90 - self._get_ang(self._ITRS - gnd_itrs, gnd_itrs)
+
+    def _azimuth(self, location) -> Stroke:
+        """Return the azimuth angle of the satellite.
 
         Parameters
         ----------
         location : GroundPosition
-            Ground location defining the center of the horizontal system.
+            Eath bound origin of the horizontal system.
+
+        Returns
+        -------
+        Stroke
+            Stroke containing the azimuth of the satellite.
+        """
+
+        # Get origin of horizontal system.
+        lat, lon, radius = location.lat, location.lon, location.radius
+        gnd_itrs = self._geo_to_itrs(np.array([[lat, lon, 0]])).reshape((3,))
+
+        los = self._ITRS - gnd_itrs
+
+        # Find Earth tangent passing through horizontal origin and z-axis.
+        st = [0, 0, radius / np.sin(np.radians(lat))] - gnd_itrs
+
+        # Find LOS projection on tangent plane.
+        los_on_n = np.sum(los * gnd_itrs)
+        los_proj = los - np.array([los_on_n * gnd_itrs[i] for i in range(3)]) / radius ** 2
+
+        # Determing azimuth.
+        direc = np.cross(st, gnd_itrs)
+        neg_ind = np.sum(los_proj * direc) < 0
+        az = self._get_ang(st, los_proj)
+        az = 360 * neg_ind - 1 * (2 * neg_ind - 1) * az
+
+        return az
+
+    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 altitude 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 of length two containing the altitude and azimuth data as
-            decimals and degrees in NumPy arrays of shape (n,).
+            Tuple containing the altitude and azimuth angles.
+
+            The angles are stroke objects if `stroke=True`. Otherwise, 1-D
+            arrays are returned.
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
+        Calculate the horizontal coordinates at (52.1579, -106.6702):
+
+        >>> julian=[2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
         >>> location = GroundPosition(52.1579, -106.6702)
-        >>> coor = Coordinate(position=position, frame="ecef", time=time)
-        >>> coor.horizontal(location=location)
+        >>> c = Coordinate(position=position, frame="ecef", julian=julian)
+        >>> c.horizontal(location=location)
         (array([-40.8786098]), array([94.73615482]))
         """
 
         if self._ITRS is None:
             self.itrs()
 
-        # Convert observer position into cartesian coordinates.
-        lat, lon, radius = location.lat, location.lon, location.radius
-        GEO_data = self._geo_to_itrs(np.array([[lat, lon, 0]]))
-        obs = np.repeat(GEO_data, self.length, 0)
-
-        # Determine line of sight vector then altitude.
-        LOS = self._ITRS - obs
-        Alt = 90 - self._get_ang(LOS, obs)
+        # Calculate the altitude and azimuth angles.
+        alt = self._altitude(location)
+        az = self._azimuth(location)
 
-        # Find surface tangent vector passing through z-axis.
-        k_hat = np.repeat(np.array([[0, 0, 1]]), self.length, axis=0)
-        beta = np.radians(lat)
-        tangent = (k_hat.T * radius / np.sin(beta)).T - obs
+        if stroke:
+            return alt, az
+        else:
+            return alt(self._julian), az(self._julian)
 
-        # Find LOS projection on tangent plane.
-        norm_proj = (obs.T * np.sum(LOS * obs, axis=1) / radius ** 2).T
-        proj_LOS = LOS - norm_proj
-
-        # Determing azimuth.
-        reference = np.cross(tangent, obs)
-        neg_ind = np.where(np.sum(proj_LOS * reference, axis=1) < 0)[0]
-        Az = self._get_ang(tangent, proj_LOS)
-        Az[neg_ind] = 360 - self._get_ang(tangent[neg_ind], proj_LOS[neg_ind])
-
-        return Alt, Az
-
-    def off_nadir(self, location: Any) -> np.ndarray:
-        """Return the off-nadir angle to a ground location.
+    def off_nadir(self, location, stroke=False) -> np.ndarray:
+        """Return off-nadir angle in decimal degrees.
 
-        This method calculates the off-nadir angle, in degrees and decimals, to
-        the input ground position at each satellite position.
+        The off-nadir angle is the angular distance of a ground location from
+        the satellites nadir.
 
         Parameters
         ----------
         location : GroundPosition
-            Calculate off-nadir angles for the specified ground location.
+            Ground location of interest.
+        stroke : bool, optional
+            Formats output as a stroke objects if true. Otherwise, the output
+            is a 1-D arrays.
 
         Returns
         -------
-        np.ndarray
-            Array of shape (n,) containing off-nadir angle data in degrees and
-            decimals.
+        Union[Stroke, np.ndarray]
+            Off-nadir angles in decimal degrees.
 
-        Notes
-        -----
-        The off-nadir angle is the acute angle measured in increasing degrees
-        from the satellites nadir to the line joining the satellite to the
-        ground location of interest.
+            The off-nadir angles are returned as a stroke object if
+            `stroke=True`. Otherwise, a 1-D array is returned.
+
+        Examples
+        --------
+        Calculate the off-nadir angle for the location (52.1579, -106.6702):
+
+        >>> toronto = GroundPosition(52.1579, -106.6702)
+
+        >>> julian = [30462.70517808 30462.70583972]
+        >>> position = [[-1148.75741119 -5527.65458256  3930.33125383]
+        ...             [-1186.47076862 -5258.03277545  4276.52904538]]
+        >>> c = Coordinate(position=position, frame="itrs", julian=julian, offset=2430000)
+        >>> c.off_nadir(location=toronto)
+        array([67.08711357 65.27748945])
         """
 
         if self._ITRS is None:
             self.itrs()
 
         lat, lon = location.lat, location.lon
-        geo_data = self._geo_to_itrs(np.array([[lat, lon, 0]]))
-        obs = np.repeat(geo_data, self.length, 0)
+        loc = self._geo_to_itrs(np.array([[lat, lon, 0]]))
+        ang = self._get_ang((self._ITRS - loc).reshape((3,)), self._ITRS)
 
-        LOS = np.subtract(self._ITRS, obs)
-        ang = self._get_ang(LOS, self._ITRS)
+        if stroke:
+            return ang
+        else:
+            return ang(self._julian)
 
-        return ang
-
-    def _WGS84_radius(self, lattitude: np.ndarray) -> np.ndarray:
-        """Calculate the Earth's geocentric radius using WGS84.
+    def _WGS84_radius(self, lattitude) -> np.ndarray:
+        """Return Earth's geocentric radius using WGS84.
 
         Parameters
         ----------
         latitude : np.ndarray
-            Array of shape (n,) representing the geocentric latitude of a
-            ground location in degrees and decimals.
+            1-D array containing latitude in decimal degrees.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n,) containing the Earth's geocentric radius in
-            kilometres.
+            1-D array containing the geocentric radius in kilometres.
 
         Notes
         -----
-        By using an Earth ellipsoid with the WGS84 parameters of
-        :math:`a=6378.137` and :math:`b=6356.7523142`, the geocentric radius
-        can be calculated using the following formulation:
-
-        .. math:: r = \sqrt{\frac{(a^2\cos(\beta))^2 + (b^2\sin(\beta))^2}{(a\cos(\beta))^2 + (b\sin(\beta))^2}}
-
-        where :math:`\beta` is the observer's latitude. [Tim18]_
+        The radius is calculated using the WGS84 Earth ellipsoid model as
+        detailed by "Earth Radius by Latitude" by Timur. [Tim18]_
 
         References
         ----------
@@ -592,64 +685,69 @@ def _WGS84_radius(self, lattitude: np.ndarray) -> np.ndarray:
         a = 6378.137
         b = 6356.752314245
 
-        c_phi, s_phi = np.cos(phi), np.sin(phi)
+        c, s = np.cos(phi), np.sin(phi)
 
-        num = (a ** 2 * c_phi) ** 2 + (b ** 2 * s_phi) ** 2
-        denom = (a * c_phi) ** 2 + (b * s_phi) ** 2
+        num = (a ** 2 * c) ** 2 + (b ** 2 * s) ** 2
+        denom = (a * c) ** 2 + (b * s) ** 2
         radius = np.sqrt(num / denom)
 
         return radius
 
-    def altitude(self) -> np.ndarray:
-        """Return the geodetic altitude above Earth's surface in kilometres.
+    def altitude(self, stroke=False) -> 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
         -------
-        np.ndarray
-            Array of shape (n,) containing satellite altitudes in kilometres.
+        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
+        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.
+        .. [KW98b] E. J. Krakiwsky and D. E. Wells. Coordinate Systems in
+           Geodesy. Jan. 1998, pp. 31–33.
 
         Examples
         --------
-        >>> time = Time(julian=np.array([2454545]))
-        >>> position = np.array([[6343.82, -2640.87, -11.26]])
-        >>> coor = Coordinate(position=position, type="itrs", time=time)
-        >>> coor.altitude()
+        >>> julian=[2454545]
+        >>> position = [[6343.82, -2640.87, -11.26]]
+        >>> c = Coordinate(position=position, type="itrs", julian=julian)
+        >>> c.altitude()
         np.array([504.1269764])
         """
 
         if self._ITRS is None:
             self.itrs()
 
-        itrs_data = self._ITRS
-        x = itrs_data[:, 0]
-        y = itrs_data[:, 1]
-        z = itrs_data[:, 2]
+        # Define WGS84 Parameters.
+        a, b = 6378.137, 6356.752314245
 
-        a = 6378.137
-        b = 6356.752314245
+        tol = 10e-10
 
-        epsilon = 10e-10
+        x, y, z = self._ITRS
 
         e = np.sqrt(1 - b ** 2 / a ** 2)
         p = np.sqrt(x ** 2 + y ** 2)
 
-        N = a
         h = np.sqrt(x ** 2 + y ** 2 + z ** 2) - np.sqrt(a * b)
-        phi = np.arctan((z / p) * (1 - (e ** 2 * N) / (N + h)) ** -1)
+        phi = np.arctan((z / p) * (1 - (e ** 2 * a) / (a + h)) ** -1)
 
         def new_vals(a, b, e, p, N, h, phi, z):
 
@@ -659,28 +757,47 @@ def new_vals(a, b, e, p, N, h, phi, z):
 
             return N, h, phi
 
-        Np, hp, phip = new_vals(a, b, e, p, N, h, phi, z)
+        Np, hp, phip = new_vals(a, b, e, p, a, h, phi, z)
 
-        while (np.mean(hp - h) > a * epsilon) and (np.mean(phip - phi) > epsilon):
+        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)
+            n, h, phi = Np, hp, phip
+            Np, hp, phip = new_vals(a, b, e, p, n, h, phi, z)
 
-        return h
+        if stroke:
+            return h
+        else:
+            return h(self._julian)
 
-    def distance(self, location: Any) -> np.ndarray:
-        """Return the distance to a ground location.
+    def distance(self, location, stroke=False) -> np.ndarray:
+        """Return distance to the ground location.
 
         Parameters
         ----------
         location : GroundPosition
-            Calculate distances to the specified ground location.
+        stroke : bool, optional
+            Formats output as a stroke object if true. Otherwise, the output
+            is a 1-D array.
 
         Returns
         -------
-        np.ndarray
-            Array of shape (n,) containing distances in kilometres between the
-            satellite and ground location for each satellite position.
+        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
+        --------
+        Calculate the satellite distances from position (52.1579, -106.6702):
+
+        >>> times = [30462.5, 30462.50069444]
+        >>> position = [[6343.81620221, -2640.87223125, -11.25541802],
+        ...             [6295.64583763, -2718.09271472, 443.08232543]]
+        >>> location = GroundPosition(52.1579, -106.6702)
+        >>> c = Coordinate(position=position, frame="itrs", julian=julian, offset=2430000)
+        >>> c.distance(location=location)
+        np.array([9070.49268746, 8776.7179543])
         """
 
         if self._ITRS is None:
@@ -688,10 +805,9 @@ def distance(self, location: Any) -> np.ndarray:
 
         lat, lon = location.lat, location.lon
         gnd_itrs = self._geo_to_itrs(np.array([[lat, lon, 0]]))
-        gnd_itrs = np.repeat(gnd_itrs, self.length, 0)
-
-        # Find LOS vector norm.
-        LOS = self._ITRS - gnd_itrs
-        distances = np.linalg.norm(LOS, axis=1)
+        distance = np.linalg.norm(self._ITRS - gnd_itrs)
 
-        return distances
+        if stroke:
+            return distance
+        else:
+            return distance(self._julian)
diff --git a/celest/satellite/nutation_precession.py b/celest/satellite/nutation_precession.py
new file mode 100644
index 0000000..2e4db25
--- /dev/null
+++ b/celest/satellite/nutation_precession.py
@@ -0,0 +1,158 @@
+
+
+from celest.satellite._astronomical_quantities import (
+    mean_obliquity, nutation_components, 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 = 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)
+
+    # Construct matrix.
+    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.
+    """
+
+    # Get time constants.
+    t = (julian - 2451545.0) / 36525
+    t2 = t * t
+    t3 = t2 * t
+
+    # Get angles.
+    delta_psi, delta_eps = nutation_components(julian)
+    eps_A = 3600 * mean_obliquity(julian) - 46.84024 * t - 0.00059 * t2 + \
+        0.001813 * t3
+
+    a1 = np.radians(eps_A / 3600)
+    a2 = - np.radians(delta_psi / 3600)
+    a3 = - np.radians(eps_A / 3600 + delta_eps / 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)
+
+    # Construct matrix.
+    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/satellite.py b/celest/satellite/satellite.py
index 4c0554b..dd55042 100644
--- a/celest/satellite/satellite.py
+++ b/celest/satellite/satellite.py
@@ -1,153 +1,93 @@
-"""Satellite orbital representations and coordinate conversions.
 
-This module contains the `Satellite` class to localize satellite-related
-information and functionality.
-"""
 
-
-from celest.core.decorators import set_module
-from celest.core.interpolation import _interpolate
-from typing import Any, Literal, Tuple
+from celest.satellite.coordinate import Coordinate
 import pandas as pd
-import numpy as np
 
 
-@set_module('celest.satellite')
-class Satellite(object):
-    """Localize satellite information and functionality.
+class Satellite(Coordinate):
+    """Satellite(position, frame, julian, offset=0)
 
-    The `Satellite` class represents a satellite, be it artificial or natural,
-    and allows for the position to be represented with time through multiple
-    representations.
+    Satellite abstraction for satellite-ground encounters.
 
-    Parameters
-    ----------
-    position : Coordinate
-        Satellite time-evolving positions.
+    `julian + offset` is the Julian time in the J2000 epoch associated with
+    input positions. If `frame=="geo"`, the `position` input can have 2
+    columns (latitude, longitude) or three columns (latitude, longitude,
+    altitude). If no geodetic altitude is provided, it is assumed zero. Other
+    frames require three columns.
 
-    Attributes
+    Parameters
     ----------
-    time : Time
-        Times associated with the satellite positions.
-    position : Coordinate
-        Position of the satellite.
-
-    Methods
-    -------
-    interpolate(windows, factor)
-        Interpolate satellite data for window times and positions.
-    save_data(fname, delimiter, times, positions)
-        Save satellite time and position data.
-    """
-
-    def __init__(self, position: Any) -> None:
-        """Initialize attributes."""
-
-        self.time = position.time
-        self.position = position
-        self.length = len(position)
-    
-    def __len__(self):
-        """Return length of satellite position and time data."""
-
-        return self.length
+    position : array_like
+        2-D array containing position coordinates.
+
+        Supported coordinate frames include the Geocentric Celestial Reference
+        System (gcrs), International Terrestrial Reference System (itrs), and
+        geographical (geo) system.
+    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.
     
-    def interpolate(self, windows: Any, factor: int=5) -> None:
-        """Interpolate satellite data for window times and positions.
-        
-        This method takes a series of windows and will interpolate the
-        satellites position and time data in the regions corresponding with
-        window times.
-
-        Parameters
-        ----------
-        windows : WindowList
-            Windows defining the regions to interpolate.
-        factor : int, optional
-            The interpolation factor will increase the number of points in
-            window regions by the factor `factor`. A strictly positive value
-            should be used.
-        
-        Notes
-        -----
-        The desire for window defined interpolation stems from the need for
-        more detailed data around encounter opportunities to get more precise
-        start and end times as well as positions for satellite orientations.
-        
-        Examples
-        --------
-        If `IMG_windows` be a series of imaging encounters, we can interpolate
-        within such regions as follows:
-
-        >>> satellite.interpolate(windows=IMG_windows, factor=5)
-        """
-        
-        indices = np.zeros((0,), dtype=int)
-        time = self.time.julian()
-        position = self.position.gcrs()
-
-        for window in windows:
-            start = window.start
-            end = window.end
-
-            temp_ind = np.where((start < time) & (time < end))[0]
-            indices = np.concatenate((indices, temp_ind))
-        
-        indices = np.split(indices, np.where(np.diff(indices) > 1)[0])
-        indices = np.array(indices, dtype=object)
+    Examples
+    --------
+    Initialize `Satellite` using gcrs positions:
+
+    >>> julian = [30462.50, 30462.50]
+    >>> position = [[-4681.50824149 -5030.09119386     0.        ]
+    ...             [-4714.35351825 -4978.74325953   454.41492765]]
+    >>> s = Satellite(position=position, frame="gcrs", julian=julian, offset=0)
+    """
 
-        time_interp = _interpolate(data=time, factor=factor, dt=2, indices=indices)
-        position_interp = _interpolate(data=position, factor=factor, dt=2, indices=indices)
+    def __init__(self, position, frame, julian, offset=0) -> None:
 
-        self.time._julian = time_interp
-        self.time._length = time_interp.size
+        super().__init__(position, frame, julian, offset)
 
-        self.position.time = self.time
-        self.position._GCRS = position_interp
-        self.position._ITRS = None
-        self.position._GEO = None
-        self.position.length = time_interp.size
+    def __len__(self):
 
-        self.length = len(self.position)
+        return self._length
 
-    def save_data(self, fname: str, delimiter: Literal[",", "\\t"], times:
-                  Tuple=("julian",), positions: Tuple=("gcrs",)) -> None:
-        """Save satellite time and position data.
+    def save_data(self, times=("julian",), positions=("gcrs",), path=None,
+                  sep=",", float_format=None) -> None:
+        """Save satellite data.
 
         Parameters
         ----------
-        fname : str
-            File name of the output file as either a .txt or .csv file.
-        delimiter : str
-            String of length 1 representing the field delimiter for the output
-            file.
         times : Tuple, optional
-            List containing the types of time data to save. Possible list
-            values include "julian", "ut1", "gmst", and "gast".
-        positions : Tuple, optional
-            List containing the types of position data to save. Possible
-            list values include "gcrs", "itrs", and "geo".
+            Tuple containing time representations to save.
 
-        Notes
-        -----
-        It is recommended to use a tab delimiter for .txt files and comma
-        delimiters for .csv files. The method will return an error if the
-        fileName already exists in the current working directory.
+            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_data(fname="data.csv", delimiter=",", positions=positions)
+        >>> finch.save_data(fname="data.csv", positions=positions, delimiter=",")
         """
 
         key_mapping = {
-            "julian": (self.time.julian, "Julian.J2000"),
-            "ut1": (self.time.ut1, "UT1.hr"),
-            "gmst": (self.time.gmst, "GMST.hr"),
-            "gast": (self.time.gast, "GAST.hr"),
-            "itrs": (self.position.itrs, "ITRS.x.km", "ITRS.y.km", "ITRS.z.km"),
-            "gcrs": (self.position.gcrs, "GCRS.x.km", "GCRS.y.km", "GCRS.z.km"),
-            "geo": (self.position.geo, "GEO.lat.deg", "GEO.lon.deg", "GEO.height.km")
+            "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"),
+            "gcrs": (self.gcrs, "GCRS.x.km", "GCRS.y.km", "GCRS.z.km"),
+            "geo": (self.geo, "GEO.lat.deg", "GEO.lon.deg", "GEO.height.km")
         }
 
         data = {}
@@ -171,7 +111,7 @@ def save_data(self, fname: str, delimiter: Literal[",", "\\t"], times:
             columns = position_information[1:]
 
             for i, column in enumerate(columns):
-                data[column] = position_data[:, i]
+                data[column] = position_data[i]
 
         df = pd.DataFrame(data)
-        df.to_csv(fname, sep=delimiter)
+        df.to_csv(path_or_buf=path, sep=sep, float_format=float_format)
diff --git a/celest/satellite/time.py b/celest/satellite/time.py
index 0257a88..c8184f7 100644
--- a/celest/satellite/time.py
+++ b/celest/satellite/time.py
@@ -1,11 +1,5 @@
-"""Time representations.
 
-The time module contains the `Time` class to allow users to input an array of
-times in Julian days and convert to various time representations.
-"""
 
-
-from celest.core.decorators import set_module
 from celest.satellite._astronomical_quantities import (
     equation_of_time, from_julian, sun_right_ascension, equation_of_equinoxes
 )
@@ -13,80 +7,82 @@
 import numpy as np
 
 
-@set_module('celest.satellite')
-class Time(object):
-    """Time representations.
+class Time:
+    """Time(julian, offset=0)
+
+    Time transformations.
 
-    The `Time` class allows a user to convert an array of Julian day times in
-    the J2000 epoch into various time representations.
+    `julian + offset` is the Julian time in J2000 epoch which can be converted
+    into different time representations useful for astronomical calculations.
 
     Parameters
     ----------
-    julian : np.ndarray
-        Array of shape (n,) containing time data in Julian days.
+    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.
-
-    Methods
-    -------
-    true_solar_time(longitude)
-        Return the true solar time (TTs) in hours and decimals.
-    mean_solar_time(longitude)
-        Return the mean solar time (MTs) in hours and decimals.
-    true_hour_angle(longitude)
-        Return the true hour angle in hours and decimals.
-    mean_hour_angle(longitude)
-        Return the mean hour angle in hours and decimals.
-    ut1()
-        Return the universal time (same as GMT) in hours and decimals.
-    julian()
-        Return Julian times.
-    datetime()
-        Return `datetime.datetime` object array.
-    gmst(longitude)
-        Return Greenwich Mean Sidereal Time in hours and decimals.
-    lmst(longitude)
-        Return Local Mean Sidereal Time in hours and decimals.
-    gast(longitude)
-        Return Greenwich Apparent Sidereal Time in hours and decimals.
-    last(longitude)
-        Return Local Apparent Sidereal Time in hours and degrees.
+        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: np.ndarray, offset: float=0) -> None:
-        """Initialize attributes."""
+    def __init__(self, julian, offset=0) -> None:
 
-        self._julian = julian + offset
+        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 data size."""
 
         return self._length
 
-    def true_solar_time(self, longitude: np.ndarray) -> np.ndarray:
-        """Return the true solar time (TTs) in hours and decimals.
+    def true_solar_time(self, longitude) -> np.ndarray:
+        """Return true solar time in decimal hours.
 
         Parameters
         ----------
-        longitude : np.ndarray
-            Array of shape (n,) containing longitude in decimal degrees.
+        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
-            Array of shape (n,) containing true solar time in hours and
-            decimals.
+            1-D array containing true solar time in decimal hours.
 
         See Also
         --------
         mean_solar_time :
-            Return the mean solar time (MTs) in hours and decimals.
+            Return mean solar time in decimal hours.
 
         Notes
         -----
-        The true solar time, :math:`TT_s` can be calculated by first converting
-        Julian data, :math:`JD`, into UTC data by the following:
+        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
 
@@ -104,47 +100,45 @@ def true_solar_time(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2456293.54167])
-        >>> longitude = np.array([147.46])
-        >>> Time(julian=julian).true_solar_time(longitude=longitude)
+        Calculate the true solar time for a constant longitude:
+
+        >>> Time(julian=[2456293.54167]).true_solar_time(longitude=147.46)
         np.array([10.773109])
         """
 
-        julian = self._julian
-
-        EoT = equation_of_time(julian=julian)
+        EoT = equation_of_time(julian=self._julian)
 
-        alpha = np.full((julian.size,), -12)
-        alpha[np.where(julian % 1 < 0.5)] = 12
+        idx = np.where(self._julian % 1 < 0.5)
 
-        UTC = 24 * (julian % 1) + alpha
-        TST = (UTC + (longitude + EoT) / 15) % 24
+        UTC = 24 * (self._julian % 1) - 12
+        UTC[idx] = UTC[idx] + 12
+        TST = (UTC + (np.array(longitude) + EoT) / 15) % 24
 
         return TST
 
-    def mean_solar_time(self, longitude: np.ndarray) -> np.ndarray:
-        """Return the mean solar time (MTs) in hours and decimals.
+    def mean_solar_time(self, longitude) -> np.ndarray:
+        """Return mean solar time in decimal hours.
 
         Parameters
         ----------
-        longitude : np.ndarray
-            Array of shape (n,) containing longitude in degrees and decimals
+        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
-            Array of shape (n,) containing mean solar time in hours and
-            decimals.
+            1-D array containing mean solar time in decimal hours.
 
         See Also
         --------
         true_solar_time :
-            Return the true solar time (TTs) in hours and decimals.
+            Return true solar time in decimal hours.
 
         Notes
         -----
-        The mean solar time, :math:`MT_s` can be calculated by first converting
-        Julian data, :math:`JD`, into UTC data by the following:
+        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
 
@@ -156,45 +150,43 @@ def mean_solar_time(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2456293.54167])
-        >>> longitude = np.array([147.46])
-        >>> Time(julian=julian).mean_solar_time(longitude=longitude)
+        Calculate the mean solar time for a constant longitude:
+
+        >>> Time(julian=2456293.54167).mean_solar_time(longitude=147.46)
         np.array([10.830667])
         """
 
-        julian = self._julian
-
-        alpha = np.full((julian.size,), -12)
-        alpha[np.where(julian % 1 < 0.5)] = 12
+        idx = np.where(self._julian % 1 < 0.5)
 
-        UTC = 24 * (julian % 1) + alpha
-        MST = (UTC + longitude / 15) % 24
+        UTC = 24 * (self._julian % 1) - 12
+        UTC[idx] = UTC[idx] + 12
+        MST = (UTC + np.array(longitude) / 15) % 24
 
         return MST
 
-    def true_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
-        """Return the true hour angle in hours and decimals.
+    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
-        the given times and its position at local solar noon. The value falls
-        within the range :math:`[0, 360)` and is measured in increasing degrees
-        past local solar noon.
+        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 : np.ndarray
-            Array of shape (n,) containing longitude in degrees and decimals.
+        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
-            Array of shape (n,) containing true solar hour angles in hours
-            and decimals.
+            1-D array containing true solar hour angles in decimal hours.
 
         See Also
         --------
         mean_hour_angle :
-            Return the mean hour angle in hours and decimals.
+            Return mean hour angle in decimal hours.
 
         Notes
         -----
@@ -213,8 +205,10 @@ def true_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85])
-        >>> lon = np.array([-105, -118.24])
+        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])
@@ -225,29 +219,29 @@ def true_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
 
         return HRA
 
-    def mean_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
-        """Return the mean hour angle in hours and decimals.
+    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
-        the given times and its position at local solar noon. The value falls
-        within the range :math:`[0, 360)` and is measured in increasing degrees
-        past local solar noon.
+        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 : np.ndarray
-            Array of shape (n,) containing longitude in degrees and decimals.
+        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
-            Array of shape (n,) containing mean solar hour angles in hours and
-            decimals.
+            1-D containing mean solar hour angles in decimal hours.
 
         See Also
         --------
         true_hour_angle :
-            Return the true hour angle in hours and decimals.
+            Return true hour angle in decimal hours.
 
         Notes
         -----
@@ -266,9 +260,9 @@ def mean_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2456293.5416666665])
-        >>> longitude = np.array([147.46])
-        >>> Time(julian=julian).mean_hour_angle(longitude=longitude)
+        Calculate the meann hour angle for a single longitude:
+
+        >>> Time(julian=2456293.5416666665).mean_hour_angle(longitude=147.46)
         np.array([22.83066666])
         """
 
@@ -278,18 +272,16 @@ def mean_hour_angle(self, longitude: np.ndarray) -> np.ndarray:
         return HRA
 
     def ut1(self) -> np.ndarray:
-        """Return the universal time (same as GMT) in hours and decimals.
+        """Return the universal time (same as GMT) in decimal hours.
 
-        This method returns a universal time by calculating the mean solar time
-        at the Greenwich prime meridian. 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.
+        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
-            Array of shape (n,) containing universal time in hours and
-            decimals.
+            1-D array containing universal time in decimal hours.
 
         Notes
         -----
@@ -304,7 +296,7 @@ def ut1(self) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).ut1()
         np.array([6.00000,
                   8.40000,
@@ -316,17 +308,16 @@ def ut1(self) -> np.ndarray:
         return ut1
 
     def julian(self) -> np.ndarray:
-        """Return Julian times.
+        """Convenience method to return Julian times.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n,) containing the Julian times in days and
-            decimals.
+            1-D array containing the Julian times in decimal days.
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).julian()
         np.array([2455368.75,
                   2459450.85,
@@ -341,11 +332,11 @@ def datetime(self) -> np.ndarray:
         Returns
         -------
         np.ndarray
-            Array of shape (n,) containing `datetime.datetime` objects.
+            1-D array containing `datetime.datetime` objects.
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> 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)
@@ -367,26 +358,18 @@ def datetime(self) -> np.ndarray:
         return datetime_arr
 
     def gmst(self) -> np.ndarray:
-        """Return Greenwich Mean Sidereal Time in hours and decimals.
+        """Return Greenwich Mean Sidereal Time in decimal hours.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n,) containing Greenwich Mean Sidereal Time in
-            hours and decimals.
+            1-D array containing Greenwich Mean Sidereal Time in hours and
+            decimals.
 
         Notes
         -----
-        The Greenwich Mean Sidereal Time can be calculated from the following
-        equation:
-
-        .. math:: GMST = \left(280.46 + 360.99(j - 2451545) + 0.000387933T^2 - \\frac{T^3}{38710000}\\right) % 360 / 15
-
-        where
-
-        .. math:: T = \\frac{JD - 2451545}{36525}
-
-        and :math:`JD` is the current Julian day. [Mee98a]_
+        The Greenwich Mean Sidereal Time is calculated using the methods
+        described in Astronomical Algorithms by Jean Meeus. [Mee98a]_
 
         References
         ----------
@@ -395,40 +378,38 @@ def gmst(self) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).gmst()
         np.array([23.955316,
                   6.589391,
                   7.723214])
         """
 
-        julData = self._julian
-
-        T = (julData - 2451545) / 36525
+        T = (self._julian - 2451545) / 36525
         T2 = T * T
         T3 = T2 * T
 
         gmst = 280.46061837
-        gmst = gmst + 360.98564736629 * (julData - 2451545)
+        gmst = gmst + 360.98564736629 * (self._julian - 2451545)
         gmst = gmst + 0.000387933 * T2
         gmst = gmst - T3 / 38710000
         gmst = gmst % 360 / 15
 
         return gmst
 
-    def lmst(self, longitude: np.ndarray) -> np.ndarray:
-        """Return Local Mean Sidereal Time in hours and decimals.
+    def lmst(self, longitude) -> np.ndarray:
+        """Return Local Mean Sidereal Time in decimal hours.
 
         Parameters
         ----------
-        longitude : np.ndarray
-            Array of shape (n,) containing longitude in degrees and decimals.
+        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
-            Array of shape (n,) containing Local Mean Sidereal Time in hours and
-            decimals.
+            1-D array containing Local Mean Sidereal Time in decimal hours.
 
         Notes
         -----
@@ -448,7 +429,7 @@ def lmst(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).lmst(longitude=150)
         np.array([9.98419514,
                   16.62892539,
@@ -463,13 +444,13 @@ def lmst(self, longitude: np.ndarray) -> np.ndarray:
         return lmst
 
     def gast(self) -> np.ndarray:
-        """Return Greenwich Apparent Sidereal Time in hours and decimals.
+        """Return Greenwich Apparent Sidereal Time in decimal hours.
 
         Returns
         -------
         np.ndarray
-            Array of shape (n,) containing Greenwich Apparent Sidereal Time in
-            hours and degrees.
+            1-D array containing Greenwich Apparent Sidereal Time in hours and
+            degrees.
 
         Notes
         -----
@@ -479,11 +460,11 @@ def gast(self) -> np.ndarray:
 
         .. math:: GAST^h = GMST^h + EoT^h
 
-        where :math:`EoT` is the equation of time in hours and decimals.
+        where :math:`EoT` is the equation of time in decimal hours.
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).gast()
         np.array([23.955596,
                   6.589146,
@@ -497,31 +478,26 @@ def gast(self) -> np.ndarray:
 
         return gast
 
-    def last(self, longitude: np.ndarray) -> np.ndarray:
+    def last(self, longitude) -> np.ndarray:
         """Return Local Apparent Sidereal Time in hours and degrees.
 
         Parameters
         ----------
-        longitude : np.ndarray
-            Array of shape (n,) containing the actual astronomical longitude
-            in degrees and decimals.
+        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
-            Array of shape (n,) containing Local Apparent Sidereal Time in
-            hours and degrees.
+            1-D array containing Local Apparent Sidereal Time in hours and
+            degrees.
 
         Notes
         -----
-        The Local Apparent Sidereal Time can be calculated using the following
-        formulation:
-
-        .. math:: LAST^h = UT1^h - 12^h + \\alpha^h_{mSun} + EoE^h + \Lambda^\circ/15
-
-        where :math:`\\alpha^h_{mSun}` is the right ascension of the mean Sun
-        position in hours, :math:`EoE^h` is the equation of equinoxes in hours,
-        and :math:`\Lambda` is the observer's longitude in  degrees. [SL13f]_
+        The Local Apparent Sidereal Time is calculated using the methods
+        described in Space-Time Reference Systems by M. Soffel and R.
+        Langhans. [SL13f]_
 
         References
         ----------
@@ -530,7 +506,7 @@ def last(self, longitude: np.ndarray) -> np.ndarray:
 
         Examples
         --------
-        >>> julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        >>> julian = [2455368.75, 2459450.85, 2456293.5416666665]
         >>> Time(julian=julian).last(longitude=150)
         np.array([9.98447567,
                   16.6286799,
@@ -541,6 +517,6 @@ def last(self, longitude: np.ndarray) -> np.ndarray:
         alpha_sun = sun_right_ascension(self._julian) / 15
         EoE = equation_of_equinoxes(self._julian) / 3600
 
-        last = (ut1 - 12 + alpha_sun + EoE + longitude / 15) % 24
+        last = (ut1 - 12 + alpha_sun + EoE + np.array(longitude) / 15) % 24
 
         return last
diff --git a/docs/source/guide/quickstart.rst b/docs/source/guide/quickstart.rst
index f48668a..f0f6fe7 100644
--- a/docs/source/guide/quickstart.rst
+++ b/docs/source/guide/quickstart.rst
@@ -14,14 +14,14 @@ Getting familiar with the library
 ---------------------------------
 
 Celest is broken primarily into two modules: (1) :mod:`celest.satellite` module
-contains all functions related to satellite dynamics such as
+contains all functions related to satellite placement and dynamics such as
 coordinate and time conversions. (2) :mod:`celest.encounter` module contains
-all functions related to the encounters between a satellite and ground location.
+all functions related to encounter calculations between a satellite and ground
+location.
 
 A typical workflow for using Celest is to begin by importing satellite time and
-position data as NumPy arrays which can then be formed into :class:`Time` and
-:class:`Coordinate` objects which are then used to initialize a
-:class:`Satellite` object. The :class:`Satellite` object can then be coupled
+position data as NumPy arrays which can then be used to initialize the
+:class:`Satellite` class. The :class:`Satellite` object can then be coupled
 with encounter information such as a ground location and encounter parameters
 to determine encounter opportunities.
 
@@ -36,16 +36,14 @@ to determine encounter opportunities.
    gcrs = np.loadtxt('gcrs.txt')
 
    # Initialize satellite representation.
-   time = Time(julian, offset=0)
-   position = Coordinate(gcrs, frame='gcrs', time=time)
-   satellite = Satellite(position=position)
+   satellite = Satellite(position=gcrs, frame='gcrs', time=julian, offset=0)
 
    # Define ground position.
    toronto = GroundPosition(latitude=43.65, longitude=-79.38)
 
    # Generate ground location windows.
-   toronto_IMG_windows = windows.generate(satellite=satellite, location=toronto, enc="image", ang=30)
-   toronto_GL_windows = windows.generate(satellite=satellite, location=toronto, enc="data link", ang=10)
+   toronto_IMG_windows = windows.generate(satellite=satellite, location=toronto, enc="image", ang=30, lighting=1)
+   toronto_GL_windows = windows.generate(satellite=satellite, location=toronto, enc="data link", ang=10, lighting=0)
 
    # Save satellite encounter windows.
    toronto_IMG_windows.save(fname="toronto_IMG_windows.csv", delimiter=",")
diff --git a/docs/source/guide/tutorials.rst b/docs/source/guide/tutorials.rst
index 276030e..d637bec 100644
--- a/docs/source/guide/tutorials.rst
+++ b/docs/source/guide/tutorials.rst
@@ -4,8 +4,8 @@ Tutorials
 Position and Time Conversions
 -----------------------------
 
-Celest allows for a variety of time and position representations that may
-prove useful in a mission planning context. Here, we will endeavour to explore
+Celest allows for various time and position representations that may
+prove useful in a mission planning context. Here, we will explore
 the basic interactions with the :class:`Time` and :class:`Coordinate` classes.
 For details on all possible conversions, refer to the :ref:`Time <Time Class>`
 and :ref:`Coordinate <Coordinate Class>` documentation.
@@ -26,9 +26,9 @@ loaded such that all array values are of a type `float`.
    julian = np.loadtxt('julian.txt')
    gcrs = np.loadtxt('gcrs.txt')
 
-We can then go ahead instantiating our :class:`Time` object. The input data
-should be Julian time data in the J2000 epoch. To incorporate data from other
-epochs, the appropriate offset can be applied using the `offset` keyword. This
+We can then instantiate the :class:`Time` class. The input data
+should be Julian time data in the J2000 epoch. The appropriate offset can be
+applied using the `offset` keyword to incorporate data from other epochs. This
 offset will be added to the input data when stored in the :class:`Time` object.
 
 .. code-block::
@@ -39,8 +39,8 @@ offset will be added to the input data when stored in the :class:`Time` object.
    # Initialize Time object using Mean Julian Date data.
    time = Time(julian, offset=2400000.5)
 
-Different time representations can then be accessed using various methods. Some
-examples are shown below.
+Different time representations can then be accessed through various methods.
+Some examples are shown below.
 
 .. code-block::
 
@@ -54,7 +54,7 @@ examples are shown below.
    datetime_data = time.datetime()
 
 The :class:`Coordinate` class can be instantiated using a data array of
-positions, a frame specifier, and a :class:`Time` object associated with the
+positions, a frame specifier, and time information associated with the
 position data. Three input position frames are currently supported by the
 :class:`Coordinate` class: the geocentric celestial reference system (GCRS),
 international terrestrial reference system (ITRS), and geodetic reference
@@ -63,16 +63,19 @@ system (GRS).
 .. code-block::
 
    # Initialize Coordinate object using GCRS data.
-   position = Coordinate(gcrs, frame='gcrs', time=time)
+   position = Coordinate(gcrs, frame='gcrs', time=julian, offset=0)
 
    # Initialize Coordinate object using ITRS data.
-   position = Coordinate(itrs, frame='itrs', time=time)
+   position = Coordinate(itrs, frame='itrs', time=julian, offset=0)
 
    # Initialize Coordinate object using GCRS data.
-   position = Coordinate(geo, frame='geo', time=time)
+   position = Coordinate(geo, frame='geo', time=julian, offset=0)
+
+The :class:`Coordinate` class inherits the :class:`Time` class; as a result,
+any time methods can be accessed using the :class:`Coordinate` class.
 
 Similar to the :class:`Time` class, different coordinate representations can be
-accessed using various methods. Some examples are shown below.
+accessed through various methods. Some examples are shown below.
 
 .. code-block::
 
@@ -88,14 +91,14 @@ accessed using various methods. Some examples are shown below.
    # Get GRS data with ISO6709 formatted output strings.
    geo_position = position.geo(iso=True)
 
-Notice that some methods require the use of a :class:`GroundPosition` object as
+Notice that some methods require using a :class:`GroundPosition` object as
 a parameter to specify a ground location. The :class:`GroundPosition` object
 can be imported from the encounter module.
 
 Window Generation Workflow
 --------------------------
 
-The basic window generation workflow can be broken down into three stages:
+The primary window generation workflow can be broken down into three stages:
 
 #. Import and prepare time and position data,
 #. Specify ground locations, and
@@ -103,7 +106,7 @@ The basic window generation workflow can be broken down into three stages:
 
 The first step is to import and prepare the time and position data. This
 includes setting up the :class:`Satellite` object that holds the necessary but
-not sufficient information to generate the desired windows.
+insufficient information to generate the desired windows.
 
 .. code-block::
 
@@ -116,11 +119,9 @@ not sufficient information to generate the desired windows.
    gcrs = np.loadtxt('gcrs.txt')
 
    # Initialize satellite representation.
-   time = Time(julian, offset=0)
-   position = Coordinate(gcrs, frame='gcrs', time=time)
-   satellite = Satellite(position=position)
+   satellite = Satellite(position=gcrs, frame='gcrs', time=julian, offset=0)
 
-Next, we specify the ground locations that we wish to generate windows for. To
+Next, we specify the ground locations for which we wish to generate windows. To
 accomplish this, we define a :class:`GroundPosition` object for each location
 we wish to encounter. If various encounter types for one location are desired,
 only one :class:`GroundPosition` object is required.
@@ -133,14 +134,14 @@ only one :class:`GroundPosition` object is required.
 
 We are now ready to generate windows. The :py:func:`windows.generate` function
 takes a satellite and ground location as an input and will populate a
-:class:`Windows` object with possible encounter opportunities for the encounter
+:class:`Windows` object with visible window times for the encounter
 defined by the `enc` and `ang` keywords.
 
 There are two encounter types that Celest currently supports: (1) imaging
 encounters where the satellite is in view of the ground location, and (2) data
 transmission encounters where the ground location is in view of the satellite.
 The `enc` keyword specifies the type of encounter as either and imaging
-(enc="image") or data transmission (enc="data link") type.
+(`enc="image"`) or data transmission (`enc="data link"`) type.
 
 The `ang` keyword defines the constraint angle that borders a
 viable/non-viable encounter region. The constraint angle type used for imaging
@@ -149,11 +150,14 @@ satellite's nadir to the ground location. Transmission encounters use the
 altitude angle of the satellite as measured in increasing degrees above the
 horizon (as seen from the ground location).
 
+The lighting conditions can also be set. For example, to image only in the
+daylight, we can set `lighting=1`.
+
 .. code-block::
 
    # Generate ground location windows.
-   toronto_IMG_windows = windows.generate(satellite=satellite, location=toronto, enc="image", ang=30)
-   toronto_GL_windows = windows.generate(satellite=satellite, location=toronto, enc="data link", ang=10)
+   toronto_IMG_windows = windows.generate(satellite=satellite, location=toronto, enc="image", ang=30, lighting=1)
+   toronto_GL_windows = windows.generate(satellite=satellite, location=toronto, enc="data link", ang=10, lighting=0)
 
    # Save satellite encounter windows.
    toronto_IMG_windows.save(fname="toronto_IMG_windows.csv", delimiter=",")
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 1dc01fb..dae998b 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -6,25 +6,23 @@
 Celest Documentation
 ====================
 
-Celest is a satellite mission planning program designed for the University of
-Toronto Aerospace Team which is applying the program to a hyperspectral
-imaging CubeSat mission. The library aims to provide the necessary tools to
-plan satellite-ground interactions such as ground target tracking, and
-calculating imaging and ground transmission opportunities.
+Celest is a satellite mission planning software designed for the University of
+Toronto Aerospace Team, which is applying the program to a hyperspectral
+imaging CubeSat mission. The library aims to provide the necessary tools for
+celestial orbital conversions and satellite-to-ground interactions such as
+determining imaging/transmission opportunities.
 
-Being subjected to a design team workflow, Celest is consistently under
-development for improvements in efficiency, functionality, and fidelity. It is 
-then expected that the library will be updated regularly.
+Celest is currently in the development phase and is subject to significant
+changes. The program aims to incorporate a high-order adaptive orbit
+propagator and state-of-the-art scheduling algorithm for agile earth
+observation satellite (AEOS) imaging encounters.
 
 What can the program do?
 ------------------------
 
-The program only requires an input data set for a satellite that contains the
-positions in either a gcrs, itrs, or geographical frame and the
-assosciated Julian times. The program can then be used to calculate various
-position and time representations, satellite to ground encounters, and
-transformations necessary for ground location tracking.
-
+Currently, the program requires satellite position data in either the gcrs,
+itrs, or geographical frame and the corresponding times. Celest can then calculate different position representations, time representations, and ground
+encounter calculations.
 
 What can be expected of the project?
 ------------------------------------
diff --git a/docs/source/modules/encounter.rst b/docs/source/modules/encounter.rst
index 079c80a..9655727 100644
--- a/docs/source/modules/encounter.rst
+++ b/docs/source/modules/encounter.rst
@@ -6,10 +6,9 @@ GroundPosition Class
 
 The :class:`GroundPosition` class is used to represent a ground location that
 can be passed into window generation functionality to produce viable
-ground-satellite encounter opportunities. Currently, the
-:class:`GroundPosition` class is fairly preliminary but will be expanded in the
-future to accommodate more specific and complicated ellipsoid models based on
-the input geographical coordinates.
+ground-satellite encounter opportunities. Currently, the :class:`GroundPosition`
+class is somewhat preliminary but will be expanded in the future to accommodate
+more complicated ellipsoid models based on the input geographical coordinates.
 
 .. autoclass:: celest.encounter.GroundPosition
    :members:
@@ -36,16 +35,45 @@ access to the data. The following details these two classes.
 Windows
 ~~~~~~~
 
-The :class:`Windows` class is interfaced through a series of interactive
-methods. The window data can be interfaced in two ways. The first is through
-the :class:`Windows.windows` attribute which holds the :class:`Window`
-objects in a Pandas :py:func:`Series` data structure indexed by the window
-start time. Window data can also be accessed by iterating over the class.
+The window data held within the :class:`Windows` class can be interfaced in two
+ways: indexing and iterating.
+
+To provide meaningful interactions with the window data, the :class:`Windows`
+is indexed by Julian start times in the J2000 epoch. Although, the start times
+are typically unknown to the user. As a result, a unique indexing scheme was
+adopted by the class to provide data access that might be convenient in
+context. The window associated with the closest start time is returned when
+indexing using an integer or float value. Similarly, when indexing with a tuple
+of values, the unique windows related to the closest start time are returned.
+When two indices are both closest to the start time of a window, the window is
+only returned once. All windows with start times falling within the slice
+(inclusive) are returned if the index is a slice.
+
+This indexing scheme allows users to easily access the data in the
+:class:`Windows` class without prior knowledge of start times. It also
+allows the user to determine the encounter closest to the desired time for
+scheduling purposes.
+
+Window data can also be accessed by iterating over the class.
+
+Examples of object interactions are seen in the following example.
 
 .. code-block::
 
-   # Window accessed by indexing by the start time 2400000.5
-   window_object = windows_object.windows[2400000.5]
+   # Assuming we have a Windows object with three windows with start times
+   # 2405795.5, 2405796.5, and 2405797.5.
+
+   # Window accessed by indexing by the start time 2405796.5.
+   window_object = windows_object[2405796.5]
+   # window_object now contains the window starting at 2405796.5.
+
+   # Window accessed by indexing with a tuple.
+   window_object = windows_object[2405796.5, 2405796.7]
+   # window_object is a list containing windows with the 2405795.5, 2405796.5 start times.
+
+   # Window accessed by indexing with a slice.
+   window_object = windows_object[2405796.5:2405796.7]
+   # window_object is a list containing windows with the 2405795.5, 2405796.5 start times.
 
    # Window accessed by iterating over the class.
    for window_object in windows_object:
diff --git a/docs/source/modules/satellite.rst b/docs/source/modules/satellite.rst
index ad49ecc..23b1578 100644
--- a/docs/source/modules/satellite.rst
+++ b/docs/source/modules/satellite.rst
@@ -5,11 +5,10 @@ Time Class
 ----------
 
 The :class:`Time` class allows for flexible and convenient time specifications.
-It takes in Julian data in the J2000 epoch or another epoch with an appropriate
-offset. This offset will be added to the input data when stored in the
-:class:`Time` object. Example usage of the :class:`Time` class to carry out
-conversions can be seen on the tutorials page
-:ref:`here <Position and Time Conversions>`.
+It takes Julian data in the J2000 epoch or another epoch with an appropriate
+offset (added to the times to achieve the J200 epoch). Example usage
+of the :class:`Time` class to carry out conversions can be seen on the
+tutorials page :ref:`here <Position and Time Conversions>`.
 
 .. autoclass:: celest.satellite.Time
    :members:
@@ -19,11 +18,11 @@ conversions can be seen on the tutorials page
 Coordinate Class
 ----------------
 
-The :class:`Coordinate` class is used to represent a satellite's position and
-can be used to carry out coordinate transformations between different
-frames. The Example usage of the :class:`Coordinate` class to carry out
-conversions can be seen on the tutorials page
-:ref:`here <Position and Time Conversions>`.
+The :class:`Coordinate` class, which inherits the :class:`Time` class, is used
+to represent a satellite's position and can be used to carry out coordinate
+transformations between different frames. The Example usage of the
+:class:`Coordinate` class to carry out conversions can be seen in the tutorials
+page :ref:`here <Position and Time Conversions>`.
 
 .. autoclass:: celest.satellite.Coordinate
    :members:
@@ -33,11 +32,11 @@ conversions can be seen on the tutorials page
 Satellite Class
 ---------------
 
-The :class:`Satellite` class is used to hold all necessary information for a
-given satellite such that it can be passed into the :py:func:window.generate`
-function for mission planning. Example usage of the :class:`Satellite` in
-a mission planning capacity can be seen on the tutorials page
-:ref:`here <Window Generation Workflow>`.
+The :class:`Satellite` class, which inherits the :class:`Coordinate` class, is
+used to hold all necessary information for a given satellite such that it can
+be passed into the :py:func:window.generate` function for mission planning.
+Example usage of the :class:`Satellite` in a mission planning capacity can be
+seen on the tutorials page :ref:`here <Window Generation Workflow>`.
 
 .. autoclass:: celest.satellite.Satellite
    :members:
diff --git a/requirements.txt b/requirements.txt
index 195cca9..a7b545e 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,34 +1,25 @@
-alabaster==0.7.12
-astropy==4.3.1
-Babel==2.9.1
-certifi==2021.5.30
-charset-normalizer==2.0.7
-coverage==6.0.2
-docutils==0.17.1
-idna==3.3
-imagesize==1.2.0
-Jinja2==3.0.2
-jplephem==2.16
+astropy==5.0.4
+attrs==21.2.0
+certifi==2021.10.8
+cffconvert==2.0.0
+docopt==0.6.2
+jplephem==2.17
+jsonschema==3.2.0
 julian==0.14
-MarkupSafe==2.0.1
-nose2==0.10.0
-numpy==1.20.3
-packaging==21.2
-pandas==1.3.3
-pyerfa==2.0.0
-pyparsing==2.4.7
+numpy==1.22.3
+packaging==21.3
+pandas==1.4.2
+polare==0.1.0
+pyerfa==2.0.0.1
+pykwalify==1.8.0
+pyparsing==3.0.8
+pyrsistent==0.18.0
 python-dateutil==2.8.2
-pytz==2021.1
-requests==2.26.0
-scipy==1.7.1
+pytz==2022.1
+PyYAML==6.0
+ruamel.yaml==0.17.17
+ruamel.yaml.clib==0.2.6
+scipy==1.8.0
 six==1.16.0
-snowballstemmer==2.1.0
-Sphinx==4.2.0
-sphinx-rtd-theme==1.0.0
-sphinxcontrib-applehelp==1.0.2
-sphinxcontrib-devhelp==1.0.2
-sphinxcontrib-htmlhelp==2.0.0
-sphinxcontrib-jsmath==1.0.1
-sphinxcontrib-qthelp==1.0.3
-sphinxcontrib-serializinghtml==1.1.5
-urllib3==1.26.7
+spinnaker-python==2.4.0.143
+wincertstore==0.2
diff --git a/setup.py b/setup.py
index 466c5a2..1f78aa2 100644
--- a/setup.py
+++ b/setup.py
@@ -1,25 +1,47 @@
+"""Celest is a package for satellite mission planning.
+
+It provides:
+
+- Diverse arrangement of time and coordinate conversions,
+- Satellite orbital and attitude calculations, and
+- Satellite-to-ground encounter calculations.
+
+Refer to project links for documentation and source code urls.
+
+All Celest wheels distributed on PyPI are BSD licensed.
+"""
+DOCLINES = (__doc__ or '').split("\n")
+
+
 import setuptools
 
-with open("README.md", "r", encoding="utf-8") as fh:
-    long_description = fh.read()
+
+__version__ = "0.3.0"
+
 
 setuptools.setup(
     name="celest",
-    version="0.2.0",
-    author="Jai Willems",
+    version=__version__,
+    author="Jai Willems et al.",
     author_email="jai52h@hotmail.com",
     description="Satellite dynamics and mission planning library.",
-    long_description=long_description,
+    long_description="\n".join(DOCLINES),
     long_description_content_type="text/markdown",
     url="https://github.com/JaiWillems/Celest",
-    license="BSD-3-Clause",
+    download_url="https://pypi.org/project/celest/",
+    project_urls={
+        "Bug Tracker": "https://github.com/JaiWillems/celest/issues",
+        "Documentation": "https://celest.readthedocs.io/en/latest/",
+        "Source Code": "https://github.com/JaiWillems/celest"
+    },
+    license="BSD",
     packages=setuptools.find_packages(include=["celest", "celest.*"]),
     install_requires=[
         "jplephem",
         "julian",
         "numpy",
         "pandas",
-        "scipy"
+        "polare"
     ],
     classifiers=[
         "Development Status :: 5 - Production/Stable",


Test patch:
diff --git a/tests/test_core/test_interpolation.py b/tests/test_core/test_interpolation.py
deleted file mode 100644
index da7d286..0000000
--- a/tests/test_core/test_interpolation.py
+++ /dev/null
@@ -1,68 +0,0 @@
-"""Testing module for the `Interpolation` class."""
-
-
-from celest.core.interpolation import _interpolate
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-class TestInterpolation(TestCase):
-
-    def test_interp(self):
-        """Test `Interpolation._interp`."""
-
-        data = np.linspace(0, 99, 100)
-
-        for factor in [0, 1, 2]:
-
-            # factor = {0, 1, 2}, dt=0, indices=None
-            out_data = np.linspace(0, 99, factor * 100)
-            calc_out = _interpolate(data=data, factor=factor)
-
-            for i in range(calc_out.size):
-                with self.subTest(i=i):
-                    self.assertAlmostEqual(out_data[i], calc_out[i], delta=0.001)
-
-            for dt in [0, 1, 2]:
-
-                if factor == 1:
-                    alpha = 1
-                else:
-                    alpha = 0
-
-                # factor = {0, 1, 2}, dt={0, 1, 2}, indices=[[10, 11, ..., 20]]
-                indices = np.array([[i for i in range(10, 21)]])
-
-                out_1 = np.linspace(0, 9 - dt, 10 - dt)
-                num = factor * (10 + 2 * dt + alpha)
-                out_2 = np.linspace(10 - dt, 20 + dt, num)
-                out_3 = np.linspace(21 + dt, 99, 79 - dt)
-                out_data = np.concatenate((out_1, out_2, out_3))
-                calc_out = _interpolate(data=data, factor=factor, dt=dt, indices=indices)
-
-                for i in range(calc_out.size):
-                    with self.subTest(i=i):
-                        self.assertAlmostEqual(out_data[i], calc_out[i], delta=0.001)
-
-                # factor = {0, 1, 2}, dt={0, 1, 2}, indices=[[10, 11, ..., 20],
-                # [40, 41, ..., 50]]
-                indices = np.array([[i for i in range(10, 21)], [i for i in range(40, 51)]])
-
-                out_1 = np.linspace(0, 9 - dt, 10 - dt)
-                num = factor * (10 + 2 * dt + alpha)
-                out_2 = np.linspace(10 - dt, 20 + dt, num)
-                out_3 = np.linspace(21 + dt, 39 - dt, 19 - 2 * dt)
-                num = factor * (10 + 2 * dt + alpha)
-                out_4 = np.linspace(40 - dt, 50 + dt, num)
-                out_5 = np.linspace(51 + dt, 99, 49 - dt)
-                out_data = np.concatenate((out_1, out_2, out_3, out_4, out_5))
-                calc_out = _interpolate(data=data, factor=factor, dt=dt, indices=indices)
-
-                for i in range(calc_out.size):
-                    with self.subTest(i=i):
-                        self.assertAlmostEqual(out_data[i], calc_out[i], delta=0.001)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_encounter/test_encounter_math_utils.py b/tests/test_encounter/test_encounter_math_utils.py
deleted file mode 100644
index 72b6f27..0000000
--- a/tests/test_encounter/test_encounter_math_utils.py
+++ /dev/null
@@ -1,62 +0,0 @@
-"""Testing module for the encounter mathematical basis."""
-
-
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._encounter_math_utils import _analytical_encounter_ind
-from celest.satellite.coordinate import Coordinate
-from celest.satellite.time import Time
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-class TestEncounterMathUtils(TestCase):
-
-    def setUp(self):
-        """Test fixure for test method execution."""
-
-        fname = "tests/test_data/coordinate_validation_set.txt"
-        data = np.loadtxt(fname=fname, delimiter="\t", skiprows=1)
-
-        times = data[:, 0]
-        itrs = data[:, 10:]
-
-        self.timeData = Time(times, 2430000)
-        self.coor = Coordinate(itrs, "itrs", self.timeData)
-        self.length = data.shape[0]
-
-    def test_anaytical_encounter_ind(self):
-        """Test `analytical_encounter_ind`."""
-
-        # Set up shared parameters.
-        lat, lon = 43.6532, -79.3832
-        GEO_data = self.coor._geo_to_itrs(np.array([[lat, lon, 0]]))
-        gnd_itrs = np.repeat(GEO_data, self.length, 0)
-        sat_itrs = self.coor.itrs()
-
-        sat_itrs = self.coor.itrs()
-
-        # Test imaging encounter.
-        location = GroundPosition(lat, lon)
-
-        altitude, _ = self.coor.horizontal(location)
-        off_nadir = self.coor.off_nadir(location)
-
-        ind = np.where((0 <= altitude) & (off_nadir < 30))[0]
-        calc_ind = _analytical_encounter_ind(sat_itrs, gnd_itrs, 30, 1)
-
-        self.assertTrue(np.array_equal(ind, calc_ind))
-
-        # Test transmission encounter.
-        location = GroundPosition(lat, lon)
-
-        altitude, _ = self.coor.horizontal(location)
-
-        ind = np.where(30 <= altitude)[0]
-        calc_ind = _analytical_encounter_ind(sat_itrs, gnd_itrs, 30, 0)
-
-        self.assertTrue(np.array_equal(ind, calc_ind))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_encounter/test_groundposition.py b/tests/test_encounter/test_groundposition.py
index f21c481..8e7ee5b 100644
--- a/tests/test_encounter/test_groundposition.py
+++ b/tests/test_encounter/test_groundposition.py
@@ -1,4 +1,3 @@
-"""Test module for the `GroundPosition` class."""
 
 
 from celest.encounter.groundposition import GroundPosition
@@ -18,7 +17,7 @@ def setUp(self):
 
     def test_radius(self):
         """Test `GroundPosition._WGS84_radius`.
-        
+
         Notes
         -----
         Validation against online calculator methodology.[1]_
@@ -43,4 +42,4 @@ def test_radius(self):
 
 
 if __name__ == "__main__":
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/tests/test_encounter/test_window_handling.py b/tests/test_encounter/test_window_handling.py
index 04f4775..8b5b9f3 100644
--- a/tests/test_encounter/test_window_handling.py
+++ b/tests/test_encounter/test_window_handling.py
@@ -1,4 +1,3 @@
-"""Testing module for the window data structures."""
 
 
 from celest.encounter.groundposition import GroundPosition
@@ -11,49 +10,68 @@
 class TestWindow(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
 
         location = GroundPosition(43.6532, 79.3832)
-        self.window = Window("", location, 21282.4, 21282.5, "image", 30, 1, 7)
-
-    def test_str(self):
-        """Test `Window.__str__`."""
-
-        print(self.window)
+        self.window = Window("", location, 21282.4, 21282.5, "image", 30, 1)
 
     def test_copy(self):
-        """Test `Window.copy`."""
 
-        window_copy = self.window.copy()
+        wind_copy = self.window.copy()
 
-        self.assertTrue(self.window.satellite == window_copy.satellite)
-        self.assertTrue(self.window.coor == window_copy.coor)
-        self.assertTrue(self.window.type == window_copy.type)
-        self.assertTrue(self.window.start == window_copy.start)
-        self.assertTrue(self.window.end == window_copy.end)
+        self.assertTrue(self.window.satellite == wind_copy.satellite)
+        self.assertTrue(self.window.coor == wind_copy.coor)
+        self.assertTrue(self.window.type == wind_copy.type)
+        self.assertTrue(self.window.start == wind_copy.start)
+        self.assertTrue(self.window.end == wind_copy.end)
 
 
 class TestWindows(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
-
 
         location = GroundPosition(43.6532, 79.3832)
         self.windows = Windows()
 
-        self.window_one = Window("", location, 1319, 1319.2, "image", 30, 1, 0)
-        self.window_two = Window("", location, 1322.3, 1322.31, "data link", 10, -1, 30)
-        self.window_three = Window("", location, 1320, 1320.2, "image", 30, 1, 0)
-        self.window_three = Window("", location, 1320, 1320.3, "image", 30, 1, 0)
-        self.window_three = Window("", location, 1320, 1320.1, "image", 30, 1, 0)
-
-        self.windows._add_window(self.window_one)
-        self.windows._add_window(self.window_two)
-        self.windows._add_window(self.window_three)
+        self.wind_one = Window("", location, 1319, 1319.2, "image", 30, 1)
+        self.wind_two = Window("", location, 1322.3, 1322.31, "data link", 10, -1)
+        self.wind_three = Window("", location, 1324, 1324.2, "image", 30, 1)
+        self.wind_four = Window("", location, 1324.3, 1325.3, "image", 30, 1)
+        self.wind_five = Window("", location, 1325, 1326.1, "image", 30, 1)
+
+        self.windows._add_window(self.wind_one)
+        self.windows._add_window(self.wind_two)
+        self.windows._add_window(self.wind_three)
+        self.windows._add_window(self.wind_four)
+        self.windows._add_window(self.wind_five)
+
+    def test_getitem(self):
+
+        # Test int/float indexing.
+        self.assertEqual(self.windows[1310], self.wind_one)
+        self.assertEqual(self.windows[1319], self.wind_one)
+        self.assertEqual(self.windows[1323], self.wind_two)
+        self.assertEqual(self.windows[1324.4], self.wind_four)
+        self.assertEqual(self.windows[1330], self.wind_five)
+
+        # Test tuple indexing with unique mapping.
+        arr = [self.wind_two, self.wind_five]
+        val_wind = np.array(arr, dtype=object)
+        self.assertTrue(np.array_equiv(self.windows[1322.5, 1325], val_wind))
+
+        # Test tuple indexing without unique mapping.
+        self.assertEqual(self.windows[1323.5, 1324.1], self.wind_three)
+
+        # Test slice indexing without boundaries cases.
+        arr = [self.wind_three, self.wind_four, self.wind_five]
+        val_wind = np.array(arr, dtype=object)
+        self.assertTrue(np.array_equiv(self.windows[1322.5:1326], val_wind))
+
+        # Test slice indexing with boundary cases.
+        arr = [self.wind_two, self.wind_three, self.wind_four, self.wind_five]
+        val_wind = np.array(arr, dtype=object)
+        self.assertTrue(np.array_equiv(self.windows[1322.3:1325], val_wind))
 
     def test_add_window(self):
-        """Test `Window._add_window`."""
 
         test_val = 0
         for window in self.windows:
@@ -62,14 +80,12 @@ def test_add_window(self):
             test_val = start
 
     def test_stats(self):
-        """Test `Window.stats`."""
 
         df = self.windows.stats()
         self.assertIsNotNone(df)
 
     def test_to_numpy(self):
-        """Test `Window.to_numpy`."""
-        
+
         np_arr = self.windows.to_numpy()
         self.assertIsInstance(np_arr, np.ndarray)
 
@@ -78,4 +94,4 @@ def test_save_encounters(self):
 
 
 if __name__ == "__main__":
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/tests/test_encounter/test_window_utils.py b/tests/test_encounter/test_window_utils.py
deleted file mode 100644
index fd87bc3..0000000
--- a/tests/test_encounter/test_window_utils.py
+++ /dev/null
@@ -1,121 +0,0 @@
-"""Testing module for the window utility functions."""
-
-
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._window_utils import _get_ang, _sun_itrs, _window_encounter_ind
-from celest.satellite.coordinate import Coordinate
-from celest.satellite.satellite import Satellite
-from celest.satellite.time import Time
-from unittest import TestCase
-import numpy as np
-import unittest
-
-
-class TestWindowUtils(TestCase):
-
-    def setUp(self):
-        """Test fixure for test method execution."""
-
-        fname = "tests/test_data/coordinate_validation_set.txt"
-        data = np.loadtxt(fname=fname, delimiter="\t", skiprows=1)
-        times, itrs = data[:, 0], data[:, 10:]
-
-        timeData = Time(times, 2430000)
-        coor = Coordinate(itrs, "itrs", timeData)
-
-        self.finch = Satellite(coor)
-
-    def test_sun_itrs(self):
-        """Test `Encounter._sun_itrs`.
-
-        Notes
-        -----
-        Test cases are generated using the `Astropy` Python package.
-        """
-        
-        from astropy.coordinates import get_body, ITRS
-        from astropy import time, units
-
-        julData = np.array([2455368.75, 2459450.85, 2456293.5416666665])
-        astropy_time = time.Time(julData, format="jd")
-        
-        sun_pos = get_body("sun", astropy_time)
-        sun_pos.representation_type="cartesian"
-
-        itrsCoor = sun_pos.transform_to(ITRS(obstime=astropy_time))
-
-        x = itrsCoor.x.to(units.km).value
-        y = itrsCoor.y.to(units.km).value
-        z = itrsCoor.z.to(units.km).value
-
-        calc_sun_pos = _sun_itrs(julData)
-        self.assertTrue(np.allclose(x, calc_sun_pos[:, 0], rtol=0.05))
-        self.assertTrue(np.allclose(y, calc_sun_pos[:, 1], rtol=0.05))
-        self.assertTrue(np.allclose(z, calc_sun_pos[:, 2], rtol=0.05))
-
-    def test_get_ang(self):
-        """Test `Encounter._get_ang`."""
-
-        vec_one = np.array([[56, 92, 76], [9238, 8479, 9387], [2, 98, 23]])
-        vec_two = np.array([[36, 29, 38], [2703, 947, 8739], [9827, 921, 1]])
-        ang = np.array([16.28, 37, 83.65])
-
-        calc_ang = _get_ang(vec_one, vec_two)
-
-        for i in range(calc_ang.size):
-            with self.subTest(i=i):
-                self.assertAlmostEqual(ang[i], calc_ang[i], delta=0.01)
-
-    def test_window_encounter_ind(self):
-        """Test `Encounter._window_encounter_ind`."""
-
-        julData = self.finch.time.julian()
-
-        # Check data link index generation.
-        location = GroundPosition(43.6532, -79.3832)
-
-        calc_enc_ind = _window_encounter_ind(self.finch, location, 30, 1, 0, 1)
-
-        altitude, _ = self.finch.position.horizontal(location=location)
-        off_nadir = self.finch.position.off_nadir(location=location)
-
-        sun_itrs = _sun_itrs(julData)
-        sun_pos = Coordinate(sun_itrs, "itrs", self.finch.time)
-        sun_alt, _ = sun_pos.horizontal(location=location)
-
-        enc_ind = np.arange(0, julData.shape[0], 1)
-        
-        alt_ind = np.where(altitude >= 0)[0]
-        off_nadir_ind = np.where(off_nadir <= 30)
-        sun_ind = np.where(sun_alt >= 0)
-        
-        enc_ind = np.intersect1d(enc_ind, alt_ind)
-        enc_ind = np.intersect1d(enc_ind, off_nadir_ind)
-        enc_ind = np.intersect1d(enc_ind, sun_ind)
-
-        self.assertTrue(np.array_equal(calc_enc_ind, enc_ind))
-
-        # Check data link index generation.
-        location = GroundPosition(43.6532, -79.3832)
-
-        calc_enc_ind = _window_encounter_ind(self.finch, location, 30, 0, 30, 0)
-
-        altitude, _ = self.finch.position.horizontal(location=location)
-
-        enc_ind = np.arange(0, julData.shape[0], 1)
-        
-        alt_ind = np.where(30 <= altitude)[0]
-
-        sat_itrs = self.finch.position.itrs()
-
-        ground_GEO = [location.lat, location.lon]
-        ground_GEO = np.repeat(np.array([ground_GEO]), julData.size, axis=0)
-        gnd_itrs = Coordinate(ground_GEO, "geo", self.finch.time).itrs()
-
-        sca_angs = _get_ang(sat_itrs - gnd_itrs, sun_itrs - gnd_itrs)
-        sca_ind = np.where(sca_angs > 30)[0]
-        
-        enc_ind = np.intersect1d(enc_ind, alt_ind)
-        enc_ind = np.intersect1d(enc_ind, sca_ind)
-
-        self.assertTrue(np.array_equal(calc_enc_ind, enc_ind))
diff --git a/tests/test_encounter/test_windows.py b/tests/test_encounter/test_windows.py
index db7cb93..11ba727 100644
--- a/tests/test_encounter/test_windows.py
+++ b/tests/test_encounter/test_windows.py
@@ -1,11 +1,7 @@
-"""Testing module for the encounter generation function."""
 
 
 from celest.encounter.groundposition import GroundPosition
-from celest.encounter.windows import generate
-from celest.encounter._window_handling import Window, Windows
-from celest.encounter._window_utils import _window_encounter_ind
-from celest.satellite.coordinate import Coordinate
+from celest.encounter.windows import _get_ang, _sun_coor, generate
 from celest.satellite.satellite import Satellite
 from celest.satellite.time import Time
 from unittest import TestCase
@@ -13,113 +9,86 @@
 import unittest
 
 
-class TestEncounter(TestCase):
+class TestWindows(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
 
         fname = "tests/test_data/coordinate_validation_set.txt"
         data = np.loadtxt(fname=fname, delimiter="\t", skiprows=1)
         times, itrs = data[:, 0], data[:, 10:]
 
-        timeData = Time(times, 2430000)
-        coor = Coordinate(itrs, "itrs", timeData)
+        self.sat = Satellite(itrs, "itrs", times, 2430000)
 
-        self.finch = Satellite(coor)
+    def test_sun_coor(self):
+
+        from astropy.coordinates import get_body, ITRS
+        from astropy import time, units
+
+        julData = np.array([2455368.75, 2456293.5416666665, 2459450.85])
+        astropy_time = time.Time(julData, format="jd")
+
+        sun_pos = get_body("sun", astropy_time)
+        sun_pos.representation_type = "cartesian"
+
+        itrsCoor = sun_pos.transform_to(ITRS(obstime=astropy_time))
+
+        x = itrsCoor.x.to(units.km).value
+        y = itrsCoor.y.to(units.km).value
+        z = itrsCoor.z.to(units.km).value
+
+        calc_x, calc_y, calc_z = _sun_coor(julData).itrs(stroke=False)
+        self.assertTrue(np.allclose(x, calc_x, rtol=0.05))
+        self.assertTrue(np.allclose(y, calc_y, rtol=0.05))
+        self.assertTrue(np.allclose(z, calc_z, rtol=0.05))
+
+    def test_get_ang(self):
+        """Test `Encounter._get_ang`."""
+
+        vec_one = np.array([[56, 92, 76], [9238, 8479, 9387], [2, 98, 23]])
+        vec_two = np.array([[36, 29, 38], [2703, 947, 8739], [9827, 921, 1]])
+        ang = [16.28, 37, 83.65]
+
+        calc_ang = _get_ang(vec_one, vec_two)
+
+        for i in range(calc_ang.size):
+            with self.subTest(i=i):
+                self.assertAlmostEqual(ang[i], calc_ang[i], delta=0.01)
 
     def test_windows(self):
-        """Test `Encounter.windows.generate.`"""
-
-        location_1 = GroundPosition(43.6532, -79.3832)
-        location_2 = GroundPosition(52.1579, -106.6702)
-
-        # Test case 1.
-        windows_1 = generate(self.finch, location_1, "image", 30, 0)
-        indices = _window_encounter_ind(self.finch, location_1, 30, 1, 0, 1)
-        indices = np.split(indices, np.where(np.diff(indices) != 1)[0] + 1)
-
-        val_windows = Windows()
-        time = self.finch.time.julian()
-        for reg in indices:
-                
-            start = time[reg[0]]
-            end = time[reg[-1]]
-                
-            window = Window(self.finch, location_1, start, end, "image", 30, 1, 0)
-            val_windows._add_window(window)
-        
-        self.assertEqual(len(val_windows), len(windows_1))
-
-        windows_val, windows_calc = val_windows.windows, windows_1.windows
-        for val, calc in zip(windows_val, windows_calc):
-            self.assertEqual(val.start, calc.start)
-            self.assertEqual(val.end, calc.end)
-
-        # Test case 2.
-        windows_2 = generate(self.finch, location_1, "data_link", 10, 0)
-        indices = _window_encounter_ind(self.finch, location_1, 10, 0, 30, 0)
-        indices = np.split(indices, np.where(np.diff(indices) != 1)[0] + 1)
-
-        val_windows = Windows()
-        time = self.finch.time.julian()
-        for reg in indices:
-                
-            start = time[reg[0]]
-            end = time[reg[-1]]
-                
-            window = Window(self.finch, location_1, start, end, "data_link", 10, 0, 30)
-            val_windows._add_window(window)
-        
-        self.assertEqual(len(val_windows), len(windows_2))
-
-        windows_val, windows_calc = val_windows.windows, windows_2.windows
-        for val, calc in zip(windows_val, windows_calc):
-            self.assertEqual(val.start, calc.start)
-            self.assertEqual(val.end, calc.end)
-
-        windows_3 = generate(self.finch, location_2, "image", 30, 0)
-        indices = _window_encounter_ind(self.finch, location_2, 30, 1, 0, 1)
-        indices = np.split(indices, np.where(np.diff(indices) != 1)[0] + 1)
-
-        val_windows = Windows()
-        time = self.finch.time.julian()
-        for reg in indices:
-                
-            start = time[reg[0]]
-            end = time[reg[-1]]
-                
-            window = Window(self.finch, location_2, start, end, "image", 30, 1, 0)
-            val_windows._add_window(window)
-        
-        self.assertEqual(len(val_windows), len(windows_3))
-
-        windows_val, windows_calc = val_windows.windows, windows_3.windows
-        for val, calc in zip(windows_val, windows_calc):
-            self.assertEqual(val.start, calc.start)
-            self.assertEqual(val.end, calc.end)
-        
-
-        windows_4 = generate(self.finch, location_2, "data_link", 10)
-        indices = _window_encounter_ind(self.finch, location_2, 10, 0, 30, 0)
-        indices = np.split(indices, np.where(np.diff(indices) != 1)[0] + 1)
-
-        val_windows = Windows()
-        time = self.finch.time.julian()
-        for reg in indices:
-                
-            start = time[reg[0]]
-            end = time[reg[-1]]
-                
-            window = Window(self.finch, location_1, start, end, "data_link", 30, 1, 0)
-            val_windows._add_window(window)
-        
-        self.assertEqual(len(val_windows), len(windows_4))
-
-        windows_val, windows_calc = val_windows.windows, windows_4.windows
-        for val, calc in zip(windows_val, windows_calc):
-            self.assertEqual(val.start, calc.start)
-            self.assertEqual(val.end, calc.end)
+
+        location = GroundPosition(43.6532, -79.3832)
+
+        altitude = self.sat._altitude(location)
+        off_nadir = self.sat.off_nadir(location, stroke=True)
+
+        # Test the window generator for day imaging case.
+
+        enc, ang, lighting, tol = "image", 30, 1, 1e-5
+        windows = generate(self.sat, location, enc, ang, lighting, tol)
+
+        for window in windows:
+            start, end = window.start, window.end
+            window_times = np.linspace(start, end, 10)
+
+            tha = Time(window_times).true_hour_angle(location.lon)
+
+            self.assertTrue(np.all(altitude(window_times) > 0))
+            self.assertTrue(np.all(off_nadir(window_times) < ang))
+            self.assertTrue(np.all((-90 < tha) & (tha < 90)))
+
+        # Test the window generator for all day data link case.
+
+        enc, ang, lighting, tol = "data_link", 10, 0, 1e-5
+        windows = generate(self.sat, location, enc, ang, lighting, tol)
+
+        for window in windows:
+            start, end = window.start, window.end
+            window_times = np.linspace(start, end, 10)
+
+            tha = Time(window_times).true_hour_angle(location.lon)
+
+            self.assertTrue(np.all(altitude(window_times) > 10))
 
 
 if __name__ == "__main__":
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/tests/test_satellite/test_angle_representations.py b/tests/test_satellite/test_angle_representations.py
index 23fc784..f25dd12 100644
--- a/tests/test_satellite/test_angle_representations.py
+++ b/tests/test_satellite/test_angle_representations.py
@@ -1,4 +1,3 @@
-"""Testing module for the angular representations functionality."""
 
 
 from celest.satellite._angle_representations import sexagesimal, _ISO6709_representation
@@ -12,32 +11,27 @@ def test_sexagesimal(self):
         """Test `Coordinate.sexagesimal`."""
 
         ang = np.array([43.6532, -79.3832, -33.2833, 149.1000])
-        true_sexagesimal = np.array(["+43\u00B039\u203211.52\u2033",
-                                "-79\u00B022\u203259.52\u2033",
-                                "-33\u00B016\u203259.88\u2033",
-                                "+149\u00B006\u203200.00\u2033"])
+        true_sexagesimal = ["+43\u00B039\u203211.52\u2033",
+                            "-79\u00B022\u203259.52\u2033",
+                            "-33\u00B016\u203259.88\u2033",
+                            "+149\u00B006\u203200.00\u2033"]
 
         calc_ang = sexagesimal(ang)
 
         for i in range(calc_ang.shape[0]):
             with self.subTest(i=i):
                 self.assertTrue(true_sexagesimal[i] == calc_ang[i])
-    
-    def test_ISO6709_representation(self):
-        """Test `Coordinate._ISO6709_representation`.
 
-        Notes
-        -----
-        Test cases generated in accordance to ISO6709 formatting standards.
-        """
+    def test_ISO6709_representation(self):
+        """Test `Coordinate._ISO6709_representation`."""
 
         geo = np.array([[43.6532, -79.3832, 430.23],
                         [-33.2833, 149.1000, 532.98]])
-        true_sexagesimal = np.array(["43\u00B039\u203211.52\u2033N 79\u00B022\u203259.52\u2033W 430.23km",
-                                "33\u00B016\u203259.88\u2033S 149\u00B006\u203200.00\u2033E 532.98km"])
+        true_sexagesimal = ["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_ang = _ISO6709_representation(geo)
 
         for i in range(calc_ang.shape[0]):
             with self.subTest(i=i):
-                self.assertTrue(true_sexagesimal[i] == calc_ang[i])
\ No newline at end of file
+                self.assertTrue(true_sexagesimal[i] == calc_ang[i])
diff --git a/tests/test_satellite/test_astronomical_quantities.py b/tests/test_satellite/test_astronomical_quantities.py
index 90541b2..acf8190 100644
--- a/tests/test_satellite/test_astronomical_quantities.py
+++ b/tests/test_satellite/test_astronomical_quantities.py
@@ -1,4 +1,3 @@
-"""Testing module for the `AstronomicalQuantities` calculations."""
 
 
 from celest.satellite._astronomical_quantities import *
@@ -194,7 +193,7 @@ def test_equation_of_equinoxes(self):
         for i in range(julian.size-1):
             with self.subTest(i=i):
                 self.assertAlmostEqual(EOE[i], true_EOE[i], delta=0.0001)
-    
+
     def test_sun_right_ascension(self):
         """Test `AstronomicalQuantities.sun_right_ascension`.
 
diff --git a/tests/test_satellite/test_coordinate.py b/tests/test_satellite/test_coordinate.py
index ebb37bb..c01d4bc 100644
--- a/tests/test_satellite/test_coordinate.py
+++ b/tests/test_satellite/test_coordinate.py
@@ -1,9 +1,8 @@
-"""Testing module for the `Coordinate` class."""
 
 
 from celest.encounter.groundposition import GroundPosition
 from celest.satellite.coordinate import Coordinate
-from celest.satellite.time import Time
+from polare import Stroke
 from unittest import TestCase
 import numpy as np
 import unittest
@@ -12,7 +11,6 @@
 class TestCoordinate(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
 
         fname = "tests/test_data/coordinate_validation_set.txt"
         data = np.loadtxt(fname=fname, delimiter="\t", skiprows=1)
@@ -23,134 +21,121 @@ def setUp(self):
         self.local_altaz = data[:, 4:7]
         self.GCRS = data[:, 7:10]
         self.ITRS = data[:, 10:]
-        self.length = data.shape[0]
+        self._length = data.shape[0]
 
-        geo = np.concatenate((self.geo, self.alt.reshape((-1, 1))), axis=1)
+        self.tri_geo = np.concatenate((self.geo, self.alt.reshape((-1, 1))), axis=1)
 
-        self.timeData = Time(self.times, 2430000)
-        self.coor_gcrs = Coordinate(self.GCRS, "gcrs", self.timeData)
-        self.coor_itrs = Coordinate(self.ITRS, "itrs", self.timeData)
-        self.coor_geo = Coordinate(geo, "geo", self.timeData)
+        self.offset = 2430000
+        self.coor_gcrs = Coordinate(self.GCRS, "gcrs", self.times, self.offset)
+        self.coor_itrs = Coordinate(self.ITRS, "itrs", self.times, self.offset)
+        self.coor_geo = Coordinate(self.tri_geo, "geo", self.times, self.offset)
 
     def test_set_base_position(self):
         """Test `Coordinate._set_base_position`."""
 
-        coor_1 = Coordinate(self.geo, "geo", self.timeData)
+        coor_1 = Coordinate(self.geo, "geo", self.times, self.offset)
 
-        self.assertEqual(coor_1._GEO.shape[1], 3)
-
-        self.assertIsNotNone(coor_1.time)
-        self.assertIsNotNone(coor_1._GEO)
+        self.assertIsNotNone(coor_1._julian)
         self.assertIsNone(coor_1._GCRS)
-        self.assertIsNone(coor_1._ITRS)
-        self.assertEqual(self.length, coor_1.length)
+        self.assertIsNotNone(coor_1._ITRS)
+        self.assertEqual(self._length, coor_1._length)
 
-        coor_2 = Coordinate(self.GCRS, "gcrs", self.timeData)
+        coor_2 = Coordinate(self.GCRS, "gcrs", self.times, self.offset)
 
-        self.assertIsNotNone(coor_2.time)
-        self.assertIsNone(coor_2._GEO)
+        self.assertIsNotNone(coor_2._julian)
         self.assertIsNotNone(coor_2._GCRS)
         self.assertIsNone(coor_2._ITRS)
-        self.assertEqual(self.length, coor_2.length)
+        self.assertEqual(self._length, coor_2._length)
 
-        coor_3 = Coordinate(self.ITRS, "itrs", self.timeData)
+        coor_3 = Coordinate(self.ITRS, "itrs", self.times, self.offset)
 
-        self.assertIsNotNone(coor_3.time)
-        self.assertIsNone(coor_3._GEO)
+        self.assertIsNotNone(coor_3._julian)
         self.assertIsNone(coor_3._GCRS)
         self.assertIsNotNone(coor_3._ITRS)
-        self.assertEqual(self.length, coor_3.length)
+        self.assertEqual(self._length, coor_3._length)
 
     def test_geo_to_itrs(self):
         """Test `Coordinate._geo_to_itrs`.
 
         Notes
         -----
-        Test cases are taken from a GMAT data set.
+        Test cases generated using a GMAT data set.
         """
 
-        geo = np.concatenate((self.geo, self.alt.reshape((-1, 1))), axis=1)
-        calc_itrs = self.coor_geo._geo_to_itrs(geo)
-
-        for i in range(self.length):
-            with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs[i, 0], self.ITRS[i, 0], delta=0.001)
-
-        for i in range(self.length):
-            with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs[i, 1], self.ITRS[i, 1], delta=0.001)
+        calc_itrs = self.coor_geo._geo_to_itrs(self.tri_geo)
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs[i, 2], self.ITRS[i, 2], delta=0.001)
+                self.assertAlmostEqual(calc_itrs[i, 0], self.ITRS[i, 0],
+                                       delta=0.001)
+                self.assertAlmostEqual(calc_itrs[i, 1], self.ITRS[i, 1],
+                                       delta=0.001)
+                self.assertAlmostEqual(calc_itrs[i, 2], self.ITRS[i, 2],
+                                       delta=0.001)
 
     def test_itrs_to_geo(self):
         """Test `Coordinate._itrs_to_geo`.
 
         Notes
         -----
-        Test cases are taken from a GMAT data set.
+        Test cases generated using a GMAT data set.
         """
 
-        calc_geo = self.coor_itrs._itrs_to_geo(self.ITRS)
-
-        for i in range(self.length):
-            with self.subTest(i=i):
-                self.assertAlmostEqual(calc_geo[i, 0], self.geo[i, 0], delta=0.18)
-
-        for i in range(self.length):
-            with self.subTest(i=i):
-                self.assertAlmostEqual(calc_geo[i, 1], self.geo[i, 1], delta=0.00001)
+        itrs = np.array([Stroke(self.times, self.ITRS[:, i]) for i in range(3)])
+        calc_lat, calc_lon = self.coor_itrs._itrs_to_geo(itrs)
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_geo[i, 2], self.alt[i], delta=0.06)
+                self.assertAlmostEqual(calc_lat(self.times[i]), self.geo[i, 0],
+                                       delta=0.18)
+                self.assertAlmostEqual(calc_lon(self.times[i]), self.geo[i, 1],
+                                       delta=0.00001)
 
     def test_geo(self):
         """Test `Coordinate.geo`."""
 
-        calc_geo_1 = self.coor_gcrs.geo()
-        calc_geo_2 = self.coor_itrs.geo()
+        lat1, lon1, alt1 = self.coor_gcrs.geo()
+        lat2, lon2, alt2 = self.coor_itrs.geo()
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
 
-                epsilon = 1
-                calc_lon, test_lon = calc_geo_1[i, 1], self.geo[i, 1]
-                condition_one = abs(180 - calc_lon) < epsilon and abs(180 + test_lon) < epsilon
-                condition_two = abs(180 + calc_lon) < epsilon and abs(180 - test_lon) < epsilon
-                if condition_one:
+                eps = 1
+                calc_lon, test_lon = lon1[i], self.geo[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 condition_two:
+                elif cond_two:
                     calc_lon = (calc_lon + 360) % 360
 
-                self.assertAlmostEqual(calc_geo_1[i, 0], self.geo[i, 0], delta=0.5)
+                self.assertAlmostEqual(lat1[i], self.geo[i, 0], delta=0.5)
                 self.assertAlmostEqual(calc_lon, test_lon, delta=0.9)
-                self.assertAlmostEqual(calc_geo_1[i, 2], self.alt[i], delta=0.1)
+                self.assertAlmostEqual(alt1[i], self.alt[i], delta=0.1)
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_geo_2[i, 0], self.geo[i, 0], delta=0.18)
-                self.assertAlmostEqual(calc_geo_2[i, 1], self.geo[i, 1], delta=0.00001)
-                self.assertAlmostEqual(calc_geo_2[i, 2], self.alt[i], delta=0.06)
+                self.assertAlmostEqual(lat2[i], self.geo[i, 0], delta=0.18)
+                self.assertAlmostEqual(lon2[i], self.geo[i, 1], delta=0.00001)
+                self.assertAlmostEqual(alt2[i], self.alt[i], delta=0.06)
 
     def test_era(self):
         """Test `Coordinate.era`.
 
         Notes
         -----
-        The `Coordinate.era` method was validated through
-        `Coordinate._gcrs_and_itrs` validation. This test method was implemented
-        using the calculated output of the `Coordinate.era` method after it was
-        shown to be correct.
+        `Coordinate.era` was validated through `Coordinate._gcrs_and_itrs`
+        validation. Test cases were generated from the `Coordinate.era` output
+        after it was shown to be correct to ensure functionality is consistent.
         """
 
         era = np.degrees(np.array([6.2360075]))
 
-        timeData = Time(np.array([2454545]))
-        basePos = np.array([[6343.82, -2640.87, -11.26]])
-        posData = Coordinate(basePos, "itrs", timeData)
-        calc_era = posData.era()
+        julian = [2454545]
+        position = [[6343.82, -2640.87, -11.26]]
+        coor = Coordinate(position, "itrs", julian)
+        calc_era = coor.era()
 
         self.assertAlmostEqual(era[0], calc_era[0], delta=0.01)
 
@@ -159,59 +144,77 @@ def test_gcrs_and_itrs(self):
 
         Notes
         -----
-        Test cases are taken from a GMAT data set.
+        Test cases generated using a GMAT data set.
         """
 
         calc_itrs = self.coor_gcrs._gcrs_and_itrs(self.GCRS, frame="gcrs")
         calc_gcrs = self.coor_itrs._gcrs_and_itrs(self.ITRS, frame="itrs")
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs[i, 0], self.ITRS[i, 0], delta=20)
-                self.assertAlmostEqual(calc_itrs[i, 1], self.ITRS[i, 1], delta=20)
-                self.assertAlmostEqual(calc_itrs[i, 2], self.ITRS[i, 2], delta=20)
-
-        for i in range(self.length):
+                self.assertAlmostEqual(calc_itrs[i, 0], self.ITRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs[i, 1], self.ITRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs[i, 2], self.ITRS[i, 2],
+                                       delta=0.35)
+
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_gcrs[i, 0], self.GCRS[i, 0], delta=20)
-                self.assertAlmostEqual(calc_gcrs[i, 1], self.GCRS[i, 1], delta=20)
-                self.assertAlmostEqual(calc_gcrs[i, 2], self.GCRS[i, 2], delta=20)
+                self.assertAlmostEqual(calc_gcrs[i, 0], self.GCRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs[i, 1], self.GCRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs[i, 2], self.GCRS[i, 2],
+                                       delta=0.35)
 
     def test_gcrs(self):
         """Test `Coordinate.gcrs`."""
 
-        calc_gcrs_1 = self.coor_itrs.gcrs()
-        calc_gcrs_2 = self.coor_geo.gcrs()
+        calc_gcrs_x1, calc_gcrs_y1, calc_gcrs_z1 = self.coor_itrs.gcrs()
+        calc_gcrs_x2, calc_gcrs_y2, calc_gcrs_z2 = self.coor_geo.gcrs()
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_gcrs_1[i, 0], self.GCRS[i, 0], delta=17)
-                self.assertAlmostEqual(calc_gcrs_1[i, 1], self.GCRS[i, 1], delta=17)
-                self.assertAlmostEqual(calc_gcrs_1[i, 2], self.GCRS[i, 2], delta=17)
-
-        for i in range(self.length):
+                self.assertAlmostEqual(calc_gcrs_x1[i], self.GCRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs_y1[i], self.GCRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs_z1[i], self.GCRS[i, 2],
+                                       delta=0.35)
+
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_gcrs_2[i, 0], self.GCRS[i, 0], delta=17)
-                self.assertAlmostEqual(calc_gcrs_2[i, 1], self.GCRS[i, 1], delta=17)
-                self.assertAlmostEqual(calc_gcrs_2[i, 2], self.GCRS[i, 2], delta=17)
+                self.assertAlmostEqual(calc_gcrs_x2[i], self.GCRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs_y2[i], self.GCRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_gcrs_z2[i], self.GCRS[i, 2],
+                                       delta=0.35)
 
     def test_itrs(self):
         """Test `Coordinate.itrs`."""
 
-        calc_itrs_1 = self.coor_gcrs.itrs()
-        calc_itrs_2 = self.coor_geo.itrs()
+        calc_itrs_x1, calc_itrs_y1, calc_itrs_z1 = self.coor_gcrs.itrs()
+        calc_itrs_x2, calc_itrs_y2, calc_itrs_z2 = self.coor_geo.itrs()
 
-        for i in range(self.length):
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs_1[i, 0], self.ITRS[i, 0], delta=17)
-                self.assertAlmostEqual(calc_itrs_1[i, 1], self.ITRS[i, 1], delta=17)
-                self.assertAlmostEqual(calc_itrs_1[i, 2], self.ITRS[i, 2], delta=17)
-
-        for i in range(self.length):
+                self.assertAlmostEqual(calc_itrs_x1[i], self.ITRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs_y1[i], self.ITRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs_z1[i], self.ITRS[i, 2],
+                                       delta=0.35)
+
+        for i in range(self._length):
             with self.subTest(i=i):
-                self.assertAlmostEqual(calc_itrs_2[i, 0], self.ITRS[i, 0], delta=5)
-                self.assertAlmostEqual(calc_itrs_2[i, 1], self.ITRS[i, 1], delta=5)
-                self.assertAlmostEqual(calc_itrs_2[i, 2], self.ITRS[i, 2], delta=5)
+                self.assertAlmostEqual(calc_itrs_x2[i], self.ITRS[i, 0],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs_y2[i], self.ITRS[i, 1],
+                                       delta=0.35)
+                self.assertAlmostEqual(calc_itrs_z2[i], self.ITRS[i, 2],
+                                       delta=0.35)
 
     def test_get_ang(self):
         """Test `Coordinate._get_ang`."""
@@ -231,7 +234,7 @@ def test_horizontal(self):
 
         Notes
         -----
-        Test cases are generated using the `Astropy` Python package.
+        Test cases generated using the `Astropy` Python package.
         """
 
         from astropy.coordinates import SkyCoord, ITRS, GCRS, EarthLocation, AltAz
@@ -243,16 +246,17 @@ def test_horizontal(self):
         loc = EarthLocation.from_geodetic(lon*u.deg, lat*u.deg)
 
         # Prepare time and position information.
-        timeData = time.Time(self.times + 2430000, format="jd")
-        x, y, z = self.GCRS[:, 0], self.GCRS[:, 1], self.GCRS[:, 2]
+        times = time.Time(self.times + self.offset, format="jd")
+        x, y, z = self.GCRS.T
 
         # Define coordinate frames.
-        gcrs = GCRS(obstime=timeData)
-        itrs = ITRS(obstime=timeData)
-        altaz = AltAz(obstime=timeData, location=loc)
+        gcrs = GCRS(obstime=times)
+        itrs = ITRS(obstime=times)
+        altaz = AltAz(obstime=times, location=loc)
 
         # Convert between frames.
-        gcrsCoor = SkyCoord(x=x, y=y, z=z, unit='km', frame=gcrs, representation_type='cartesian')
+        gcrsCoor = SkyCoord(x=x, y=y, z=z, unit='km',
+                            frame=gcrs, representation_type='cartesian')
         itrsCoor = gcrsCoor.transform_to(itrs)
         altazCoor = itrsCoor.transform_to(altaz)
 
@@ -262,16 +266,13 @@ def test_horizontal(self):
         az = altazCoor.az.degree
 
         # Get Celest results.
-        coor = Coordinate(itrsData, "itrs", self.timeData)
+        coor = Coordinate(itrsData, "itrs", self.times, self.offset)
         groundPos = GroundPosition(lat, lon)
         calc_alt, calc_az = coor.horizontal(groundPos)
 
-        for i in range(calc_alt.size-5000):
+        for i in range(calc_alt.size):
             with self.subTest(i=i):
                 self.assertAlmostEqual(alt[i], calc_alt[i], delta=0.25)
-
-        for i in range(calc_az.size-5000):
-            with self.subTest(i=i):
                 self.assertAlmostEqual(az[i], calc_az[i], delta=7.5)
 
     def test_off_nadir(self):
@@ -279,27 +280,27 @@ def test_off_nadir(self):
 
         Notes
         -----
-        Test cases generated from geometric model designed by Mingde Yin.
+        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 number of significant figures and coarse parameter
-        selection of the geometric model used to generated test cases.
+        due to the low precision in the model's output and parameters used to
+        generate the test cases.
         """
 
         off_nadir = np.array([66.88, 65.09, 63.90, 63.22, 62.46, 61.67, 58.42,
                               38.27, 23.73, 56.29])
 
-        groundPos = GroundPosition(52.1579, -106.6702)
+        location = GroundPosition(52.1579, -106.6702)
 
-        timeData = Time(self.times[210:220], 2430000)
-        coor = Coordinate(self.ITRS[210:220], "itrs", timeData)
-        calc_off_nadir = coor.off_nadir(groundPos)
+        julian = self.times[210:220]
+        coor = Coordinate(self.ITRS[210:220], "itrs", julian, self.offset)
+        calc_off_nadir = coor.off_nadir(location)
 
         for i in range(10):
             with self.subTest(i=i):
                 self.assertAlmostEqual(off_nadir[i], calc_off_nadir[i], delta=0.3)
 
-    def test_wgs84_radius(self):
+    def test_WGS84_radius(self):
         """Test `Coordinate._WGS84_radius`.
 
         Notes
@@ -312,19 +313,21 @@ def test_wgs84_radius(self):
            https://planetcalc.com/7721/.
         """
 
-        a = 6378.1370
-        b = 6356.7523142
+        a, b = 6378.1370, 6356.7523142
         lat = np.radians(self.geo[:, 0])
         clat, slat = np.cos(lat), np.sin(lat)
         num = (a ** 2 * clat) ** 2 + (b ** 2 * slat) ** 2
         denom = (a * clat) ** 2 + (b * slat) ** 2
         radius = np.sqrt(num / denom)
 
-        calc_radius = self.coor_geo._WGS84_radius(self.geo[:, 0])
+        lat = Stroke(self.times, self.geo[:, 0], "cubic")
 
-        for i in range(5):
+        calc_radius = self.coor_geo._WGS84_radius(lat)
+
+        for i in range(radius.size):
             with self.subTest(i=i):
-                self.assertAlmostEqual(radius[i], calc_radius[i], delta=0.001)
+                self.assertAlmostEqual(radius[i], calc_radius(self.times[i]),
+                                       delta=0.001)
 
     def test_altitude(self):
         """Test `Coordinate.altitude`.
@@ -343,21 +346,20 @@ def test_altitude(self):
     def test_distance(self):
         """Test `Coordinate.distance`."""
 
-        times = np.array([30462.5, 30462.50069444, 30462.50171802,
-                          30462.50278711, 30462.50386162])
-        position = np.array([[6343.81620221, -2640.87223125, -11.25541802],
-                             [6295.64583763, -2718.09271472, 443.08232543],
-                             [6173.04658005, -2808.91831102, 1108.5854422],
-                             [5980.91111229, -2872.96257946, 1792.17964249],
-                             [5724.3020284, -2904.09809986, 2460.25799377]])
-        dist = np.array([9070.49268746, 8776.7179543, 8330.99543153,
-                         7851.70082642, 7359.09189844])
-
-        groundPos = GroundPosition(52.1579, -106.6702)
-
-        timeData = Time(times, 2430000)
-        coor = Coordinate(position, "itrs", timeData)
-        calc_dist = coor.distance(groundPos)
+        julian = [30462.5, 30462.50069444, 30462.50171802, 30462.50278711,
+                  30462.50386162]
+        position = [[6343.81620221, -2640.87223125, -11.25541802],
+                    [6295.64583763, -2718.09271472, 443.08232543],
+                    [6173.04658005, -2808.91831102, 1108.5854422],
+                    [5980.91111229, -2872.96257946, 1792.17964249],
+                    [5724.3020284, -2904.09809986, 2460.25799377]]
+        dist = [9070.49268746, 8776.7179543, 8330.99543153, 7851.70082642,
+                7359.09189844]
+
+        location = GroundPosition(52.1579, -106.6702)
+
+        coor = Coordinate(position, "itrs", julian, self.offset)
+        calc_dist = coor.distance(location)
 
         for i in range(5):
             with self.subTest(i=i):
diff --git a/tests/test_satellite/test_satellite.py b/tests/test_satellite/test_satellite.py
index 2236edb..e58b864 100644
--- a/tests/test_satellite/test_satellite.py
+++ b/tests/test_satellite/test_satellite.py
@@ -1,11 +1,6 @@
-"""Testing module for the `Satellite` class."""
 
 
-from celest.encounter.groundposition import GroundPosition
-from celest.encounter._window_handling import Window, Windows
-from celest.satellite.coordinate import Coordinate
 from celest.satellite.satellite import Satellite
-from celest.satellite.time import Time
 from unittest import TestCase
 import numpy as np
 import unittest
@@ -14,7 +9,6 @@
 class TestSatellite(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
 
         fname = "tests/test_data/coordinate_validation_set.txt"
         data = np.loadtxt(fname=fname, delimiter="\t", skiprows=1)
@@ -22,59 +16,48 @@ def setUp(self):
         self.times = data[:, 0]
         self.ITRS = data[:, 10:]
 
-        self.timeData = Time(self.times, 2430000)
-        self.coor = Coordinate(self.ITRS, "itrs", self.timeData)
+        self.offset = 2430000
+        self.finch = Satellite(self.ITRS, "itrs", self.times, self.offset)
 
-        self.finch = Satellite(self.coor)
-
-    def test_interpolate(self):
-        """Test `Satellite.interpolate`."""
-
-        location = GroundPosition(0, 0)
-
-        window = Window(None, location, self.times[100] + 2430000, self.times[150] + 2430000, None, None, None, None)
-        window_list = Windows()
-        window_list._add_window(window)
-
-        len_i = len(self.finch)
-
-        self.finch.interpolate(windows=window_list, factor=5)
-
-        len_f = len(self.finch)
-
-        self.assertGreater(len_f, len_i)
+        self.julian = self.finch.julian()
+        self.ut1 = self.finch.ut1()
+        self.gmst = self.finch.gmst()
+        self.gast = self.finch.gast()
+        self.lat, self.lon, self.alt = self.finch.geo()
+        self.gcrs_x, self.gcrs_y, self.gcrs_z = self.finch.gcrs()
+        self.itrs_x, self.itrs_y, self.itrs_z = self.finch.itrs()
 
     def test_save_data(self):
         """Test `Satellite.save_data`."""
 
         times = ("julian", "ut1", "gmst", "gast")
         positions = ("geo", "gcrs", "itrs")
-        self.finch.save_data("test_data.csv", ",", times, positions)
+        self.finch.save_data(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_geo = data[:, 5:8]
+        load_lat = data[:, 5]
+        load_lon = data[:, 6]
+        load_alt = data[:, 7]
         load_gcrs = data[:, 8:11]
         load_itrs = data[:, 11:]
 
-        julian = self.finch.time.julian()
-        ut1 = self.finch.time.ut1()
-        gmst = self.finch.time.gmst()
-        gast = self.finch.time.gast()
-        geo = self.finch.position.geo()
-        gcrs = self.finch.position.gcrs()
-        itrs = self.finch.position.itrs()
-
-        self.assertTrue(np.array_equal(load_julian, julian))
-        self.assertTrue(np.array_equal(load_ut1, ut1))
-        self.assertTrue(np.array_equal(load_gmst, gmst))
-        self.assertTrue(np.array_equal(load_gast, gast))
-        self.assertTrue(np.array_equal(load_geo, geo))
-        self.assertTrue(np.array_equal(load_gcrs, gcrs))
-        self.assertTrue(np.array_equal(load_itrs, itrs))
+        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.lat))
+        self.assertTrue(np.array_equal(load_lon, self.lon))
+        self.assertTrue(np.array_equal(load_alt, self.alt))
+        self.assertTrue(np.array_equal(load_gcrs[:, 0], self.gcrs_x))
+        self.assertTrue(np.array_equal(load_gcrs[:, 1], self.gcrs_y))
+        self.assertTrue(np.array_equal(load_gcrs[:, 2], self.gcrs_z))
+        self.assertTrue(np.array_equal(load_itrs[:, 0], self.itrs_x))
+        self.assertTrue(np.array_equal(load_itrs[:, 1], self.itrs_y))
+        self.assertTrue(np.array_equal(load_itrs[:, 2], self.itrs_z))
 
         import os
         os.remove("test_data.csv")
diff --git a/tests/test_satellite/test_time.py b/tests/test_satellite/test_time.py
index 206c4b7..4b58bbf 100644
--- a/tests/test_satellite/test_time.py
+++ b/tests/test_satellite/test_time.py
@@ -1,4 +1,3 @@
-"""Testing module for the `Time` class."""
 
 
 from astropy import coordinates, time
@@ -12,9 +11,8 @@
 class TestTime(TestCase):
 
     def setUp(self):
-        """Test fixure for test method execution."""
 
-        self.julData = np.array([2455368.75, 2459450.85, 2456293.5416666665])
+        self.julData = [2455368.75, 2459450.85, 2456293.5416666665]
         self.astropy_time = time.Time(self.julData, format="jd")
 
     def test_true_solar_time(self):
@@ -30,11 +28,12 @@ def test_true_solar_time(self):
         .. [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.
+        .. [2] Time Conventions. url:
+           http://www.bom.gov.au/climate/data-services/solar/content/data-time.html.
         """
 
-        julian = np.array([2455368.75, 2459450.85, 2456293.5416666665])
-        lon = np.array([-105, -118.24, 147.46])
+        julian = [2455368.75, 2459450.85, 2456293.5416666665]
+        lon = [-105, -118.24, 147.46]
         tst = np.array([23.0715, 0.378059, 10.76556])
 
         calc_tst = Time(julian=julian).true_solar_time(longitude=lon)
@@ -52,11 +51,12 @@ def test_mean_solar_time(self):
 
         References
         ----------
-        .. [1] Time Conventions. url: http://www.bom.gov.au/climate/data-services/solar/content/data-time.html.
+        .. [1] Time Conventions. url:
+           http://www.bom.gov.au/climate/data-services/solar/content/data-time.html.
         """
 
-        julData = np.array([2456293.5416666665])
-        lon = np.array([147.46])
+        julData = 2456293.5416666665
+        lon = 147.46
         tst = np.array([10.83056])
 
         calc_tst = Time(julian=julData).mean_solar_time(longitude=lon)
@@ -79,15 +79,16 @@ def test_true_hour_angle(self):
            https://gml.noaa.gov/grad/solcalc/.
         """
 
-        julData = np.array([2455368.75, 2459450.85, 2456293.5416666665])
-        lon = np.array([-105, -118.24, 147.46])
+        julData = [2455368.75, 2459450.85, 2456293.5416666665]
+        lon = [-105, -118.24, 147.46]
         hour_angle = np.array([166.0734, 185.671, 341.53]) / 15 % 24
 
         calc_hour_angle = Time(julian=julData).true_hour_angle(longitude=lon)
 
         for i in range(calc_hour_angle.size):
             with self.subTest(i=i):
-                self.assertAlmostEqual(hour_angle[i], calc_hour_angle[i], delta=1.51)
+                self.assertAlmostEqual(hour_angle[i], calc_hour_angle[i],
+                                       delta=1.51)
 
     def test_mean_hour_angle(self):
         """Test `Time.true_hour_angle`.
@@ -98,15 +99,16 @@ def test_mean_hour_angle(self):
         angle.
         """
 
-        julData = np.array([2456293.5416666665])
-        lon = np.array([147.46])
+        julData = 2456293.5416666665
+        lon = 147.46
         hour_angle = np.array([342.4584]) / 15 % 24
 
         calc_hour_angle = Time(julian=julData).mean_hour_angle(longitude=lon)
 
         for i in range(calc_hour_angle.size):
             with self.subTest(i=i):
-                self.assertAlmostEqual(hour_angle[i], calc_hour_angle[i], delta=0.01)
+                self.assertAlmostEqual(hour_angle[i], calc_hour_angle[i],
+                                       delta=0.01)
 
     def test_ut1(self):
         """Test `Time.UT1`.
@@ -149,7 +151,8 @@ def test_datetime(self):
                 self.assertEqual(calc_datetime[i].month, dt.month)
                 self.assertEqual(calc_datetime[i].day, dt.day)
                 self.assertEqual(calc_datetime[i].second, dt.second)
-                self.assertAlmostEqual(calc_datetime[i].microsecond, dt.microsecond, delta=1)
+                self.assertAlmostEqual(calc_datetime[i].microsecond,
+                                       dt.microsecond, delta=1)
 
     def test_gmst(self):
         """Test `Time.gmst`.