From 8aa39f714252298f8e6b7a7da5855f9400ff9bae Mon Sep 17 00:00:00 2001
From: thetazero <lesha.seletskiy@gmail.com>
Date: Wed, 10 May 2023 15:39:50 -0400
Subject: [PATCH] Replace Array with AbstractArray

---
 src/astrodynamics.jl                |  12 ++--
 src/attitude.jl                     |  30 ++++----
 src/coordinates.jl                  |  90 ++++++++++++------------
 src/earth_environment/nrlmsise00.jl |  66 +++++++++---------
 src/orbit_dynamics.jl               | 102 ++++++++++++++--------------
 src/sgp_models.jl                   |   2 +-
 src/simulation/integrators.jl       |  12 ++--
 src/simulation/propagators.jl       |  20 +++---
 src/universe.jl                     |   2 +-
 test/simulation/test_integrators.jl |   2 +-
 10 files changed, 169 insertions(+), 169 deletions(-)

diff --git a/src/astrodynamics.jl b/src/astrodynamics.jl
index b96bdf7..86efdea 100644
--- a/src/astrodynamics.jl
+++ b/src/astrodynamics.jl
@@ -183,14 +183,14 @@ The osculating elements are assumed to be (in order):
 6. _M_, Mean anomaly [rad]
 
 Arguments:
-- x_oe `x::Array{<:Real, 1}`: Osculating orbital elements. See above for desription of the elements and their required order.
+- x_oe `x::AbstractArray{<:Real, 1}`: Osculating orbital elements. See above for desription of the elements and their required order.
 - `use_degrees:Bool`: If `true` interpret input will be interpreted as being in degrees, and output will be returned in degrees.
 - `GM::Real`: Gravitational constant of central body. Defaults to `SatelliteDynamics.GM_EARTH` if none is provided.
 
 # Returns
-- x `x::Array{<:Real, 1}`: Cartesean inertial state. Returns position and velocity. [m; m/s]
+- x `x::AbstractArray{<:Real, 1}`: Cartesean inertial state. Returns position and velocity. [m; m/s]
 """
-function sOSCtoCART(x_oe::Array{<:Real, 1}; use_degrees::Bool=false, GM::Real=GM_EARTH)
+function sOSCtoCART(x_oe::AbstractArray{<:Real, 1}; use_degrees::Bool=false, GM::Real=GM_EARTH)
 
     if use_degrees == true
         # Copy and convert input from degrees to radians if necessary
@@ -238,14 +238,14 @@ The osculating elements are assumed to be (in order):
 6. _M_, Mean anomaly [rad]
 
 Arguments:
-- x `x::Array{<:Real, 1}`: Cartesean inertial state. Returns position and velocity. [m; m/s]
+- x `x::AbstractArray{<:Real, 1}`: Cartesean inertial state. Returns position and velocity. [m; m/s]
 - `use_degrees:Bool`: If `true` interpret input will be interpreted as being in degrees, and output will be returned in degrees.
 - `GM::Real`: Gravitational constant of central body. Defaults to `SatelliteDynamics.GM_EARTH` if none is provided.
 
 # Returns
-- x_oe `x::Array{<:Real, 1}`: Osculating orbital elements. See above for desription of the elements and their required order.
+- x_oe `x::AbstractArray{<:Real, 1}`: Osculating orbital elements. See above for desription of the elements and their required order.
 """
-function sCARTtoOSC(x::Array{<:Real, 1}; use_degrees::Bool=false, GM::Real=GM_EARTH)
+function sCARTtoOSC(x::AbstractArray{<:Real, 1}; use_degrees::Bool=false, GM::Real=GM_EARTH)
 
     # Initialize Cartesian Polistion and Velocity
     r = x[1:3]
diff --git a/src/attitude.jl b/src/attitude.jl
index 7c9e39f..80e4d4b 100644
--- a/src/attitude.jl
+++ b/src/attitude.jl
@@ -11,7 +11,7 @@ Arguments:
 - `use_degrees:Bool`: If `true` interpret input as being in degrees.
 
 Returns:
-- `r::Array{<:Real, 2}`: Rotation matrix
+- `r::AbstractArray{<:Real, 2}`: Rotation matrix
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
@@ -39,7 +39,7 @@ Arguments:
 - `use_degrees:Bool`: If `true` interpret input as being in degrees.
 
 Returns:
-- `r::Array{<:Real, 2}`: Rotation matrix
+- `r::AbstractArray{<:Real, 2}`: Rotation matrix
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
@@ -67,7 +67,7 @@ Arguments:
 - `use_degrees:Bool`: If `true` interpret input as being in degrees.
 
 Returns:
-- `r::Array{<:Real, 2}`: Rotation matrix
+- `r::AbstractArray{<:Real, 2}`: Rotation matrix
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.27.
@@ -156,16 +156,16 @@ representation.
 
 Data members:
 - `theta::Float64`: Angle of rotation
-- `vec::Array{Float64, 1}`: Axis of rotation
+- `vec::AbstractArray{Float64, 1}`: Axis of rotation
 
 References:
 1. J. Diebel, _Representing attitude: Euler angles, unit quaternions, and rotation vectors._ Matrix 58(15-16) (2006).
 """
 mutable struct EulerAxis
     angle::Float64
-    axis::Array{Float64, 1}
+    axis::AbstractArray{Float64, 1}
 
-    function EulerAxis(angle::Real, axis::Array{<:Real, 1})
+    function EulerAxis(angle::Real, axis::AbstractArray{<:Real, 1})
         if length(axis) != 3
             throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(axis)), Required size: (3,)"))
         end
@@ -179,7 +179,7 @@ end
 ##############
 
 # Quaternion Constructors 
-function Quaternion(vec::Array{<:Real, 1})
+function Quaternion(vec::AbstractArray{<:Real, 1})
     if length(vec) != 4
         throw(ArgumentError("Invalid array for Quaternion initialization. Input length: $(length(vec)), Required length: 4"))
     end
@@ -187,7 +187,7 @@ function Quaternion(vec::Array{<:Real, 1})
     Quaternion(vec...)
 end
 
-function Quaternion(mat::Array{<:Real, 2})
+function Quaternion(mat::AbstractArray{<:Real, 2})
     if size(mat) == (1,4)
         # Actually vector initialization. so it and return early
         return Quaternion(mat...)
@@ -377,7 +377,7 @@ Arguments:
 - `q::Quaternion`: Quaternion
 
 Returns:
-- `vec::Array{Float64, 1}`: Quaternion as a (4,) vector
+- `vec::AbstractArray{Float64, 1}`: Quaternion as a (4,) vector
 """
 function as_vector(q::Quaternion)
     return q[:]
