apexpy.apex

Complexity

Total Complexity

91

Size/Duplication

Total Lines	1306
Duplicated Lines	5.05 %

Importance

Changes

0

27 Methods

Rating	Name	Duplication	Size	Complexity
D	Apex.__init__()	0	67	13
F	Apex.convert()	0	102	19
A	Apex.set_epoch()	0	16	1
A	Apex.qd2geo()	0	39	1
A	Apex.basevectors_qd()	0	61	2
A	Apex.geo2qd()	0	28	1
A	Apex.mlt2mlon()	0	40	1
C	Apex._map_EV_to_height()	0	75	11
A	Apex.map_V_to_height()	0	30	1
A	Apex.mlon2mlt()	0	42	1
B	Apex.basevectors_apex()	0	149	7
A	Apex.set_refh()	0	15	1
B	Apex.__eq__()	0	42	7
A	Apex.map_E_to_height()	0	30	1
A	Apex.qd2apex()	0	30	1
A	Apex.apex2geo()	0	39	1
A	Apex.__repr__()	0	8	1
A	Apex.geo2apex()	0	35	3
A	Apex.get_babs()	0	23	1
A	Apex._get_babs_nonvectorized()	0	28	1
A	Apex.__str__()	0	12	1
A	Apex.get_apex()	0	29	2
A	Apex._apex2qd_nonvectorized()	44	45	4
A	Apex._qd2apex_nonvectorized()	22	48	4
A	Apex.apex2qd()	0	31	1
A	Apex.bvectors_apex()	0	71	1
A	Apex.map_to_height()	0	55	3

# -*- coding: utf-8 -*-
"""
Classes that make up the core of apexpy.
"""

import datetime as dt
import numpy as np
import os
import warnings

from apexpy import helpers

# Below try..catch required for autodoc to work on readthedocs
try:
    from apexpy import fortranapex as fa
except ImportError:
    warnings.warn("".join(["fortranapex module could not be imported, so ",
                           "apexpy probably won't work.  Make sure you have ",
                           "a gfortran compiler."]))

# Make sure invalid warnings are always shown
warnings.filterwarnings('always', message='.*set to NaN where*',
                        module='apexpy.apex')


class ApexHeightError(ValueError):
    """Specialized ValueError definition, to be used when apex height is wrong.
    """
    pass


