diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..0535e07
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,15 @@
+# Change Log
+All notable changes to this project will be documented in this file.
+This project adheres to [Semantic Versioning](http://semver.org/).
+
+## [0.?.?] - 2020-07-03
+- Added coords from pysat.utils
+
+## [0.0.3] - 2020-06-15
+- pypi compatibility
+
+## [0.0.2] - 2020-05-13
+- zenodo link
+
+## [0.0.1] - 2020-05-13
+- Alpha release
diff --git a/pysatMadrigal/__init__.py b/pysatMadrigal/__init__.py
index 0b73a26..c1a628f 100644
--- a/pysatMadrigal/__init__.py
+++ b/pysatMadrigal/__init__.py
@@ -1 +1,2 @@
 from pysatMadrigal import instruments  # noqa F401
+from pysatMadrigal import utils  # noqa F401
diff --git a/pysatMadrigal/instruments/jro_isr.py b/pysatMadrigal/instruments/jro_isr.py
index 149be12..3f1ec34 100644
--- a/pysatMadrigal/instruments/jro_isr.py
+++ b/pysatMadrigal/instruments/jro_isr.py
@@ -38,13 +38,14 @@
 from __future__ import absolute_import
 import datetime as dt
 import functools
+import logging
 import numpy as np
 
-from pysatMadrigal.instruments.methods import madrigal as mad_meth
 from pysat.instruments.methods import general as mm_gen
-from pysat.utils import coords
 
-import logging
+from pysatMadrigal.instruments.methods import madrigal as mad_meth
+from pysatMadrigal.utils import coords
+
 logger = logging.getLogger(__name__)
 
 
diff --git a/pysatMadrigal/tests/test_utils_coords.py b/pysatMadrigal/tests/test_utils_coords.py
new file mode 100644
index 0000000..4b2549e
--- /dev/null
+++ b/pysatMadrigal/tests/test_utils_coords.py
@@ -0,0 +1,304 @@
+"""
+tests the pysat coords area
+"""
+import numpy as np
+
+import pytest
+
+from pysatMadrigal.utils import coords
+
+
+class TestGeodeticGeocentric():
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        self.val = {'lat': 45.0, 'lon': 8.0, 'az': 52.0, 'el': 63.0}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+    @pytest.mark.parametrize("kwargs,target",
+                             [({}, [44.8075768, 8.0, 6367.48954386]),
+                              ({'inverse': False},
+                               [44.8075768, 8.0, 6367.48954386]),
+                              ({'inverse': True},
+                               [45.1924232, 8.0, 6367.3459085])])
+    def test_geodetic_to_geocentric(self, kwargs, target):
+        """Test conversion from geodetic to geocentric coordinates"""
+
+        self.out = coords.geodetic_to_geocentric(self.val['lat'],
+                                                 lon_in=self.val['lon'],
+                                                 **kwargs)
+
+        for i, self.loc in enumerate(self.out):
+            assert np.all(abs(self.loc - target[i]) < 1.0e-6)
+            if isinstance(self.loc, np.ndarray):
+                assert self.loc.shape == self.val['lat'].shape
+
+    def test_geodetic_to_geocentric_and_back(self):
+        """Tests the reversibility of geodetic to geocentric conversions"""
+
+        self.out = coords.geodetic_to_geocentric(self.val['lat'],
+                                                 lon_in=self.val['lon'],
+                                                 inverse=False)
+        self.loc = coords.geodetic_to_geocentric(self.out[0],
+                                                 lon_in=self.out[1],
+                                                 inverse=True)
+        assert np.all(abs(self.val['lat'] - self.loc[0]) < 1.0e-6)
+        assert np.all(abs(self.val['lon'] - self.loc[1]) < 1.0e-6)
+
+    @pytest.mark.parametrize("kwargs,target",
+                             [({}, [44.8075768, 8.0, 6367.48954386,
+                                    51.7037677, 62.8811403]),
+                              ({'inverse': False},
+                               [44.8075768, 8.0, 6367.48954386,
+                                51.7037677, 62.8811403]),
+                              ({'inverse': True},
+                               [45.1924232, 8.0, 6367.3459085,
+                                52.2989610, 63.1180720])])
+    def test_geodetic_to_geocentric_horizontal(self, kwargs, target):
+        """Test conversion from geodetic to geocentric coordinates"""
+
+        self.out = coords.geodetic_to_geocentric_horizontal(self.val['lat'],
+                                                            self.val['lon'],
+                                                            self.val['az'],
+                                                            self.val['el'],
+                                                            **kwargs)
+
+        for i, self.loc in enumerate(self.out):
+            assert np.all(abs(self.loc - target[i]) < 1.0e-6)
+            if isinstance(self.loc, np.ndarray):
+                assert self.loc.shape == self.val['lat'].shape
+
+    def test_geodetic_to_geocentric_horizontal_and_back(self):
+        """Tests the reversibility of geodetic to geocentric horiz conversions
+
+        Note:  inverse of az and el angles currently non-functional
+
+        """
+
+        self.out = coords.geodetic_to_geocentric_horizontal(self.val['lat'],
+                                                            self.val['lon'],
+                                                            self.val['az'],
+                                                            self.val['el'],
+                                                            inverse=False)
+        self.loc = coords.geodetic_to_geocentric_horizontal(self.out[0],
+                                                            self.out[1],
+                                                            self.out[3],
+                                                            self.out[4],
+                                                            inverse=True)
+
+        assert np.all(abs(self.val['lat'] - self.loc[0]) < 1.0e-6)
+        assert np.all(abs(self.val['lon'] - self.loc[1]) < 1.0e-6)
+        assert np.all(abs(self.val['az'] - self.loc[3]) < 1.0e-6)
+        assert np.all(abs(self.val['el'] - self.loc[4]) < 1.0e-6)
+
+
+class TestGeodeticGeocentricArray(TestGeodeticGeocentric):
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        arr = np.ones(shape=(10,), dtype=float)
+        self.val = {'lat': 45.0 * arr,
+                    'lon': 8.0 * arr,
+                    'az': 52.0 * arr,
+                    'el': 63.0 * arr}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+
+class TestSphereCartesian():
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        self.val = {'az': 45.0, 'el': 30.0, 'r': 1.0,
+                    'x': 0.6123724356957946,
+                    'y': 0.6123724356957946,
+                    'z': 0.5}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+    @pytest.mark.parametrize("kwargs,input,target",
+                             [({}, ['az', 'el', 'r'],
+                               ['x', 'y', 'z']),
+                              ({'inverse': False}, ['az', 'el', 'r'],
+                               ['x', 'y', 'z']),
+                              ({'inverse': True}, ['x', 'y', 'z'],
+                               ['az', 'el', 'r'])])
+    def test_spherical_to_cartesian(self, kwargs, input, target):
+        """Test conversion from spherical to cartesian coordinates"""
+
+        self.out = coords.spherical_to_cartesian(self.val[input[0]],
+                                                 self.val[input[1]],
+                                                 self.val[input[2]],
+                                                 **kwargs)
+
+        for i, self.loc in enumerate(self.out):
+            assert np.all(abs(self.loc - self.val[target[i]]) < 1.0e-6)
+            if isinstance(self.loc, np.ndarray):
+                assert self.loc.shape == self.val[input[0]].shape
+
+    def test_spherical_to_cartesian_and_back(self):
+        """Tests the reversibility of spherical to cartesian conversions"""
+
+        self.out = coords.spherical_to_cartesian(self.val['x'], self.val['y'],
+                                                 self.val['z'], inverse=True)
+        self.out = coords.spherical_to_cartesian(self.out[0], self.out[1],
+                                                 self.out[2], inverse=False)
+
+        assert np.all(abs(self.val['x'] - self.out[0]) < 1.0e-6)
+        assert np.all(abs(self.val['y'] - self.out[1]) < 1.0e-6)
+        assert np.all(abs(self.val['z'] - self.out[2]) < 1.0e-6)
+
+
+class TestSphereCartesianArray(TestSphereCartesian):
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        arr = np.ones(shape=(10,), dtype=float)
+        self.val = {'az': 45.0 * arr, 'el': 30.0 * arr, 'r': 1.0 * arr,
+                    'x': 0.6123724356957946 * arr,
+                    'y': 0.6123724356957946 * arr,
+                    'z': 0.5 * arr}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+
+class TestGlobalLocal():
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        self.val = {'x': 7000.0, 'y': 8000.0, 'z': 9000.0,
+                    'lat': 37.5, 'lon': 289.0, 'rad': 6380.0}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+    @pytest.mark.parametrize("kwargs,target",
+                             [({},
+                               [-9223.1752649, -2239.8352784, 11323.1268511]),
+                              ({'inverse': False},
+                               [-9223.1752649, -2239.8352784, 11323.1268511]),
+                              ({'inverse': True},
+                               [-5709.804677, -4918.397556, 15709.5775005])])
+    def test_global_to_local_cartesian(self, kwargs, target):
+        """Test conversion from global to local cartesian coordinates"""
+
+        self.out = coords.global_to_local_cartesian(self.val['x'],
+                                                    self.val['y'],
+                                                    self.val['z'],
+                                                    self.val['lat'],
+                                                    self.val['lon'],
+                                                    self.val['rad'],
+                                                    **kwargs)
+
+        for i, self.loc in enumerate(self.out):
+            assert np.all(abs(self.loc - target[i]) < 1.0e-6)
+            if isinstance(self.loc, np.ndarray):
+                assert self.loc.shape == self.val['x'].shape
+
+    def test_global_to_local_cartesian_and_back(self):
+        """Tests the reversibility of the global to loc cartesian transform"""
+
+        self.out = coords.global_to_local_cartesian(self.val['x'],
+                                                    self.val['y'],
+                                                    self.val['z'],
+                                                    self.val['lat'],
+                                                    self.val['lon'],
+                                                    self.val['rad'],
+                                                    inverse=False)
+        self.out = coords.global_to_local_cartesian(self.out[0], self.out[1],
+                                                    self.out[2],
+                                                    self.val['lat'],
+                                                    self.val['lon'],
+                                                    self.val['rad'],
+                                                    inverse=True)
+        assert np.all(abs(self.val['x'] - self.out[0]) < 1.0e-6)
+        assert np.all(abs(self.val['y'] - self.out[1]) < 1.0e-6)
+        assert np.all(abs(self.val['z'] - self.out[2]) < 1.0e-6)
+
+
+class TestGlobalLocalArray(TestGlobalLocal):
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        arr = np.ones(shape=(10,), dtype=float)
+        self.val = {'x': 7000.0 * arr, 'y': 8000.0 * arr, 'z': 9000.0 * arr,
+                    'lat': 37.5 * arr, 'lon': 289.0 * arr, 'rad': 6380.0 * arr}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+
+class TestLocalHorizontalGlobal():
+    """Tests for local horizontal to global geo and back """
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        self.val = {'az': 30.0, 'el': 45.0, 'dist': 1000.0,
+                    'lat': 45.0, 'lon': 0.0, 'alt': 400.0}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
+
+    @pytest.mark.parametrize("kwargs,target",
+                             [({}, [50.4190376, -7.6940084, 7172.1548652]),
+                              ({'geodetic': True},
+                               [50.4190376, -7.6940084, 7172.1548652]),
+                              ({'geodetic': False},
+                               [50.4143159, -7.6855552, 7185.6983666])])
+    def test_local_horizontal_to_global_geo(self, kwargs, target):
+        """Tests the conversion of the local horizontal to global geo"""
+
+        self.out = coords.local_horizontal_to_global_geo(self.val['az'],
+                                                         self.val['el'],
+                                                         self.val['dist'],
+                                                         self.val['lat'],
+                                                         self.val['lon'],
+                                                         self.val['alt'],
+                                                         **kwargs)
+
+        for i, self.loc in enumerate(self.out):
+            assert np.all(abs(self.loc - target[i]) < 1.0e-6)
+            if isinstance(self.loc, np.ndarray):
+                assert self.loc.shape == self.val['lat'].shape
+
+
+class TestLocalHorizontalGlobalArray(TestLocalHorizontalGlobal):
+    """Tests for local horizontal to global geo and back """
+
+    def setup(self):
+        """Runs before every method to create a clean testing setup."""
+        arr = np.ones(shape=(10,), dtype=float)
+        self.val = {'az': 30.0 * arr, 'el': 45.0 * arr, 'dist': 1000.0 * arr,
+                    'lat': 45.0 * arr, 'lon': 0.0 * arr, 'alt': 400.0 * arr}
+        self.out = None
+        self.loc = None
+
+    def teardown(self):
+        """Runs after every method to clean up previous testing."""
+        del self.val, self.out, self.loc
diff --git a/pysatMadrigal/utils/__init__.py b/pysatMadrigal/utils/__init__.py
new file mode 100644
index 0000000..0854b8d
--- /dev/null
+++ b/pysatMadrigal/utils/__init__.py
@@ -0,0 +1 @@
+from pysatMadrigal.utils import coords  # noqa F401
diff --git a/pysatMadrigal/utils/coords.py b/pysatMadrigal/utils/coords.py
new file mode 100644
index 0000000..ff2a0db
--- /dev/null
+++ b/pysatMadrigal/utils/coords.py
@@ -0,0 +1,352 @@
+"""
+pysatMadrigal.coords - coordinate transformations for radar instruments
+=========================================
+
+pysat.coords contains a number of coordinate-transformation
+functions used throughout the pysat package.
+"""
+
+import numpy as np
+
+
+def geodetic_to_geocentric(lat_in, lon_in=None, inverse=False):
+    """Converts position from geodetic to geocentric or vice-versa.
+
+    Parameters
+    ----------
+    lat_in : float
+        latitude in degrees.
+    lon_in : float or NoneType
+        longitude in degrees.  Remains unchanged, so does not need to be
+        included. (default=None)
+    inverse : bool
+        False for geodetic to geocentric, True for geocentric to geodetic.
+        (default=False)
+
+    Returns
+    -------
+    lat_out : float
+        latitude [degree] (geocentric/detic if inverse=False/True)
+    lon_out : float or NoneType
+        longitude [degree] (geocentric/detic if inverse=False/True)
+    rad_earth : float
+        Earth radius [km] (geocentric/detic if inverse=False/True)
+
+    Notes
+    -----
+    Uses WGS-84 values
+
+    References
+    ----------
+    Based on J.M. Ruohoniemi's geopack and R.J. Barnes radar.pro
+
+    """
+    rad_eq = 6378.1370  # WGS-84 semi-major axis
+    flat = 1.0 / 298.257223563  # WGS-84 flattening
+    rad_pol = rad_eq * (1.0 - flat)  # WGS-84 semi-minor axis
+
+    # The ratio between the semi-major and minor axis is used several times
+    rad_ratio_sq = (rad_eq / rad_pol)**2
+
+    # Calculate the square of the second eccentricity (e')
+    eprime_sq = rad_ratio_sq - 1.0
+
+    # Calculate the tangent of the input latitude
+    tan_in = np.tan(np.radians(lat_in))
+
+    # If converting from geodetic to geocentric, take the inverse of the
+    # radius ratio
+    if not inverse:
+        rad_ratio_sq = 1.0 / rad_ratio_sq
+
+    # Calculate the output latitude
+    lat_out = np.degrees(np.arctan(rad_ratio_sq * tan_in))
+
+    # Calculate the Earth radius at this latitude
+    rad_earth = rad_eq / np.sqrt(1.0 + eprime_sq
+                                 * np.sin(np.radians(lat_out))**2)
+
+    # longitude remains unchanged
+    lon_out = lon_in
+
+    return lat_out, lon_out, rad_earth
+
+
+def geodetic_to_geocentric_horizontal(lat_in, lon_in, az_in, el_in,
+                                      inverse=False):
+    """Converts from local horizontal coordinates in a geodetic system to local
+    horizontal coordinates in a geocentric system
+
+    Parameters
+    ----------
+    lat_in : float
+        latitude in degrees of the local horizontal coordinate system center
+    lon_in : float
+        longitude in degrees of the local horizontal coordinate system center
+    az_in : float
+        azimuth in degrees within the local horizontal coordinate system
+    el_in : float
+        elevation in degrees within the local horizontal coordinate system
+    inverse : bool
+        False for geodetic to geocentric, True for inverse (default=False)
+
+    Returns
+    -------
+    lat_out : float
+        latitude in degrees of the converted horizontal coordinate system
+        center
+    lon_out : float
+        longitude in degrees of the converted horizontal coordinate system
+        center
+    rad_earth : float
+        Earth radius in km at the geocentric/detic (False/True) location
+    az_out : float
+        azimuth in degrees of the converted horizontal coordinate system
+    el_out : float
+        elevation in degrees of the converted horizontal coordinate system
+
+    References
+    ----------
+    Based on J.M. Ruohoniemi's geopack and R.J. Barnes radar.pro
+
+    """
+
+    az = np.radians(az_in)
+    el = np.radians(el_in)
+
+    # Transform the location of the local horizontal coordinate system center
+    lat_out, lon_out, rad_earth = geodetic_to_geocentric(lat_in, lon_in,
+                                                         inverse=inverse)
+
+    # Calcualte the deviation from vertical in radians
+    dev_vert = np.radians(lat_in - lat_out)
+
+    # Calculate cartesian coordinated in local system
+    x_local = np.cos(el) * np.sin(az)
+    y_local = np.cos(el) * np.cos(az)
+    z_local = np.sin(el)
+
+    # Now rotate system about the x axis to align local vertical vector
+    # with Earth radial vector
+    x_out = x_local
+    y_out = y_local * np.cos(dev_vert) + z_local * np.sin(dev_vert)
+    z_out = -y_local * np.sin(dev_vert) + z_local * np.cos(dev_vert)
+
+    # Transform the azimuth and elevation angles
+    az_out = np.degrees(np.arctan2(x_out, y_out))
+    el_out = np.degrees(np.arctan(z_out / np.sqrt(x_out**2 + y_out**2)))
+
+    return lat_out, lon_out, rad_earth, az_out, el_out
+
+
+def spherical_to_cartesian(az_in, el_in, r_in, inverse=False):
+    """Convert a position from spherical to cartesian, or vice-versa
+
+    Parameters
+    ----------
+    az_in : float
+        azimuth/longitude in degrees or cartesian x in km (inverse=False/True)
+    el_in : float
+        elevation/latitude in degrees or cartesian y in km (inverse=False/True)
+    r_in : float
+        distance from origin in km or cartesian z in km (inverse=False/True)
+    inverse : boolian
+        False to go from spherical to cartesian and True for the inverse
+
+    Returns
+    -------
+    x_out : float
+        cartesian x in km or azimuth/longitude in degrees (inverse=False/True)
+    y_out : float
+        cartesian y in km or elevation/latitude in degrees (inverse=False/True)
+    z_out : float
+        cartesian z in km or distance from origin in km (inverse=False/True)
+
+    Notes
+    ------
+    This transform is the same for local or global spherical/cartesian
+    transformations.
+
+    Returns elevation angle (angle from the xy plane) rather than zenith angle
+    (angle from the z-axis)
+
+    """
+
+    if inverse:
+        # Cartesian to Spherical
+        xy_sq = az_in**2 + el_in**2
+        z_out = np.sqrt(xy_sq + r_in**2)  # This is r
+        y_out = np.degrees(np.arctan2(np.sqrt(xy_sq), r_in))  # This is zenith
+        y_out = 90.0 - y_out  # This is the elevation
+        x_out = np.degrees(np.arctan2(el_in, az_in))  # This is azimuth
+    else:
+        # Spherical coordinate system uses zenith angle (degrees from the
+        # z-axis) and not the elevation angle (degrees from the x-y plane)
+        zen_in = np.radians(90.0 - el_in)
+
+        # Spherical to Cartesian
+        x_out = r_in * np.sin(zen_in) * np.cos(np.radians(az_in))
+        y_out = r_in * np.sin(zen_in) * np.sin(np.radians(az_in))
+        z_out = r_in * np.cos(zen_in)
+
+    return x_out, y_out, z_out
+
+
+def global_to_local_cartesian(x_in, y_in, z_in, lat_cent, lon_cent, rad_cent,
+                              inverse=False):
+    """Converts a position from global to local cartesian or vice-versa
+
+    Parameters
+    ----------
+    x_in : float
+        global or local cartesian x in km (inverse=False/True)
+    y_in : float
+        global or local cartesian y in km (inverse=False/True)
+    z_in : float
+        global or local cartesian z in km (inverse=False/True)
+    lat_cent : float
+        geocentric latitude in degrees of local cartesian system origin
+    lon_cent : float
+        geocentric longitude in degrees of local cartesian system origin
+    rad_cent : float
+        distance from center of the Earth in km of local cartesian system
+        origin
+    inverse : bool
+        False to convert from global to local cartesian coodiantes, and True
+        for the inverse (default=False)
+
+    Returns
+    -------
+    x_out : float
+        local or global cartesian x in km (inverse=False/True)
+    y_out : float
+        local or global cartesian y in km (inverse=False/True)
+    z_out : float
+        local or global cartesian z in km (inverse=False/True)
+
+    Notes
+    -------
+    The global cartesian coordinate system has its origin at the center of the
+    Earth, while the local system has its origin specified by the input
+    latitude, longitude, and radius.  The global system has x intersecting
+    the equatorial plane and the prime meridian, z pointing North along the
+    rotational axis, and y completing the right-handed coodinate system.
+    The local system has z pointing up, y pointing North, and x pointing East.
+
+    """
+
+    # Get the global cartesian coordinates of local origin
+    x_cent, y_cent, z_cent = spherical_to_cartesian(lon_cent, lat_cent,
+                                                    rad_cent)
+
+    # Get the amount of rotation needed to align the x-axis with the
+    # Earth's rotational axis
+    ax_rot = np.radians(90.0 - lat_cent)
+
+    # Get the amount of rotation needed to align the global x-axis with the
+    # prime meridian
+    mer_rot = np.radians(lon_cent - 90.0)
+
+    if inverse:
+        # Rotate about the x-axis to align the z-axis with the Earth's
+        # rotational axis
+        xrot = x_in
+        yrot = y_in * np.cos(ax_rot) - z_in * np.sin(ax_rot)
+        zrot = y_in * np.sin(ax_rot) + z_in * np.cos(ax_rot)
+
+        # Rotate about the global z-axis to get the global x-axis aligned
+        # with the prime meridian and translate the local center to the
+        # global origin
+        x_out = xrot * np.cos(mer_rot) - yrot * np.sin(mer_rot) + x_cent
+        y_out = xrot * np.sin(mer_rot) + yrot * np.cos(mer_rot) + y_cent
+        z_out = zrot + z_cent
+    else:
+        # Translate global origin to the local origin
+        xtrans = x_in - x_cent
+        ytrans = y_in - y_cent
+        ztrans = z_in - z_cent
+
+        # Rotate about the global z-axis to get the local x-axis pointing East
+        xrot = xtrans * np.cos(-mer_rot) - ytrans * np.sin(-mer_rot)
+        yrot = xtrans * np.sin(-mer_rot) + ytrans * np.cos(-mer_rot)
+        zrot = ztrans
+
+        # Rotate about the x-axis to get the z-axis pointing up
+        x_out = xrot
+        y_out = yrot * np.cos(-ax_rot) - zrot * np.sin(-ax_rot)
+        z_out = yrot * np.sin(-ax_rot) + zrot * np.cos(-ax_rot)
+
+    return x_out, y_out, z_out
+
+
+def local_horizontal_to_global_geo(az, el, dist, lat_orig, lon_orig, alt_orig,
+                                   geodetic=True):
+    """ Convert from local horizontal coordinates to geodetic or geocentric
+    coordinates
+
+    Parameters
+    ----------
+    az : float
+        Azimuth (angle from North) of point in degrees
+    el : float
+        Elevation (angle from ground) of point in degrees
+    dist : float
+        Distance from origin to point in km
+    lat_orig : float
+        Latitude of origin in degrees
+    lon_orig : float
+        Longitude of origin in degrees
+    alt_orig : float
+        Altitude of origin in km from the surface of the Earth
+    geodetic : bool
+        True if origin coordinates are geodetic, False if they are geocentric.
+        Will return coordinates in the same system as the origin input.
+        (default=True)
+
+    Returns
+    -------
+    lat_pnt : float
+        Latitude of point in degrees
+    lon_pnt : float
+        Longitude of point in degrees
+    rad_pnt : float
+        Distance to the point from the centre of the Earth in km
+
+    References
+    ----------
+    Based on J.M. Ruohoniemi's geopack and R.J. Barnes radar.pro
+
+    """
+
+    # If the data are in geodetic coordiantes, convert to geocentric
+    if geodetic:
+        (glat, glon, rearth, gaz, gel) = \
+            geodetic_to_geocentric_horizontal(lat_orig, lon_orig, az, el,
+                                              inverse=False)
+        grad = rearth + alt_orig
+    else:
+        glat = lat_orig
+        glon = lon_orig
+        grad = alt_orig + 6371.0  # Add the mean earth radius in km
+        gaz = az
+        gel = el
+
+    # Convert from local horizontal to local cartesian coordiantes
+    x_loc, y_loc, z_loc = spherical_to_cartesian(gaz, gel, dist, inverse=False)
+
+    # Convert from local to global cartesian coordiantes
+    x_glob, y_glob, z_glob = global_to_local_cartesian(x_loc, y_loc, z_loc,
+                                                       glat, glon, grad,
+                                                       inverse=True)
+
+    # Convert from global cartesian to geocentric coordinates
+    lon_pnt, lat_pnt, rad_pnt = spherical_to_cartesian(x_glob, y_glob, z_glob,
+                                                       inverse=True)
+
+    # Convert from geocentric to geodetic, if desired
+    if geodetic:
+        lat_pnt, lon_pnt, rearth = geodetic_to_geocentric(lat_pnt, lon_pnt,
+                                                          inverse=True)
+        rad_pnt = rearth + rad_pnt - 6371.0
+
+    return lat_pnt, lon_pnt, rad_pnt
diff --git a/setup.cfg b/setup.cfg
index d1c319b..29eec56 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -2,6 +2,9 @@
 omit =
   */instruments/templates/
 
+[flake8]
+max-line-length = 80
+
 [tool:pytest]
 markers =
     all_inst: tests all instruments