@@ -392,7 +392,7 @@ Arguments:
 - `q::Quaternion`: Quaternion
 
 Returns:
-- `mat::Array{Float64, 2}`: Rotation Matrix on SO(3).
+- `mat::AbstractArray{Float64, 2}`: Rotation Matrix on SO(3).
 """
 function as_matrix(q::Quaternion)
     # initialize Empty Matrix
@@ -591,7 +591,7 @@ end
 # EulerAngle #
 ##############
 
-function EulerAngle(seq::Integer, vec::Array{<:Real, 1})
+function EulerAngle(seq::Integer, vec::AbstractArray{<:Real, 1})
     if length(vec) != 3
         throw(ArgumentError("Invalid array for EulerAngle initialization. Input length: $(length(vec)), Required length: 3"))
     end
@@ -599,7 +599,7 @@ function EulerAngle(seq::Integer, vec::Array{<:Real, 1})
     EulerAngle(seq, vec...)
 end
 
-function EulerAngle(seq::Integer, mat::Array{<:Real, 2})
+function EulerAngle(seq::Integer, mat::AbstractArray{<:Real, 2})
     if size(mat) != (3,3)
         throw(ArgumentError("Invalid array for Quaternion initialization. Input size: $(size(mat)), Required size: (3,3)"))
     end
@@ -714,7 +714,7 @@ Arguments:
 - `e::EulerAngle` Euler Angle
 
 Returns:
-- `evec::Array{Float64, 1}` Euler angles components in vector form.
+- `evec::AbstractArray{Float64, 1}` Euler angles components in vector form.
 """
 function as_vector(e::EulerAngle)
     return [e.phi, e.theta, e.psi]
@@ -741,7 +741,7 @@ function EulerAxis(angle::Real, v1::Real, v2::Real, v3::Real)
     return EulerAxis(angle, [v1, v2, v3])
 end
 
-function EulerAxis(vec::Array{<:Real, 1})
+function EulerAxis(vec::AbstractArray{<:Real, 1})
     if length(vec) != 4
         throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(vec)), Required size: (4,)"))
     end
@@ -749,7 +749,7 @@ function EulerAxis(vec::Array{<:Real, 1})
     return EulerAxis(vec[1], [vec[2], vec[3], vec[4]])
 end
 
-function EulerAxis(mat::Array{<:Real, 2})
+function EulerAxis(mat::AbstractArray{<:Real, 2})
     if size(mat) != (3,3)
         throw(ArgumentError("Invalid array for EulerAxis initialization. Input size: $(size(mat)), Required size: (3,3)"))
     end
diff --git a/src/coordinates.jl b/src/coordinates.jl
index 9444992..6db286f 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -14,13 +14,13 @@ export sGEOCtoECEF
 Convert geocentric position to equivalent Earth-fixed position.
 
 Arguments:
-- `geoc::Array{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
+- `geoc::AbstractArray{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
 - `use_degrees:Bool`: If `true` interpret input as being in degrees.
 
 Returns:
-- `ecef::Array{<:Real, 1}`: Earth-fixed coordinates [m]
+- `ecef::AbstractArray{<:Real, 1}`: Earth-fixed coordinates [m]
 """
-function sGEOCtoECEF(geoc::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sGEOCtoECEF(geoc::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Extract lat and lon
     lon = geoc[1]
     lat = geoc[2]
@@ -57,13 +57,13 @@ export sECEFtoGEOC
 Convert Earth-fixed position to geocentric location.
 
 Arguments:
-- `ecef::Array{<:Real, 1}`: Earth-fixed coordinated [m]
+- `ecef::AbstractArray{<:Real, 1}`: Earth-fixed coordinated [m]
 - `use_degrees:Bool`: If `true` returns result in units of degrees
 
 Returns:
 - `geoc`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
 """
-function sECEFtoGEOC(ecef::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sECEFtoGEOC(ecef::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Expand ECEF coordinates
     x, y, z = ecef
 
@@ -90,13 +90,13 @@ export sGEODtoECEF
 Convert geodetic position to equivalent Earth-fixed position.
 
 Arguments:
-- `geod::Array{<:Real, 1}`: Geodetic coordinates (lon, lat, altitude) [rad] / [deg]
+- `geod::AbstractArray{<:Real, 1}`: Geodetic coordinates (lon, lat, altitude) [rad] / [deg]
 - `use_degrees:Bool`: If `true` interpret input as being in degrees.
 
 Returns:
-- `ecef::Array{<:Real, 1}`: Earth-fixed coordinates [m]
+- `ecef::AbstractArray{<:Real, 1}`: Earth-fixed coordinates [m]
 """
-function sGEODtoECEF(geod::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sGEODtoECEF(geod::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Extract lat and lon
     lon = geod[1]
     lat = geod[2]
@@ -133,13 +133,13 @@ export sECEFtoGEOD
 Convert geodetic coordinaties to Earth-fixed position
 
 Arguments:
-- `ecef::Array{<:Real, 1}`: Earth-fixed position [m]
+- `ecef::AbstractArray{<:Real, 1}`: Earth-fixed position [m]
 - `use_degrees:Bool`: If `true` returns result in units of degrees
 
 Returns:
-- `geod::Array{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
+- `geod::AbstractArray{<:Real, 1}`: Geocentric coordinates (lon, lat, altitude) [rad] / [deg]
 """
-function sECEFtoGEOD(ecef::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sECEFtoGEOD(ecef::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Expand ECEF coordinates
     x, y, z = ecef
 
@@ -190,13 +190,13 @@ Compute the rotation matrix from the Earth-fixed to the East-North-Up
 coorindate basis.
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Real, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Real, 2}`: Topocentric rotation matrix
 """
-function rECEFtoENZ(ecef::Array{<:Real, 1} ; conversion::String="geodetic")   
+function rECEFtoENZ(ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")   
     if length(ecef) < 3
         throw(ArgumentError("Input coordinates must be length 3."))
     end
@@ -228,13 +228,13 @@ Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
 coorindate basis.
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function rENZtoECEF(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
+function rENZtoECEF(ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")
     # Check input coordinates
     if length(ecef) < 3
         throw(ArgumentError("Input coordinates must be length 3."))
@@ -249,14 +249,14 @@ Compute the coordinates of an object in the topocentric frame of an
 Earth-fixed frame
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
-- `ecef::Array{<:Real, 1}`: Coordinates of the object in Earth-fixed station
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `ecef::AbstractArray{<:Real, 1}`: Coordinates of the object in Earth-fixed station
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function sECEFtoENZ(station_ecef::Array{<:Real, 1}, ecef::Array{<:Real, 1} ; conversion::String="geodetic")
+function sECEFtoENZ(station_ecef::AbstractArray{<:Real, 1}, ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")
     # Check input sizes
     if length(ecef) < 3
         throw(ArgumentError("Input ecef state must be at least length 3."))
@@ -295,14 +295,14 @@ Compute the coordinates of an object in the topocentric frame of an
 Earth-fixed frame
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
-- `sez::Array{<:Real, 1}`: SEZ coordinates of the object
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `sez::AbstractArray{<:Real, 1}`: SEZ coordinates of the object
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function sENZtoECEF(station_ecef::Array{<:Real, 1}, enz::Array{<:Real, 1} ; conversion::String="geodetic")
+function sENZtoECEF(station_ecef::AbstractArray{<:Real, 1}, enz::AbstractArray{<:Real, 1} ; conversion::String="geodetic")
     # Check input sizes
     if length(enz) < 3
         throw(ArgumentError("Input ENZ state must be at least length 3."))
@@ -345,13 +345,13 @@ Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
 coorindate basis.
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function rECEFtoSEZ(ecef::Array{<:Real, 1} ; conversion::String="geodetic")    
+function rECEFtoSEZ(ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")    
     if length(ecef) < 3
         throw(ArgumentError("Input coordinates must be length 3."))
     end
@@ -384,13 +384,13 @@ Compute the rotation matrix from the Earth-fixed to the South-East-Zenith
 coorindate basis.
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function rSEZtoECEF(ecef::Array{<:Real, 1} ; conversion::String="geodetic")
+function rSEZtoECEF(ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")
     # Check input coordinates
     if length(ecef) < 3
         throw(ArgumentError("Input coordinates must be length 3."))
@@ -406,14 +406,14 @@ Compute the coordinates of an object in the topocentric frame of an
 Earth-fixed frame
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
-- `ecef::Array{<:Real, 1}`: Coordinates of the object in Earth-fixed station
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `ecef::AbstractArray{<:Real, 1}`: Coordinates of the object in Earth-fixed station
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function sECEFtoSEZ(station_ecef::Array{<:Real, 1}, ecef::Array{<:Real, 1} ; conversion::String="geodetic")    
+function sECEFtoSEZ(station_ecef::AbstractArray{<:Real, 1}, ecef::AbstractArray{<:Real, 1} ; conversion::String="geodetic")    
     # Check input sizes
     if length(ecef) < 3
         throw(ArgumentError("Input ecef state must be at least length 3."))
@@ -452,14 +452,14 @@ Compute the coordinates of an object in the topocentric frame of an
 Earth-fixed frame
 
 Arguments:
-- `station_ecef::Array{<:Real, 1}`: Earth-fixed cartesian station coordinates
-- `sez::Array{<:Real, 1}`: SEZ coordinates of the object
+- `station_ecef::AbstractArray{<:Real, 1}`: Earth-fixed cartesian station coordinates
+- `sez::AbstractArray{<:Real, 1}`: SEZ coordinates of the object
 - `conversion::Bool`: Conversion type to use. Can be "geocentric" or "geodetic"
 
 Returns:
-- `E::Array{Float64, 2}`: Topocentric rotation matrix
+- `E::AbstractArray{Float64, 2}`: Topocentric rotation matrix
 """
-function sSEZtoECEF(station_ecef::Array{<:Real, 1}, sez::Array{<:Real, 1} ; conversion::String="geodetic")
+function sSEZtoECEF(station_ecef::AbstractArray{<:Real, 1}, sez::AbstractArray{<:Real, 1} ; conversion::String="geodetic")
     # Check input sizes
     if length(sez) < 3
         throw(ArgumentError("Input SEZ state must be at least length 3."))
@@ -501,13 +501,13 @@ export sENZtoAZEL
 Convert East-North-Zenith topocentric state to azimuth, elevation, and range.
 
 Arguments:
-- `x::Array{<:Real, 1}`: East-North-Up state
+- `x::AbstractArray{<:Real, 1}`: East-North-Up state
 - `use_degrees:Bool`: If `true` returns result in units of degrees
 
 Returns:
-- `azel::Array{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
+- `azel::AbstractArray{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
 """
-function sENZtoAZEL(x::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sENZtoAZEL(x::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Check inputs
     if !(length(x) == 3 || length(x) == 6)
         throw(ArgumentError("Input ENZ state must be length 3 or 6."))
@@ -582,13 +582,13 @@ export sSEZtoAZEL
 Convert South-East-Zenith topocentric state to azimuth, elevation, and range.
 
 Arguments:
-- `x::Array{<:Real, 1}`: South-East-Zenith state
+- `x::AbstractArray{<:Real, 1}`: South-East-Zenith state
 - `use_degrees:Bool`: If `true` returns result in units of degrees
 
 Returns:
-- `azel::Array{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
+- `azel::AbstractArray{<:Real, 1}`: Azimuth, elevation and range [rad; rad; m]
 """
-function sSEZtoAZEL(x::Array{<:Real, 1} ; use_degrees::Bool=false)
+function sSEZtoAZEL(x::AbstractArray{<:Real, 1} ; use_degrees::Bool=false)
     # Check inputs
     if !(length(x) == 3 || length(x) == 6)
         throw(ArgumentError("Input rECEFtoSEZ state must be length 3 or 6."))
diff --git a/src/earth_environment/nrlmsise00.jl b/src/earth_environment/nrlmsise00.jl
index dc7561f..a6472a7 100644
--- a/src/earth_environment/nrlmsise00.jl
+++ b/src/earth_environment/nrlmsise00.jl
@@ -17,13 +17,13 @@ include("nrlmsise00_data.jl")
 Internal data structure used for setting model parameters
 """
 mutable struct NRLMSISE_Flags
-    switches::Array{<:Int, 1}
-    sw::Array{<:Real, 1}
-    swc::Array{<:Real, 1}
+    switches::AbstractArray{<:Int, 1}
+    sw::AbstractArray{<:Real, 1}
+    swc::AbstractArray{<:Real, 1}
 
-    function NRLMSISE_Flags(switches::Array{<:Int, 1}=zeros(Int, 24),
-            sw::Array{<:Real, 1}=zeros(Float32, 24),
-            swc::Array{<:Real, 1}=zeros(Float32, 24),)
+    function NRLMSISE_Flags(switches::AbstractArray{<:Int, 1}=zeros(Int, 24),
+            sw::AbstractArray{<:Real, 1}=zeros(Float32, 24),
+            swc::AbstractArray{<:Real, 1}=zeros(Float32, 24),)
         new(switches, sw, swc)
     end
 end
@@ -51,12 +51,12 @@ mutable struct NRLMSISE_Input
     f107A::Float64  # 81 day average of F10.7cm flux (centered on day)
     f107::Float64   # Daily F10.7cm flux for previous day
     ap::Float64     # Magnetic index (daily)
-    ap_array::Array{<:Real, 1} # Magnetic index array
+    ap_array::AbstractArray{<:Real, 1} # Magnetic index array
 
     function NRLMSISE_Input(year::Int=2000, doy::Int=1, sec::Float64=0.0, 
                 alt::Float64=0.0, g_lat::Float64=0.0, g_lon::Float64=0.0,
                 lst::Float64=0.0, f107A::Float64=0.0, f107::Float64=0.0, 
-                ap::Float64=0.0, ap_array::Array{<:Real, 1}=zeros(Float64, 7))
+                ap::Float64=0.0, ap_array::AbstractArray{<:Real, 1}=zeros(Float64, 7))
         new(year, doy, sec, alt, g_lat, g_lon, lst, f107A, f107, ap, ap_array)
     end
 end
@@ -67,10 +67,10 @@ d = zeros(Float64, 9)
 t = zeros(Float64, 2)
 """
 mutable struct NRLMSISE_Output
-    d::Array{Float64, 1}
-    t::Array{Float64, 1}
+    d::AbstractArray{Float64, 1}
+    t::AbstractArray{Float64, 1}
 
-    function NRLMSISE_Output(d::Array{<:Real, 1}=zeros(Float64, 9), t::Array{<:Real, 1}=zeros(Float64, 2))
+    function NRLMSISE_Output(d::AbstractArray{<:Real, 1}=zeros(Float64, 9), t::AbstractArray{<:Real, 1}=zeros(Float64, 2))
         new(d, t)
     end
 end
@@ -222,15 +222,15 @@ end
 Integrate cubic spline from xa[1] to x
 
 Arguments:
--`xa::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
--`ya::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
--`y2a::Array{Real, 1}` Array of second derivatives
+-`xa::AbstractArray{Real, 1}` Array of tabulated function inputs in ascending order by x
+-`ya::AbstractArray{Real, 1}` Array of tabulated function outputs in ascending order by x
+-`y2a::AbstractArray{Real, 1}` Array of second derivatives
 -`x::Real` Array output point
 
 Returns:
 -`y::Real` Array output value
 """
-function splini(xa::Array{<:Real, 1}, ya::Array{<:Real, 1}, y2a::Array{<:Real, 1}, n::Int, x::Real)
+function splini(xa::AbstractArray{<:Real, 1}, ya::AbstractArray{<:Real, 1}, y2a::AbstractArray{<:Real, 1}, n::Int, x::Real)
     yi  = 0
     klo = 1
     khi = 2
@@ -261,15 +261,15 @@ end
 Interpolate cubic spline from to x
 
 Arguments:
--`xa::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
--`ya::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
--`y2a::Array{Real, 1}` Array of second derivates
+-`xa::AbstractArray{Real, 1}` Array of tabulated function inputs in ascending order by x
+-`ya::AbstractArray{Real, 1}` Array of tabulated function outputs in ascending order by x
+-`y2a::AbstractArray{Real, 1}` Array of second derivates
 -`x::Real` Array output point
 
 Returns:
 -`y::Real` Array output value
 """
-function splint(xa::Array{<:Real, 1}, ya::Array{<:Real, 1}, y2a::Array{<:Real, 1}, n::Int, x::Real)
+function splint(xa::AbstractArray{<:Real, 1}, ya::AbstractArray{<:Real, 1}, y2a::AbstractArray{<:Real, 1}, n::Int, x::Real)
     klo = 1
     khi = n
     k   = 0
@@ -300,16 +300,16 @@ end
 Interpolate cubic spline from to x
 
 Arguments:
--`x::Array{Real, 1}` Array of tabulated function inputs in ascending order by x
--`y::Array{Real, 1}` Array of tabulated function outputs in ascending order by x
+-`x::AbstractArray{Real, 1}` Array of tabulated function inputs in ascending order by x
+-`y::AbstractArray{Real, 1}` Array of tabulated function outputs in ascending order by x
 -`yp1::Real` Specified derivative at x[1]   (yp1 >= 1e30 signals second derivative is zero)
 -`ypn::Real` Specified derivative at x[end] (ypn >= 1e30 signals second derivative is zero)
 -`x::Real` Array output point
 
 Returns:
--`y::Array{Float64, 1}` Array second derivatives
+-`y::AbstractArray{Float64, 1}` Array second derivatives
 """
-function spline(x::Array{<:Real, 1}, y::Array{<:Real, 1}, n::Int, yp1::Real, ypn::Real)
+function spline(x::AbstractArray{<:Real, 1}, y::AbstractArray{<:Real, 1}, n::Int, yp1::Real, ypn::Real)
     u  = zeros(Float64, n) # Need the + 1 here because C implementation already allocates extra
     y2 = zeros(Float64, n)
 
@@ -382,8 +382,8 @@ Returns:
 - `meso_tgn2`
 """
 function densm!(alt::Real, d0::Real, xm::Real, tz::Real, 
-                mn3::Int, zn3::Array{<:Real, 1}, tn3::Array{<:Real, 1}, tgn3::Array{<:Real, 1}, 
-                mn2::Int, zn2::Array{<:Real, 1}, tn2::Array{<:Real, 1}, tgn2::Array{<:Real, 1};
+                mn3::Int, zn3::AbstractArray{<:Real, 1}, tn3::AbstractArray{<:Real, 1}, tgn3::AbstractArray{<:Real, 1}, 
+                mn2::Int, zn2::AbstractArray{<:Real, 1}, tn2::AbstractArray{<:Real, 1}, tgn2::AbstractArray{<:Real, 1};
                 gsurf::Real, re::Real)
     
 
@@ -530,7 +530,7 @@ end
 
 function densu!(alt::Real, dlb::Real, tinf::Real, tlb::Real, xm::Real, 
                alpha::Real, tz::Real, zlb::Real, s2::Real, 
-               mn1::Int, zn1::Array{<:Real, 1}, tn1::Array{<:Real, 1}, tgn1::Array{<:Real, 1};
+               mn1::Int, zn1::AbstractArray{<:Real, 1}, tn1::AbstractArray{<:Real, 1}, tgn1::AbstractArray{<:Real, 1};
                gsurf::Real, re::Real)
 
     x          = 0
@@ -693,7 +693,7 @@ function densu!(alt::Real, dlb::Real, tinf::Real, tlb::Real, xm::Real,
 end
 
 # Equation A24a
-@inline function g0(a::Real, p::Array{<:Real, 1})
+@inline function g0(a::Real, p::AbstractArray{<:Real, 1})
     return (a - 4.0 + (p[26] - 1.0) * (a - 4.0 + (exp(-sqrt(p[25]*p[25]) * (a - 4.0)) - 1.0) / sqrt(p[25]*p[25])))
 end
 
@@ -703,7 +703,7 @@ end
 end
 
 # Equation A24a
-@inline function sg0(ex::Real, p::Array{<:Real, 1}, ap::Array{<:Real, 1})
+@inline function sg0(ex::Real, p::AbstractArray{<:Real, 1}, ap::AbstractArray{<:Real, 1})
     return (g0(ap[2],p) + (g0(ap[3],p)*ex + g0(ap[4],p)*ex*ex + 
             g0(ap[5],p)*ex^3.0	+ (g0(ap[6],p)*ex^4.0 + 
             g0(ap[7],p)*ex^12.0)*(1.0- ex^8.0)/(1.0-ex)))/sumex(ex)
@@ -716,7 +716,7 @@ end
 """
 Calculate G(L) function
 """
-function globe7(p::Array{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags)
+function globe7(p::AbstractArray{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags)
     # Working variables
     t = zeros(Float64, 15)
 
@@ -940,8 +940,8 @@ end
 """
 Calculate G(L) function for lower atmosphere
 """
-function glob7s(p::Array{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags; 
-                dfa::Real, plg::Array{<:Array{<:Real,1},1}, ctloc::Real, stloc::Real, c2tloc::Real, s2tloc::Real, s3tloc::Real, c3tloc::Real)
+function glob7s(p::AbstractArray{<:Real, 1}, input::NRLMSISE_Input, flags::NRLMSISE_Flags; 
+                dfa::Real, plg::AbstractArray{<:Array{<:Real,1},1}, ctloc::Real, stloc::Real, c2tloc::Real, s2tloc::Real, s3tloc::Real, c3tloc::Real)
     # Working variables
     pset = 2.0
     t    = zeros(Float64, 14)
@@ -1593,7 +1593,7 @@ Computes the local atmospheric density using the NRLMSISE00 atmosphere model.
 
 Arguments:
 - `epc::Epoch`: Epoch of computation. Used to lookup space weather data
-- `x::Array{<:Real, 1}`: Satellite state in geodetic coordinates [lon, lat, alt]
+- `x::AbstractArray{<:Real, 1}`: Satellite state in geodetic coordinates [lon, lat, alt]
 - `use_degrees:Bool`: If `true` interprets geodetic inputs as being in degrees
 
 Returns:
@@ -1605,7 +1605,7 @@ Notes:
 References:
 1. _Picone, JM, et al._ NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues _Journal of Geophysical Research: Space Physics_
 """
-function density_nrlmsise00(epc::Epoch, x::Array{<:Real, 1}; use_degrees::Bool=false)
+function density_nrlmsise00(epc::Epoch, x::AbstractArray{<:Real, 1}; use_degrees::Bool=false)
     # Create NRLMSISE00 model input
     input  = NRLMSISE_Input()
     flags  = NRLMSISE_Flags()
diff --git a/src/orbit_dynamics.jl b/src/orbit_dynamics.jl
index ae5d988..18d5ab5 100644
--- a/src/orbit_dynamics.jl
+++ b/src/orbit_dynamics.jl
@@ -52,15 +52,15 @@ of the central body. Returns the acceleration vector of the satellite.
 Assumes the satellite is much, much less massive than the central body.
 
 Arguments:
-- `r_sat::Array{<:Real, 1}`: satellite position in a commonn inertial frame [m]
-- `r_body::Array{<:Real, 1}`: position of body in a commonn inertial frame [m]
-- `GM::Array{<:Real, 1}`: gravitational coeffient of attracting body [m^3/s^2] Default: SatelliteDynamics.GM_EARTH)
+- `r_sat::AbstractArray{<:Real, 1}`: satellite position in a commonn inertial frame [m]
+- `r_body::AbstractArray{<:Real, 1}`: position of body in a commonn inertial frame [m]
+- `GM::AbstractArray{<:Real, 1}`: gravitational coeffient of attracting body [m^3/s^2] Default: SatelliteDynamics.GM_EARTH)
 (Default: SatelliteDynamics.GM_EARTH
 
 Return:
-- `a::Array{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
 """
-function accel_point_mass(r_sat::Array{<:Real, 1}, r_body::Array{<:Real, 1}, gm_body::Real=GM_EARTH)
+function accel_point_mass(r_sat::AbstractArray{<:Real, 1}, r_body::AbstractArray{<:Real, 1}, gm_body::Real=GM_EARTH)
     # Restrict inputs to position only
     r_sat  = r_sat[1:3]
     r_body = r_body[1:3]
@@ -79,14 +79,14 @@ Computes the acceleration on a satellite caused by a point-mass approximation
 of a massive body. Returns the acceleration vector of the satellite.
 
 Arguments:
-- `r_sat::Array{<:Real, 1}`: satellite position in the inertial frame [m]
-- `GM::Array{<:Real, 1}`: gravitational coeffient of attracting body [m^3/s^2] Default: SatelliteDynamics.GM_EARTH)
+- `r_sat::AbstractArray{<:Real, 1}`: satellite position in the inertial frame [m]
+- `GM::AbstractArray{<:Real, 1}`: gravitational coeffient of attracting body [m^3/s^2] Default: SatelliteDynamics.GM_EARTH)
 (Default: SatelliteDynamics.GM_EARTH
 
 Return:
-- `a::Array{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
 """
-function accel_point_mass(x::Array{<:Real, 1}, gm_body::Real=GM_EARTH)
+function accel_point_mass(x::AbstractArray{<:Real, 1}, gm_body::Real=GM_EARTH)
     # Restrict inputs to position only. Considered in body frame
     r  = x[1:3]
 
@@ -101,8 +101,8 @@ Compute the gravitational acceleration at a model given a spherical harmonic
 gravity field model.
 
 Arguments:
-- `r::Array{<:Real, 1}`: Position of the point in the body (field) fixed frame. [m]
-- `coef::Array{<:Real, 2}`: Gravity coefficients stored in dense matrix form. C_nm terns are stored along rows indexed from C_00 in coef[1, 1] to C_nm in coef[n+1, m+1]. S_nm terms are stored along matrix columns with S_n0 stored in coef[0, n+1]
+- `r::AbstractArray{<:Real, 1}`: Position of the point in the body (field) fixed frame. [m]
+- `coef::AbstractArray{<:Real, 2}`: Gravity coefficients stored in dense matrix form. C_nm terns are stored along rows indexed from C_00 in coef[1, 1] to C_nm in coef[n+1, m+1]. S_nm terms are stored along matrix columns with S_n0 stored in coef[0, n+1]
 - `n_max::Integer`: Maximum degree coefficient to use in expansion
 - `m_max::Integer`: Maximum order coefficient to use in the expansion. Must be less than the degree.
 - `r_ref::Real`: Reference distance of the gravity field.
@@ -110,9 +110,9 @@ Arguments:
 - `normralized::Bool`: Whether the input gravity field coefficients are normalized coefficients (Default: true)
 
 Returns:
-- `a::Array{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
 """
-function spherical_harmonic_gravity(r::Array{<:Real, 1}, coef::Array{<:Real, 2}, n_max::Integer, m_max::Integer, r_ref::Real, GM::Real; normalized::Bool=true)
+function spherical_harmonic_gravity(r::AbstractArray{<:Real, 1}, coef::AbstractArray{<:Real, 2}, n_max::Integer, m_max::Integer, r_ref::Real, GM::Real; normalized::Bool=true)
     # Intermediate computations
     r_sqr = dot(r, r)
     rho   = r_ref^2/r_sqr
@@ -206,18 +206,18 @@ Computes the accleration caused by Earth gravity as modeled by a spherical
 harmonic gravity field.
 
 Arguments:
-- `r_sat::Array{<:Real, 1}`: Satellite position in the inertial frame [m]
-- `R_eci_ecef::Array{<:Real, 2}`: Rotation matrix transforming a vector from the inertial to body-fixed reference frames. 
+- `r_sat::AbstractArray{<:Real, 1}`: Satellite position in the inertial frame [m]
+- `R_eci_ecef::AbstractArray{<:Real, 2}`: Rotation matrix transforming a vector from the inertial to body-fixed reference frames. 
 - `n_max::Integer`: Maximum degree coefficient to use in expansion
 - `m_max::Integer`: Maximum order coefficient to use in the expansion. Must be less than the degree.
     
 Return:
-- `a::Array{<:Real, 1}`: Gravitational acceleration in X, Y, and Z inertial directions [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Gravitational acceleration in X, Y, and Z inertial directions [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.56-68.
 """
-function accel_gravity(x::Array{<:Real, 1}, R_eci_ecef::Array{<:Real, 2}, n_max::Int=20, m_max::Int=20)
+function accel_gravity(x::AbstractArray{<:Real, 1}, R_eci_ecef::AbstractArray{<:Real, 2}, n_max::Int=20, m_max::Int=20)
     
     # Check Limits of Gravity Field
     if n_max > GRAVITY_MODEL.n_max
@@ -258,7 +258,7 @@ Argument:
 - `epc::Epoch`: Epoch
 
 Returns:
-- `r_sun::Array{<:Real, 1}`: Position vector of the Sun in the Earth-centered inertial fame.
+- `r_sun::AbstractArray{<:Real, 1}`: Position vector of the Sun in the Earth-centered inertial fame.
 
 Notes:
 1. The EME2000 inertial frame is for most purposes equivalent to the GCRF frame.
@@ -294,7 +294,7 @@ Argument:
 - `epc::Epoch`: Epoch
 
 Returns:
-- `r_moon::Array{<:Real, 1}`: Position vector of the Moon in the Earth-centered inertial fame.
+- `r_moon::AbstractArray{<:Real, 1}`: Position vector of the Moon in the Earth-centered inertial fame.
 
 Notes:
 1. The EME2000 inertial frame is for most purposes equivalent to the GCRF frame.
@@ -352,23 +352,23 @@ Computes the acceleration of a satellite in the inertial frame due to the
 gravitational attraction of the Sun.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
-- `r_sun::Array{<:Real, 1}`: Position of sun in inertial frame.
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `r_sun::AbstractArray{<:Real, 1}`: Position of sun in inertial frame.
 
 Return:
-- `a::Array{<:Real, 1}`: Acceleration due to the Sun's gravity in the inertial frame [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration due to the Sun's gravity in the inertial frame [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.69-70.
 """
-function accel_thirdbody_sun(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
+function accel_thirdbody_sun(x::AbstractArray{<:Real, 1}, r_sun::AbstractArray{<:Real, 1})
     # Acceleration due to sun point mass
     a_sun = accel_point_mass(x[1:3], r_sun, GM_SUN)
 
     return a_sun
 end
 
-function accel_thirdbody_sun(epc::Epoch, x::Array{<:Real, 1})
+function accel_thirdbody_sun(epc::Epoch, x::AbstractArray{<:Real, 1})
     # Compute solar position
     r_sun = sun_position(epc)
 
@@ -384,23 +384,23 @@ Computes the acceleration of a satellite in the inertial frame due to the
 gravitational attraction of the Moon.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
-- `r_moon::Array{<:Real, 1}`: Position of moon in inertial frame.
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `r_moon::AbstractArray{<:Real, 1}`: Position of moon in inertial frame.
 
 Returns:
-- `a::Array{<:Real, 1}`: Acceleration due to the Moon's gravity in the inertial frame [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration due to the Moon's gravity in the inertial frame [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.69-70.
 """
-function accel_thirdbody_moon(x::Array{<:Real, 1}, r_moon::Array{<:Real, 1})
+function accel_thirdbody_moon(x::AbstractArray{<:Real, 1}, r_moon::AbstractArray{<:Real, 1})
     # Acceleration due to moon point mass
     a_moon = accel_point_mass(x[1:3], r_moon, GM_MOON)
 
     return a_moon
 end
 
-function accel_thirdbody_moon(epc::Epoch, x::Array{<:Real, 1})
+function accel_thirdbody_moon(epc::Epoch, x::AbstractArray{<:Real, 1})
     # Compute solar position
     r_moon = moon_position(epc)
 
@@ -419,8 +419,8 @@ export density_harris_priester
 Computes the local density using the Harris-Priester density model.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
-- `r_sun::Array{<:Real, 1}`: Position of sun in inertial frame.
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `r_sun::AbstractArray{<:Real, 1}`: Position of sun in inertial frame.
 
 Returns:
 - `rho:Float64`: Local atmospheric density [kg/m^3]
@@ -428,7 +428,7 @@ Returns:
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.89-91.
 """
-function density_harris_priester(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
+function density_harris_priester(x::AbstractArray{<:Real, 1}, r_sun::AbstractArray{<:Real, 1})
     # Harris-Priester Constants
     hp_upper_limit =   1000.0          # Upper height limit [km]
     hp_lower_limit =    100.0          # Lower height limit [km]
@@ -518,7 +518,7 @@ function density_harris_priester(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
     return density
 end
 
-function density_harris_priester(epc::Epoch, x::Array{<:Real, 1})
+function density_harris_priester(epc::Epoch, x::AbstractArray{<:Real, 1})
     r_sun = sun_position(epc)
     return density_harris_priester(x, r_sun)
 end
@@ -530,20 +530,20 @@ drag assuming that the ballistic properties of the spacecraft are captured by
 the coefficient of drag.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
 - `rho::Real`: atmospheric density [kg/m^3]
 - `mass::Real`: Spacecraft mass [kg]
 - `area::Real`: Wind-facing cross-sectional area [m^2]
 - `Cd::Real`: coefficient of drag [dimensionless]
-- `T::Array{<:Real, 2}`: Rotation matrix from the inertial to the true-of-date frame
+- `T::AbstractArray{<:Real, 2}`: Rotation matrix from the inertial to the true-of-date frame
 
 Return:
-- `a::Array{<:Real, 1}`: Acceleration due to drag in the X, Y, and Z inertial directions. [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Acceleration due to drag in the X, Y, and Z inertial directions. [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.83-86.
 """
-function accel_drag(x::Array{<:Real, 1}, rho::Real, mass::Real, area::Real, Cd::Real, T::Array{<:Real, 2})
+function accel_drag(x::AbstractArray{<:Real, 1}, rho::Real, mass::Real, area::Real, Cd::Real, T::AbstractArray{<:Real, 2})
 
     # Constants
     omega = [0, 0, OMEGA_EARTH]
@@ -573,8 +573,8 @@ Computes the illumination fraction of a satellite in Earth orbit using a
 cylindrical Earth shadow model.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
-- `r_sun::Array{<:Real, 1}`: Position of sun in inertial frame.
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `r_sun::AbstractArray{<:Real, 1}`: Position of sun in inertial frame.
 
 Return:
 - `nu::Float64`: Illumination fraction (0 <= nu <= 1). nu = 0 means spacecraft in complete shadow, nu = 1 mean spacecraft fully illuminated by sun.
@@ -582,7 +582,7 @@ Return:
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.80-83.
 """
-function eclipse_cylindrical(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
+function eclipse_cylindrical(x::AbstractArray{<:Real, 1}, r_sun::AbstractArray{<:Real, 1})
     # Satellite inertial position
     r = x[1:3]
     
@@ -601,7 +601,7 @@ function eclipse_cylindrical(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
     return nu
 end
 
-function eclipse_cylindrical(epc::Epoch, x::Array{<:Real, 1})
+function eclipse_cylindrical(epc::Epoch, x::AbstractArray{<:Real, 1})
     r_sun = sun_position(epc)
     return eclipse_cylindrical(x, r_sun)
 end
@@ -612,8 +612,8 @@ Computes the illumination fraction of a satellite in Earth orbit using a
 conical Earth shadow model.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
-- `r_sun::Array{<:Real, 1}`: Position of sun in inertial frame.
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `r_sun::AbstractArray{<:Real, 1}`: Position of sun in inertial frame.
 
 Return:
 - `nu::Float64`: Illumination fraction (0 <= nu <= 1). nu = 0 means spacecraft in complete shadow, nu = 1 mean spacecraft fully illuminated by sun.
@@ -621,7 +621,7 @@ Return:
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.80-83.
 """
-function eclipse_conical(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
+function eclipse_conical(x::AbstractArray{<:Real, 1}, r_sun::AbstractArray{<:Real, 1})
     # Satellite inertial position
     r = x[1:3]
 
@@ -654,7 +654,7 @@ function eclipse_conical(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1})
 end
 
 
-function eclipse_conical(epc::Epoch, x::Array{<:Real, 1})
+function eclipse_conical(epc::Epoch, x::AbstractArray{<:Real, 1})
     r_sun = sun_position(epc)
     return eclipse_conical(x, r_sun)
 end
@@ -665,15 +665,15 @@ pressure assuming the reflecting surface is a flat plate pointed directly at
 the Sun.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
 
 Returns:
-- `a::Array{<:Real, 1}`: Satellite acceleration due to solar radiation pressure [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Satellite acceleration due to solar radiation pressure [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.77-79.
 """
-function accel_srp(x::Array{<:Real, 1}, r_sun::Array{<:Real, 1}, mass::Real=0, area::Real=0, CR::Real=1.8, p0::Real=P_SUN, au::Real=AU)
+function accel_srp(x::AbstractArray{<:Real, 1}, r_sun::AbstractArray{<:Real, 1}, mass::Real=0, area::Real=0, CR::Real=1.8, p0::Real=P_SUN, au::Real=AU)
     # Spacecraft position vector
     r = x[1:3]
 
@@ -696,15 +696,15 @@ export accel_relativity
 due to the combined effects of special and general relativity.
 
 Arguments:
-- `x::Array{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
+- `x::AbstractArray{<:Real, 1}`: Satellite Cartesean state in the inertial reference frame [m; m/s]
 
 Returns:
-- `a::Array{<:Real, 1}`: Satellite acceleration due to relativity. [m/s^2]
+- `a::AbstractArray{<:Real, 1}`: Satellite acceleration due to relativity. [m/s^2]
 
 References:
 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.110-112.
 """
-function accel_relativity(x::Array{<:Real, 1})
+function accel_relativity(x::AbstractArray{<:Real, 1})
     # Extract state variables
     r = x[1:3]
     v = x[4:6]
diff --git a/src/sgp_models.jl b/src/sgp_models.jl
index e68d2cb..09e1e7a 100644
--- a/src/sgp_models.jl
+++ b/src/sgp_models.jl
@@ -2414,7 +2414,7 @@ Arguments:
 - `opscode::Char` Operatitonal mode of propagators
 
 Returns:
-- `rv::Array{Float64, 1}` Position and velocity as output by TLE propagator. 
+- `rv::AbstractArray{Float64, 1}` Position and velocity as output by TLE propagator. 
     Units [m m/s]
 """
 function state(tle::TLE, epc::Epoch)
diff --git a/src/simulation/integrators.jl b/src/simulation/integrators.jl
index 35a031f..29d1dac 100644
--- a/src/simulation/integrators.jl
+++ b/src/simulation/integrators.jl
@@ -28,7 +28,7 @@ Internal funciton used to compute the variational matrix.
 Arguments:
 - `rk4::RK4` 4-th order Runge-Kutta numerical integrator object
 - `epc::Union{Real, Epoch}` Absolute time of start of integration step
-- `x::Array{<:Real, 1}` State vector for linearizing perturbations around
+- `x::AbstractArray{<:Real, 1}` State vector for linearizing perturbations around
 - `apert::Union{Real, Array{<:Real,1}}` Absolute perturbations to state elements
 used to numerically calculate the partial derivatives through a forward difference
 method.
@@ -38,7 +38,7 @@ Returns:
 matrix of partial derivatives of the dynamics function _f_ with respect to the 
 input state _x_.
 """
-function varmat(rk4::RK4, epc::Union{Real,Epoch}, x::Array{<:Real, 1}; apert=1.0::Union{Real, Array{<:Real,1}})
+function varmat(rk4::RK4, epc::Union{Real,Epoch}, x::AbstractArray{<:Real, 1}; apert=1.0::Union{Real, Array{<:Real,1}})
     # Get state length
     n = length(x)
 
@@ -73,7 +73,7 @@ function varmat(rk4::RK4, epc::Union{Real,Epoch}, x::Array{<:Real, 1}; apert=1.0
     return V
 end
 
-function istep(rk4::RK4, epc::Union{Real,Epoch}, dt::Real, x::Array{<:Real, 1})
+function istep(rk4::RK4, epc::Union{Real,Epoch}, dt::Real, x::AbstractArray{<:Real, 1})
     # Compute State coefficients
 
     # Perform internal steps
@@ -95,10 +95,10 @@ Arguments:
 - `rk4::RK4` 4-th order Runge-Kutta numerical integrator object
 - `epc::Union{Real, Epoch}` Absolute time of start of integration step
 - `dt::Real` Integration step size
-- `x::Array{<:Real, 1}` State vector
-- `phi::Array{<:Real, 2}` State transition matrix at start of integration step
+- `x::AbstractArray{<:Real, 1}` State vector
+- `phi::AbstractArray{<:Real, 2}` State transition matrix at start of integration step
 """
-function istep(rk4::RK4, epc::Union{Real,Epoch}, dt::Real, x::Array{<:Real, 1}, phi::Array{<:Real, 2}; apert=1.0::Union{Real, Array{<:Real,1}})
+function istep(rk4::RK4, epc::Union{Real,Epoch}, dt::Real, x::AbstractArray{<:Real, 1}, phi::AbstractArray{<:Real, 2}; apert=1.0::Union{Real, Array{<:Real,1}})
     # Compute State coefficients
 
     # Perform state integration
diff --git a/src/simulation/propagators.jl b/src/simulation/propagators.jl
index 2c839c2..02e8bb0 100644
--- a/src/simulation/propagators.jl
+++ b/src/simulation/propagators.jl
@@ -8,7 +8,7 @@ Compute the state derivative.
 
 Arguments:
 - `epc::Epoch`: Current epoch
-- `x::Array{<:Real, 1}`: Satellite state vector
+- `x::AbstractArray{<:Real, 1}`: Satellite state vector
 - `mass::Real`: Satellite mass [kg]
 - `area_drag`: Velocity-facing area affected by drag. [m^2]
 - `coef_drag`: Coefficient of drag [dimensionless]
@@ -23,9 +23,9 @@ Arguments:
 - `relativity::Bool`: Include relativistic effects in force model (Default: `true`)
 
 Returns:
-- `dx::Array{<:Float64, 1}`: Satellite state derivative, velocity and accelerations [m; m/s]
+- `dx::AbstractArray{<:Float64, 1}`: Satellite state derivative, velocity and accelerations [m; m/s]
 """
-function fderiv_earth_orbit(epc::Epoch, x::Array{<:Real} ;
+function fderiv_earth_orbit(epc::Epoch, x::AbstractArray{<:Real} ;
              mass::Real=1.0, area_drag::Real=1.0, coef_drag::Real=2.3, 
              area_srp::Real=1.0, coef_srp::Real=1.8, 
              n_grav::Integer=20, m_grav::Integer=20, 
@@ -102,7 +102,7 @@ Attributes:
 the state vector is requested to propagate to a time smaller than this step size,
 which it will do.
 - `epc::Epoch` Epoch of state
-- `x::Array{Float64, 1}` State vector. Earth-centered inertial Cartesian state.
+- `x::AbstractArray{Float64, 1}` State vector. Earth-centered inertial Cartesian state.
 - `phi::Union{Nothing, Array{Float64, 2}}` State transition matrix, or the matrix
 of partial derivatives of the state at the current time with respect to the 
 start of propagation.
@@ -128,11 +128,11 @@ mutable struct EarthInertialState
     rk4::RK4
     dt::Real
     epc::Epoch
-    x::Array{Float64, 1}
+    x::AbstractArray{Float64, 1}
     phi::Union{Nothing, Array{Float64, 2}}
 end
 
-function EarthInertialState(epc::Epoch, x::Array{<:Real, 1}, phi::Union{Nothing, Array{Float64, 2}}=nothing; dt::Real=30.0, kwargs...)
+function EarthInertialState(epc::Epoch, x::AbstractArray{<:Real, 1}, phi::Union{Nothing, Array{Float64, 2}}=nothing; dt::Real=30.0, kwargs...)
     rk4 = RK4(fderiv_earth_orbit; kwargs...)
     return EarthInertialState(rk4, dt, epc, x, phi)
 end
@@ -206,10 +206,10 @@ Arguments:
 a real number to advance the state by or the Epoch
 
 Returns:
-- `t::Array{Float64, 1}` Elapsed time as a scalar from the initial simulation epoch 
-- `epc::Array{Epoch, 1}` Epoch at each timestep 
-- `x::Array{Float64, 2}` State vectors at each time step. Time is along second axis
-- `Phi::Array{Float64, 2}` Stacked array of state transition matrices
+- `t::AbstractArray{Float64, 1}` Elapsed time as a scalar from the initial simulation epoch 
+- `epc::AbstractArray{Epoch, 1}` Epoch at each timestep 
+- `x::AbstractArray{Float64, 2}` State vectors at each time step. Time is along second axis
+- `Phi::AbstractArray{Float64, 2}` Stacked array of state transition matrices
 """
 function sim!(state::EarthInertialState, time::Union{Real, Epoch}=0.0)
     if typeof(time) <: Real
diff --git a/src/universe.jl b/src/universe.jl
index db2c966..ae63101 100644
--- a/src/universe.jl
+++ b/src/universe.jl
@@ -344,7 +344,7 @@ struct GravModel
     GM::Float64
     n_max::Int64
     m_max::Int64
-    data::Array{Float64, 2}
+    data::AbstractArray{Float64, 2}
 end
 
 function GravModel(filepath::String)
diff --git a/test/simulation/test_integrators.jl b/test/simulation/test_integrators.jl
index 618e05a..be4aeec 100644
--- a/test/simulation/test_integrators.jl
+++ b/test/simulation/test_integrators.jl
@@ -12,7 +12,7 @@ end
 """
 Simple test function for spherical earth dynamics
 """
-function point_earth(epc::Epoch, x::Array{<:Real, 1})
+function point_earth(epc::Epoch, x::AbstractArray{<:Real, 1})
     # Restrict inputs to position only. Considered in body frame
     r  = x[1:3]
     v  = x[4:6]