class Apex(object):
    """Calculates coordinate conversions, field-line mapping, and base vectors.

    Parameters
    ----------
    date : NoneType, float, :class:`dt.date`, or :class:`dt.datetime`, optional
        Determines which IGRF coefficients are used in conversions. Uses
        current date as default.  If float, use decimal year.  If None, uses
        current UTC. (default=None)
    refh : float, optional
        Reference height in km for apex coordinates, the field lines are mapped
        to this height. (default=0)
    datafile : str or NoneType, optional
        Path to custom coefficient file, if None uses `apexsh.dat` file
        (default=None)
    fortranlib : str or NoneType, optional
        Path to Fortran Apex CPython library, if None uses linked library file
        (default=None)

    Attributes
    ----------
    year : float
        Decimal year used for the IGRF model
    RE : float
        Earth radius in km, defaults to mean Earth radius
    refh : float
        Reference height in km for apex coordinates
    datafile : str
        Path to coefficient file
    fortranlib : str
        Path to Fortran Apex CPython library
    igrf_fn : str
        IGRF coefficient filename

    Notes
    -----
    The calculations use IGRF-13 with coefficients from 1900 to 2025 [1]_.

    The geodetic reference ellipsoid is WGS84.

    References
    ----------

    .. [1] Thébault, E. et al. (2015), International Geomagnetic Reference
           Field: the 12th generation, Earth, Planets and Space, 67(1), 79,
           :doi:`10.1186/s40623-015-0228-9`.

    """

    # ------------------------
    # Define the magic methods

    def __init__(self, date=None, refh=0, datafile=None, fortranlib=None):

        if datafile is None:
            datafile = os.path.join(os.path.dirname(__file__), 'apexsh.dat')

        if fortranlib is None:
            fortranlib = fa.__file__

        self.RE = 6371.009  # Mean Earth radius in km
        self.set_refh(refh)  # Reference height in km

        if date is None:
            self.year = helpers.toYearFraction(dt.datetime.utcnow())
        else:
            try:
                # Convert date/datetime object to decimal year
                self.year = helpers.toYearFraction(date)
            except AttributeError:
                # Failed while finding datetime attribute, so
                # date is probably an int or float; use directly
                self.year = date

        if not os.path.isfile(datafile):
            raise IOError('Data file does not exist: {}'.format(datafile))

        if not os.path.isfile(fortranlib):
            raise IOError('Fortran library does not exist: {}'.format(
                fortranlib))

        self.datafile = datafile
        self.fortranlib = fortranlib

        # Set the IGRF coefficient text file name
        self.igrf_fn = os.path.join(os.path.dirname(__file__),
                                    'igrf13coeffs.txt')
        if not os.path.exists(self.igrf_fn):
            raise OSError("File {} does not exist".format(self.igrf_fn))

        # Update the Fortran epoch using the year defined above
        self.set_epoch(self.year)

        # Vectorize the fortran functions
        self._geo2qd = np.frompyfunc(
            lambda glat, glon, height: fa.apxg2q(glat, (glon + 180) % 360 - 180,
                                                 height, 0)[:2], 3, 2)
        self._geo2apex = np.frompyfunc(
            lambda glat, glon, height: fa.apxg2all(glat, (glon + 180) % 360
                                                   - 180, height, self.refh,
                                                   0)[2:4], 3, 2)
        self._geo2apexall = np.frompyfunc(
            lambda glat, glon, height: fa.apxg2all(glat, (glon + 180) % 360
                                                   - 180, height, self.refh,
                                                   1), 3, 14)
        self._qd2geo = np.frompyfunc(
            lambda qlat, qlon, height, precision: fa.apxq2g(qlat, (qlon + 180)
                                                            % 360 - 180, height,
                                                            precision), 4, 3)
        self._basevec = np.frompyfunc(
            lambda glat, glon, height: fa.apxg2q(glat, (glon + 180) % 360 - 180,
                                                 height, 1)[2:4], 3, 2)

        # Vectorize other nonvectorized functions
        self._apex2qd = np.frompyfunc(self._apex2qd_nonvectorized, 3, 2)
        self._qd2apex = np.frompyfunc(self._qd2apex_nonvectorized, 3, 2)
        self._get_babs = np.frompyfunc(self._get_babs_nonvectorized, 3, 1)

        return

    def __repr__(self):
        """Produce an evaluatable representation of the Apex class."""
        rstr = "".join(["apexpy.Apex(date=", self.year.__repr__(), ", refh=",
                        self.refh.__repr__(), ", datafile=",
                        self.datafile.__repr__(), ", fortranlib=",
                        self.fortranlib.__repr__(), ")"])

        return rstr

    def __str__(self):
        """Produce a user-friendly representation of the Apex class."""

        out_str = '\n'.join(["Apex class conversions performed with",
                             "-------------------------------------",
                             "Decimal year: {:.8f}".format(self.year),
                             "Reference height: {:.3f} km".format(self.refh),
                             "Earth radius: {:.3f} km".format(self.RE), '',
                             "Coefficient file: '{:s}'".format(self.datafile),
                             "Cython Fortran library: '{:s}'".format(
                                 self.fortranlib)])
        return out_str

    def __eq__(self, comp_obj):
        """Performs equivalency evaluation.

        Parameters
        ----------
        comp_obj
            Object of any time to be compared to the class object

        Returns
        -------
        bool or NotImplemented
            True if self and comp_obj are identical, False if they are not,
            and NotImplemented if the classes are not the same

        """

        if isinstance(comp_obj, self.__class__):
            # Objects are the same if all the attributes are the same
            for apex_attr in ['year', 'refh', 'RE', 'datafile', 'fortranlib',
                              'igrf_fn']:
                bad_attr = False
                if hasattr(self, apex_attr):
                    aval = getattr(self, apex_attr)

                    if hasattr(comp_obj, apex_attr):
                        cval = getattr(comp_obj, apex_attr)

                        if cval != aval:
                            bad_attr = True
                    else:
                        bad_attr = True
                else:
                    bad_attr = True

                if bad_attr:
                    return False

            # If we arrive here, all expected attributes exist in both classes
            # and have the same values
            return True

        return NotImplemented

    # -------------------------
    # Define the hidden methods

    def _apex2qd_nonvectorized(self, alat, alon, height):

This code seems to be duplicated in your project.

        """Convert from apex to quasi-dipole (not-vectorised)

        Parameters
        -----------
        alat : (float)
            Apex latitude in degrees
        alon : (float)
            Apex longitude in degrees
        height : (float)
            Height in km

        Returns
        ---------
        qlat : (float)
            Quasi-dipole latitude in degrees
        qlon : (float)
            Quasi-diplole longitude in degrees

        """
        # Evaluate the latitude
        alat = helpers.checklat(alat, name='alat')

        # Convert modified apex to quasi-dipole, longitude is the same
        qlon = alon

        # Get the apex height
        h_apex = self.get_apex(alat)

        if h_apex < height:
            if np.isclose(h_apex, height, rtol=0, atol=1e-5):
                # Allow for values that are close
                h_apex = height
            else:
                estr = ''.join(['height {:.3g} is > '.format(np.max(height)),
                                'apex height {:.3g} for alat {:.3g}'.format(
                                    h_apex, alat)])
                raise ApexHeightError(estr)

        # Convert the latitude
        salat = np.sign(alat) if alat != 0 else 1
        qlat = salat * np.degrees(np.arccos(np.sqrt((self.RE + height) /
                                                    (self.RE + h_apex))))

        return qlat, qlon

    def _qd2apex_nonvectorized(self, qlat, qlon, height):

This code seems to be duplicated in your project.

        """Converts quasi-dipole to modified apex coordinates.

        Parameters
        ----------
        qlat : float
            Quasi-dipole latitude
        qlon : float
            Quasi-dipole longitude
        height : float
            Altitude in km

        Returns
        -------
        alat : float
            Modified apex latitude
        alon : float
            Modified apex longitude

        Raises
        ------
        ApexHeightError
            if apex height < reference height

        """
        # Evaluate the latitude
        qlat = helpers.checklat(qlat, name='qlat')

        # Get the longitude and apex height
        alon = qlon
        h_apex = self.get_apex(qlat, height)

        if h_apex < self.refh:
            if np.isclose(h_apex, self.refh, rtol=0, atol=1e-5):
                # Allow for values that are close
                h_apex = self.refh
            else:
                estr = ''.join(['apex height ({:.3g}) is < '.format(h_apex),
                                'reference height ({:.3g})'.format(self.refh),
                                ' for qlat {:.3g}'.format(qlat)])
                raise ApexHeightError(estr)

        # Convert the latitude
        sqlat = np.sign(qlat) if qlat != 0 else 1
        alat = sqlat * np.degrees(np.arccos(np.sqrt((self.RE + self.refh) /
                                                    (self.RE + h_apex))))

        return alat, alon

    def _map_EV_to_height(self, alat, alon, height, newheight, data, ev_flag):
        """Maps electric field related values to a desired height

        Parameters
        ----------
        alat : array-like
            Apex latitude in degrees.
        alon : array-like
            Apex longitude in degrees.
        height : array-like
            Current altitude in km.
        new_height : array-like
            Desired altitude to which EV values will be mapped in km.
        data : array-like
            3D value(s) for the electric field or electric drift
        ev_flag : str
           Specify if value is an electric field ('E') or electric drift ('V')

        Returns
        -------
        data_mapped : array-like
            Data mapped along the magnetic field from the old height to the
            new height.

        Raises
        ------
        ValueError
            If an unknown `ev_flag` or badly shaped data input is supplied.

        """
        # Ensure the E-V flag is correct
        ev_flag = ev_flag.upper()
        if ev_flag not in ['E', 'V']:
            raise ValueError('unknown electric field/drift flag')

        # Make sure data is array of correct shape
        if(not (np.ndim(data) == 1 and np.size(data) == 3)
           and not (np.ndim(data) == 2 and np.shape(data)[0] == 3)):
            # Raise ValueError because if passing e.g. a (6,) ndarray the
            # reshape below will work even though the input is invalid
            raise ValueError('{:} must be (3, N) or (3,) ndarray'.format(
                ev_flag))

        data = np.reshape(data, (3, np.size(data) // 3))

        # Get the necessary base vectors at the current and new height
        v1 = list()
        v2 = list()
        for i, alt in enumerate([height, newheight]):
            _, _, _, _, _, _, d1, d2, _, e1, e2, _ = self.basevectors_apex(
                alat, alon, alt, coords='apex')

            if ev_flag == 'E' and i == 0 or ev_flag == 'V' and i == 1:
                v1.append(e1)
                v2.append(e2)
            else:
                v1.append(d1)
                v2.append(d2)

            # Make sure v1 and v2 have shape (3, N)
            v1[-1] = np.reshape(v1[-1], (3, v1[-1].size // 3))
            v2[-1] = np.reshape(v2[-1], (3, v2[-1].size // 3))

        # Take the dot product between the data value and each base vector
        # at the current height
        data1 = np.sum(data * v1[0], axis=0)
        data2 = np.sum(data * v2[0], axis=0)
        print("BI", data, v1[0], np.sum(data * v1[0], axis=0), data1)

        # Map the data to the new height, removing any axes of length 1
        # after the calculation
        data_mapped = np.squeeze(data1[np.newaxis, :] * v1[1]
                                 + data2[np.newaxis, :] * v2[1])

        return data_mapped

    def _get_babs_nonvectorized(self, glat, glon, height):
        """Get the absolute value of the B-field in Tesla

        Parameters
        ----------
        glat : float
            Geodetic latitude in degrees
        glon : float
            Geodetic longitude in degrees
        height : float
            Altitude in km

        Returns
        -------
        babs : float
            Absolute value of the magnetic field in Tesla

        """
        # Evaluate the latitude
        glat = helpers.checklat(glat, name='qlat')

        # Get the magnetic field output: North, East, Down, Absolute value
        bout = fa.feldg(1, glat, glon, height)

        # Convert the absolute value from Gauss to Tesla
        babs = bout[-1] / 10000.0

        return babs

    # -----------------------
    # Define the user methods

    def convert(self, lat, lon, source, dest, height=0, datetime=None,
                precision=1e-10, ssheight=50 * 6371):
        """Converts between geodetic, modified apex, quasi-dipole and MLT.

        Parameters
        ----------
        lat : array_like
            Latitude in degrees
        lon : array_like
            Longitude in degrees or MLT in hours
        source : str
            Input coordinate system, accepts 'geo', 'apex', 'qd', or 'mlt'
        dest : str
            Output coordinate system, accepts 'geo', 'apex', 'qd', or 'mlt'
        height : array_like, optional
            Altitude in km
        datetime : :class:`datetime.datetime`
            Date and time for MLT conversions (required for MLT conversions)
        precision : float, optional
            Precision of output (degrees) when converting to geo. A negative
            value of this argument produces a low-precision calculation of
            geodetic lat/lon based only on their spherical harmonic
            representation.
            A positive value causes the underlying Fortran routine to iterate
            until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces
            the input QD lat/lon to within the specified precision (all
            coordinates being converted to geo are converted to QD first and
            passed through APXG2Q).
        ssheight : float, optional
            Altitude in km to use for converting the subsolar point from
            geographic to magnetic coordinates. A high altitude is used
            to ensure the subsolar point is mapped to high latitudes, which
            prevents the South-Atlantic Anomaly (SAA) from influencing the MLT.

        Returns
        -------
        lat : ndarray or float
            Converted latitude (if converting to MLT, output latitude is apex)
        lon : ndarray or float
            Converted longitude or MLT

        Raises
        ------
        ValueError
            For unknown source or destination coordinate system or a missing
            or badly formed latitude or datetime input

        """
        # Test the input values
        source = source.lower()
        dest = dest.lower()

        if source not in ['geo', 'apex', 'qd', 'mlt'] \
           or dest not in ['geo', 'apex', 'qd', 'mlt']:
            estr = 'Unknown coordinate transformation: {} -> {}'.format(source,
                                                                        dest)
            raise ValueError(estr)

        if datetime is None and ('mlt' in [source, dest]):
            raise ValueError('datetime must be given for MLT calculations')

        lat = helpers.checklat(lat)

        if source == dest:
            # The source and destination are the same, no conversion necessary
            return lat, lon
        else:
            # Select the correct conversion method
            if source == 'geo':
                if dest == 'qd':
                    lat, lon = self.geo2qd(lat, lon, height)
                else:
                    lat, lon = self.geo2apex(lat, lon, height)
                    if dest == 'mlt':
                        lon = self.mlon2mlt(lon, datetime, ssheight=ssheight)
            elif source == 'apex':
                if dest == 'geo':
                    lat, lon, _ = self.apex2geo(lat, lon, height,
                                                precision=precision)
                elif dest == 'qd':
                    lat, lon = self.apex2qd(lat, lon, height=height)
                elif dest == 'mlt':
                    lon = self.mlon2mlt(lon, datetime, ssheight=ssheight)
            elif source == 'qd':
                if dest == 'geo':
                    lat, lon, _ = self.qd2geo(lat, lon, height,
                                              precision=precision)
                else:
                    lat, lon = self.qd2apex(lat, lon, height=height)
                    if dest == 'mlt':
                        lon = self.mlon2mlt(lon, datetime, ssheight=ssheight)
            elif source == 'mlt':
                # From MLT means that the input latitude is assumed to be apex,
                # so we don't need to update latitude for apex conversions.
                lon = self.mlt2mlon(lon, datetime, ssheight=ssheight)
                if dest == 'geo':
                    lat, lon, _ = self.apex2geo(lat, lon, height,
                                                precision=precision)
                elif dest == 'qd':
                    lat, lon = self.apex2qd(lat, lon, height=height)

        return lat, lon

    def geo2apex(self, glat, glon, height):
        """Converts geodetic to modified apex coordinates.

        Parameters
        ----------
        glat : array_like
            Geodetic latitude
        glon : array_like
            Geodetic longitude
        height : array_like
            Altitude in km

        Returns
        -------
        alat : ndarray or float
            Modified apex latitude
        alon : ndarray or float
            Modified apex longitude

        """

        glat = helpers.checklat(glat, name='glat')

        alat, alon = self._geo2apex(glat, glon, height)

        if np.any(alat == -9999):
            warnings.warn(''.join(['Apex latitude set to NaN where undefined ',
                                   '(apex height may be < reference height)']))
            if np.isscalar(alat):
                alat = np.nan
            else:
                alat[alat == -9999] = np.nan

        # If array is returned, dtype is object, so convert to float
        return np.float64(alat), np.float64(alon)

    def apex2geo(self, alat, alon, height, precision=1e-10):
        """Converts modified apex to geodetic coordinates.

        Parameters
        ----------
        alat : array_like
            Modified apex latitude
        alon : array_like
            Modified apex longitude
        height : array_like
            Altitude in km
        precision : float, optional
            Precision of output (degrees). A negative value of this argument
            produces a low-precision calculation of geodetic lat/lon based only
            on their spherical harmonic representation. A positive value causes
            the underlying Fortran routine to iterate until feeding the output
            geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to
            within the specified precision.

        Returns
        -------
        glat : ndarray or float
            Geodetic latitude
        glon : ndarray or float
            Geodetic longitude
        error : ndarray or float
            The angular difference (degrees) between the input QD coordinates
            and the qlat/qlon produced by feeding the output glat and glon
            into geo2qd (APXG2Q)

        """
        # Evaluate the latitude
        alat = helpers.checklat(alat, name='alat')

        # Perform the two-step convertion to geodetic coordinates
        qlat, qlon = self.apex2qd(alat, alon, height=height)
        glat, glon, error = self.qd2geo(qlat, qlon, height, precision=precision)

        return glat, glon, error

    def geo2qd(self, glat, glon, height):
        """Converts geodetic to quasi-dipole coordinates.

        Parameters
        ----------
        glat : array_like
            Geodetic latitude
        glon : array_like
            Geodetic longitude
        height : array_like
            Altitude in km

        Returns
        -------
        qlat : ndarray or float
            Quasi-dipole latitude
        qlon : ndarray or float
            Quasi-dipole longitude

        """
        # Evaluate the latitude
        glat = helpers.checklat(glat, name='glat')

        # Convert to quasi-dipole coordinates
        qlat, qlon = self._geo2qd(glat, glon, height)

        # If array is returned, dtype is object, so convert to float
        return np.float64(qlat), np.float64(qlon)

    def qd2geo(self, qlat, qlon, height, precision=1e-10):
        """Converts quasi-dipole to geodetic coordinates.

        Parameters
        ----------
        qlat : array_like
            Quasi-dipole latitude
        qlon : array_like
            Quasi-dipole longitude
        height : array_like
            Altitude in km
        precision : float, optional
            Precision of output (degrees). A negative value of this argument
            produces a low-precision calculation of geodetic lat/lon based only
            on their spherical harmonic representation. A positive value causes
            the underlying Fortran routine to iterate until feeding the output
            geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to
            within the specified precision.

        Returns
        -------
        glat : ndarray or float
            Geodetic latitude
        glon : ndarray or float
            Geodetic longitude
        error : ndarray or float
            The angular difference (degrees) between the input QD coordinates
            and the qlat/qlon produced by feeding the output glat and glon
            into geo2qd (APXG2Q)

        """
        # Evaluate the latitude
        qlat = helpers.checklat(qlat, name='qlat')

        # Convert to geodetic coordinates
        glat, glon, error = self._qd2geo(qlat, qlon, height, precision)

        # If array is returned, dtype is object, so convert to float
        return np.float64(glat), np.float64(glon), np.float64(error)

    def apex2qd(self, alat, alon, height):
        """Converts modified apex to quasi-dipole coordinates.

        Parameters
        ----------
        alat : array_like
            Modified apex latitude
        alon : array_like
            Modified apex longitude
        height : array_like
            Altitude in km

        Returns
        -------
        qlat : ndarray or float
            Quasi-dipole latitude
        qlon : ndarray or float
            Quasi-dipole longitude

        Raises
        ------
        ApexHeightError
            if `height` > apex height

        """
        # Convert the latitude to apex, latitude check performed in the
        # hidden method
        qlat, qlon = self._apex2qd(alat, alon, height)

        # If array is returned, the dtype is object, so convert to float
        return np.float64(qlat), np.float64(qlon)

    def qd2apex(self, qlat, qlon, height):
        """Converts quasi-dipole to modified apex coordinates.

        Parameters
        ----------
        qlat : array_like
            Quasi-dipole latitude
        qlon : array_like
            Quasi-dipole longitude
        height : array_like
            Altitude in km

        Returns
        -------
        alat : ndarray or float
            Modified apex latitude
        alon : ndarray or float
            Modified apex longitude

        Raises
        ------
        ApexHeightError
            if apex height < reference height

        """
        # Perform the conversion from quasi-dipole to apex coordinates
        alat, alon = self._qd2apex(qlat, qlon, height)

        # If array is returned, the dtype is object, so convert to float
        return np.float64(alat), np.float64(alon)

    def mlon2mlt(self, mlon, dtime, ssheight=318550):
        """Computes the magnetic local time at the specified magnetic longitude
        and UT.

        Parameters
        ----------
        mlon : array_like
            Magnetic longitude (apex and quasi-dipole longitude are always
            equal)
        dtime : :class:`datetime.datetime`
            Date and time
        ssheight : float, optional
            Altitude in km to use for converting the subsolar point from
            geographic to magnetic coordinates. A high altitude is used
            to ensure the subsolar point is mapped to high latitudes, which
            prevents the South-Atlantic Anomaly (SAA) from influencing the MLT.
            The current default is  50 * 6371, roughly 50 RE. (default=318550)

        Returns
        -------
        mlt : ndarray or float
            Magnetic local time in hours [0, 24)

        Notes
        -----
        To compute the MLT, we find the apex longitude of the subsolar point at
        the given time. Then the MLT of the given point will be computed from
        the separation in magnetic longitude from this point (1 hour = 15
        degrees).

        """
        # Get the subsolar location
        ssglat, ssglon = helpers.subsol(dtime)

        # Convert the subsolar location to apex coordinates
        _, ssalon = self.geo2apex(ssglat, ssglon, ssheight)

        # Calculate the magnetic local time (0-24 h range) from apex longitude.
        # np.float64 will ensure lists are converted to arrays
        mlt = (180 + np.float64(mlon) - ssalon) / 15 % 24

        return mlt

    def mlt2mlon(self, mlt, dtime, ssheight=318550):
        """Computes the magnetic longitude at the specified MLT and UT.

        Parameters
        ----------
        mlt : array_like
            Magnetic local time
        dtime : :class:`datetime.datetime`
            Date and time
        ssheight : float, optional
            Altitude in km to use for converting the subsolar point from
            geographic to magnetic coordinates. A high altitude is used
            to ensure the subsolar point is mapped to high latitudes, which
            prevents the South-Atlantic Anomaly (SAA) from influencing the MLT.
            The current default is  50 * 6371, roughly 50 RE. (default=318550)

        Returns
        -------
        mlon : ndarray or float
            Magnetic longitude [0, 360) (apex and quasi-dipole longitude are
            always equal)

        Notes
        -----
        To compute the magnetic longitude, we find the apex longitude of the
        subsolar point at the given time. Then the magnetic longitude of the
        given point will be computed from the separation in magnetic local time
        from this point (1 hour = 15 degrees).
        """
        # Get the location of the subsolar point at this time
        ssglat, ssglon = helpers.subsol(dtime)

        # Convert the location of the subsolar point to apex coordinates
        _, ssalon = self.geo2apex(ssglat, ssglon, ssheight)

        # Calculate the magnetic longitude (0-360 h range) from MLT.
        # np.float64 will ensure lists are converted to arrays
        mlon = (15 * np.float64(mlt) - 180 + ssalon + 360) % 360

        return mlon

    def map_to_height(self, glat, glon, height, newheight, conjugate=False,
                      precision=1e-10):
        """Performs mapping of points along the magnetic field to the closest
        or conjugate hemisphere.

        Parameters
        ----------
        glat : array_like
            Geodetic latitude
        glon : array_like
            Geodetic longitude
        height : array_like
            Source altitude in km
        newheight : array_like
            Destination altitude in km
        conjugate : bool, optional
            Map to `newheight` in the conjugate hemisphere instead of the
            closest hemisphere
        precision : float, optional
            Precision of output (degrees). A negative value of this argument
            produces a low-precision calculation of geodetic lat/lon based only
            on their spherical harmonic representation. A positive value causes
            the underlying Fortran routine to iterate until feeding the output
            geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to
            within the specified precision.

        Returns
        -------
        newglat : ndarray or float
            Geodetic latitude of mapped point
        newglon : ndarray or float
            Geodetic longitude of mapped point
        error : ndarray or float
            The angular difference (degrees) between the input QD coordinates
            and the qlat/qlon produced by feeding the output glat and glon
            into geo2qd (APXG2Q)

        Notes
        -----
        The mapping is done by converting glat/glon/height to modified apex
        lat/lon, and converting back to geographic using newheight (if
        conjugate, use negative apex latitude when converting back)

        """

        alat, alon = self.geo2apex(glat, glon, height)
        if conjugate:
            alat = -alat
        try:
            newglat, newglon, error = self.apex2geo(alat, alon, newheight,
                                                    precision=precision)
        except ApexHeightError:
            raise ApexHeightError("input height is > apex height")

        return newglat, newglon, error

    def map_E_to_height(self, alat, alon, height, newheight, edata):
        """Performs mapping of electric field along the magnetic field.

        It is assumed that the electric field is perpendicular to B.

        Parameters
        ----------
        alat : (N,) array_like or float
            Modified apex latitude
        alon : (N,) array_like or float
            Modified apex longitude
        height : (N,) array_like or float
            Source altitude in km
        newheight : (N,) array_like or float
            Destination altitude in km
        edata : (3,) or (3, N) array_like
            Electric field (at `alat`, `alon`, `height`) in geodetic east,
            north, and up components

        Returns
        -------
        out : (3, N) or (3,) ndarray
            The electric field at `newheight` (geodetic east, north, and up
            components)

        """
        # Call hidden mapping method with flag for the electric field
        out = self._map_EV_to_height(alat, alon, height, newheight, edata, 'E')

        return out

    def map_V_to_height(self, alat, alon, height, newheight, vdata):
        """Performs mapping of electric drift velocity along the magnetic field.

        It is assumed that the electric field is perpendicular to B.

        Parameters
        ----------
        alat : (N,) array_like or float
            Modified apex latitude
        alon : (N,) array_like or float
            Modified apex longitude
        height : (N,) array_like or float
            Source altitude in km
        newheight : (N,) array_like or float
            Destination altitude in km
        vdata : (3,) or (3, N) array_like
            Electric drift velocity (at `alat`, `alon`, `height`) in geodetic
            east, north, and up components

        Returns
        -------
        out : (3, N) or (3,) ndarray
            The electric drift velocity at `newheight` (geodetic east, north,
            and up components)

        """
        # Call hidden mapping method with flag for the electric drift velocities
        out = self._map_EV_to_height(alat, alon, height, newheight, vdata, 'V')

        return out

    def basevectors_qd(self, lat, lon, height, coords='geo', precision=1e-10):
        """Get quasi-dipole base vectors f1 and f2 at the specified coordinates.

        Parameters
        ----------
        lat : (N,) array_like or float
            Latitude
        lon : (N,) array_like or float
            Longitude
        height : (N,) array_like or float
            Altitude in km
        coords : {'geo', 'apex', 'qd'}, optional
            Input coordinate system
        precision : float, optional
            Precision of output (degrees) when converting to geo. A negative
            value of this argument produces a low-precision calculation of
            geodetic lat/lon based only on their spherical harmonic
            representation.
            A positive value causes the underlying Fortran routine to iterate
            until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces
            the input QD lat/lon to within the specified precision (all
            coordinates being converted to geo are converted to QD first and
            passed through APXG2Q).

        Returns
        -------
        f1 : (2, N) or (2,) ndarray
        f2 : (2, N) or (2,) ndarray

        Notes
        -----
        The vectors are described by Richmond [1995] [2]_ and
        Emmert et al. [2010] [3]_.  The vector components are geodetic east and
        north.

        References
        ----------
        .. [2] Richmond, A. D. (1995), Ionospheric Electrodynamics Using
               Magnetic Apex Coordinates, Journal of geomagnetism and
               geoelectricity, 47(2), 191–212, :doi:`10.5636/jgg.47.191`.

        .. [3] Emmert, J. T., A. D. Richmond, and D. P. Drob (2010),
               A computationally compact representation of Magnetic-Apex
               and Quasi-Dipole coordinates with smooth base vectors,
               J. Geophys. Res., 115(A8), A08322, :doi:`10.1029/2010JA015326`.

        """
        # Convert from current coordinates to geodetic coordinates
        glat, glon = self.convert(lat, lon, coords, 'geo', height=height,
                                  precision=precision)

        # Get the geodetic base vectors
        f1, f2 = self._basevec(glat, glon, height)

        # If inputs are not scalar, each vector is an array of arrays,
        # so reshape to a single array
        if f1.dtype == object:
            f1 = np.vstack(f1).T
            f2 = np.vstack(f2).T

        return f1, f2

    def basevectors_apex(self, lat, lon, height, coords='geo', precision=1e-10):
        """Returns base vectors in quasi-dipole and apex coordinates.

        Parameters
        ----------
        lat : array_like or float
            Latitude in degrees; input must be broadcastable with `lon` and
            `height`.
        lon : array_like or float
            Longitude in degrees; input must be broadcastable with `lat` and
            `height`.
        height : array_like or float
            Altitude in km; input must be broadcastable with `lon` and `lat`.
        coords : str, optional
            Input coordinate system, expects one of 'geo', 'apex', or 'qd'
            (default='geo')
        precision : float, optional
            Precision of output (degrees) when converting to geo. A negative
            value of this argument produces a low-precision calculation of
            geodetic lat/lon based only on their spherical harmonic
            representation.
            A positive value causes the underlying Fortran routine to iterate
            until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces
            the input QD lat/lon to within the specified precision (all
            coordinates being converted to geo are converted to QD first and
            passed through APXG2Q).

        Returns
        -------
        f1 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to e1, tangent to contours of
            constant lambdaA
        f2 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to e2
        f3 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to e3, tangent to contours of
            PhiA
        g1 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to d1
        g2 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to d2
        g3 : (3, N) or (3,) ndarray
            Quasi-dipole base vector equivalent to d3
        d1 : (3, N) or (3,) ndarray
            Apex base vector normal to contours of constant PhiA
        d2 : (3, N) or (3,) ndarray
            Apex base vector that completes the right-handed system
        d3 : (3, N) or (3,) ndarray
            Apex base vector normal to contours of constant lambdaA
        e1 : (3, N) or (3,) ndarray
            Apex base vector tangent to contours of constant V0
        e2 : (3, N) or (3,) ndarray
            Apex base vector that completes the right-handed system
        e3 : (3, N) or (3,) ndarray
            Apex base vector tangent to contours of constant PhiA

        Notes
        -----
        The vectors are described by Richmond [1995] [4]_ and
        Emmert et al. [2010] [5]_.  The vector components are geodetic east,
        north, and up (only east and north for `f1` and `f2`).

        `f3`, `g1`, `g2`, and `g3` are not part of the Fortran code
        by Emmert et al. [2010] [5]_. They are calculated by this
        Python library according to the following equations in
        Richmond [1995] [4]_:

        * `g1`: Eqn. 6.3
        * `g2`: Eqn. 6.4
        * `g3`: Eqn. 6.5
        * `f3`: Eqn. 6.8

        References
        ----------
        .. [4] Richmond, A. D. (1995), Ionospheric Electrodynamics Using
               Magnetic Apex Coordinates, Journal of geomagnetism and
               geoelectricity, 47(2), 191–212, :doi:`10.5636/jgg.47.191`.

        .. [5] Emmert, J. T., A. D. Richmond, and D. P. Drob (2010),
               A computationally compact representation of Magnetic-Apex
               and Quasi-Dipole coordinates with smooth base vectors,
               J. Geophys. Res., 115(A8), A08322, :doi:`10.1029/2010JA015326`.

        """
        # Convert to geodetic coordinates from current coordinate system
        glat, glon = self.convert(lat, lon, coords, 'geo', height=height,
                                  precision=precision)

        # Retrieve the desired magnetic locations and base vectors
        returnvals = list(self._geo2apexall(glat, glon, height))
        qlat = np.float64(returnvals[0])
        alat = np.float64(returnvals[2])
        bvec_ind = [4, 5, 7, 8, 9, 11, 12, 13]

        # If inputs are not scalar, each vector is an array of arrays,
        # so reshape to a single array
        if returnvals[4].dtype == object:
            for i in bvec_ind:
                returnvals[i] = np.vstack(returnvals[i]).T

        # Make sure the base vector arrays are 2D
        for i in bvec_ind:
            if i in [4, 5]:
                rsize = (2, returnvals[i].size // 2)
            else:
                rsize = (3, returnvals[i].size // 3)
            returnvals[i] = returnvals[i].reshape(rsize)

        # Assign the reshaped base vectors
        f1, f2 = returnvals[4:6]
        d1, d2, d3 = returnvals[7:10]
        e1, e2, e3 = returnvals[11:14]

        # Compute f3, g1, g2, g3 (outstanding quasi-dipole base vectors)
        #
        # Start by calculating the D equivalent, F (F_scalar)
        f1_stack = np.vstack((f1, np.zeros_like(f1[0])))
        f2_stack = np.vstack((f2, np.zeros_like(f2[0])))
        F_scalar = np.cross(f1_stack.T, f2_stack.T).T[-1]

        # Get the cosine of the magnetic inclination
        cos_mag_inc = helpers.getcosIm(alat)

        # Define the k base vectors
        k_unit = np.array([0, 0, 1], dtype=np.float64).reshape((3, 1))

        # Calculate the remaining quasi-dipole base vectors
        g1 = ((self.RE + np.float64(height))
              / (self.RE + self.refh)) ** (3 / 2) * d1 / F_scalar
        g2 = -1.0 / (2.0 * F_scalar * np.tan(np.radians(qlat))) * (
            k_unit + ((self.RE + np.float64(height))
                      / (self.RE + self.refh)) * d2 / cos_mag_inc)
        g3 = k_unit * F_scalar
        f3 = np.cross(g1.T, g2.T).T

        # Reshape the output
        out = [f1, f2, f3, g1, g2, g3, d1, d2, d3, e1, e2, e3]

        if np.any(alat == -9999):
            warnings.warn(''.join(['Base vectors g, d, e, and f3 set to NaN ',
                                   'where apex latitude is undefined (apex ',
                                   'height may be < reference height)']))
            mask = alat == -9999
            for i in range(len(out) - 2):
                out[i + 2] = np.where(mask, np.nan, out[i + 2])

        out = tuple(np.squeeze(bvec) for bvec in out)

        return out

    def get_apex(self, lat, height=None):
        """ Calculate apex height

        Parameters
        ----------
        lat : float
            Apex latitude in degrees
        height : float or NoneType
            Height above the surface of the Earth in km or NoneType to use
            reference height (default=None)

        Returns
        -------
        apex_height : float
            Height of the field line apex in km

        """
        # Check the latitude
        lat = helpers.checklat(lat, name='alat')

        # Ensure the height is set
        if height is None:
            height = self.refh

        # Calculate the apex height
        cos_lat_squared = np.cos(np.radians(lat)) ** 2
        apex_height = (self.RE + height) / cos_lat_squared - self.RE

        return apex_height

    def set_epoch(self, year):
        """Updates the epoch for all subsequent conversions.

        Parameters
        ----------
        year : float
            Decimal year

        """
        # Set the year and load the data file
        self.year = np.float64(year)
        fa.loadapxsh(self.datafile, self.year)

        # Call the Fortran routine to set time
        fa.cofrm(self.year, self.igrf_fn)
        return

    def set_refh(self, refh):
        """Updates the apex reference height for all subsequent conversions.

        Parameters
        ----------
        refh : float
            Apex reference height in km

        Notes
        -----
        The reference height is the height to which field lines will be mapped,
        and is only relevant for conversions involving apex (not quasi-dipole).

        """
        self.refh = refh

    def get_babs(self, glat, glon, height):
        """Returns the magnitude of the IGRF magnetic field in tesla.

        Parameters
        ----------
        glat : array_like
            Geodetic latitude in degrees
        glon : array_like
            Geodetic longitude in degrees
        height : array_like
            Altitude in km

        Returns
        -------
        babs : ndarray or float
            Magnitude of the IGRF magnetic field in Tesla

        """
        # Get the absolute value of the magnetic field at the desired location
        babs = self._get_babs(glat, glon, height)

        # If array is returned, the dtype is object, so convert to float
        return np.float64(babs)

    def bvectors_apex(self, lat, lon, height, coords='geo', precision=1e-10):
        """Returns the magnetic field vectors in apex coordinates.

        Parameters
        ----------
        lat : (N,) array_like or float
            Latitude
        lon : (N,) array_like or float
            Longitude
        height : (N,) array_like or float
            Altitude in km
        coords : {'geo', 'apex', 'qd'}, optional
            Input coordinate system
        precision : float, optional
            Precision of output (degrees) when converting to geo. A negative
            value of this argument produces a low-precision calculation of
            geodetic lat/lon based only on their spherical harmonic
            representation.
            A positive value causes the underlying Fortran routine to iterate
            until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces
            the input QD lat/lon to within the specified precision (all
            coordinates being converted to geo are converted to QD first and
            passed through APXG2Q).

        Returns
        -------
        main_mag_e3: (1, N) or (1,) ndarray
           IGRF magnitude divided by a scaling factor, D (d_scale) to give
           the main B field magnitude along the e3 base vector
        e3 : (3, N) or (3,) ndarray
           Base vector tangent to the contours of constant V_0 and Phi_A
        main_mag_d3: (1, N) or (1,) ndarray
           IGRF magnitude multiplied by a scaling factor, D (d_scale) to give
           the main B field magnitudee along the d3 base vector
        d3 : (3, N) or (3,) ndarray
           Base vector equivalent to the scaled main field unit vector

        Notes
        -----
        See Richmond, A. D. (1995) [4]_ equations 3.8-3.14

        The apex magnetic field vectors described by Richmond [1995] [4]_ and
        Emmert et al. [2010] [5]_, specfically the Be3 (main_mag_e3) and Bd3
        (main_mag_d3) components. The vector components are geodetic east,
        north, and up.

        References
        ----------
        Richmond, A. D. (1995) [4]_
        Emmert, J. T. et al. (2010) [5]_

        """
        # Convert the current coordinates to geodetic coordinates
        glat, glon = self.convert(lat, lon, coords, 'geo', height=height,
                                  precision=precision)

        # Get the magnitude of the magnetic field at the desired location
        babs = self.get_babs(glat, glon, height)

        # Retrieve the necessary base vectors
        _, _, _, _, _, _, d1, d2, d3, _, _, e3 = self.basevectors_apex(
            glat, glon, height, coords='geo')

        # Perform the calculations described in [4]
        d1_cross_d2 = np.cross(d1.T, d2.T).T
        d_scale = np.sqrt(np.sum(d1_cross_d2 ** 2, axis=0))  # D in [4] 3.13

        main_mag_e3 = babs / d_scale  # Solve for b0 in [4] 3.13
        main_mag_d3 = babs * d_scale  # Solve for b0 in [4] 3.10

        return main_mag_e3, e3, main_mag_d3, d3