# Source code for spacepy.ae9ap9

    #!/usr/bin/env python
    # -*- coding: utf-8 -*-
    """
    Module for reading and dealing with AE9AP9 data files.

    See https://www.vdl.afrl.af.mil/programs/ae9ap9/ to download the model.
    This is not a AE9AP9 runner.

    Authors: Brian Larsen, Steve Morley
    Institution: Los Alamos National Laboratory
    Contact: balarsen@lanl.gov

    Copyright 2015 Los Alamos National Security, LLC.

    This module provides a convenient class for handling data from
    AE9/AP9 (and legacy models provided by the software).

    .. rubric:: Class
    .. autosummary::
        :template: clean_class.rst
        :toctree: autosummary

        Ae9Data

    Though the class is derived from SpacePy's SpaceData, the class also provides several methods
    targeted at the AE9/AP9 output. Additional functions for working with the data are provided.

    .. rubric:: Functions
    .. autosummary::
        :template: clean_function.rst
        :toctree: autosummary

        readFile
        parseHeader

    """

    import datetime as dt
    import gzip
    import os
    import re
    import warnings

    from dateutil import relativedelta
    import numpy as np

    import spacepy.coordinates as spc
    import spacepy.datamodel as dm
    import spacepy.plot as splot
    import spacepy.time as spt
    import spacepy.toolbox as tb

    __contact__ = 'Brian Larsen, balarsen@lanl.gov'

    class Ae9Data(dm.SpaceData):
        """Dictionary-like container for AE9/AP9/SPM data, derived from SpacePy's datamodel

        To inspect the variables within this class, use the tree method.
        To export the data to a CDF, HDF5 or JSON-headed ASCII file use the relevant "to" method
        (toCDF, toHDF5, toJSONheadedASCII).

        .. autosummary::

            ~Ae9Data.getLm
            ~Ae9Data.plotOrbit
            ~Ae9Data.plotSummary
            ~Ae9Data.plotSpectrogram
            ~Ae9Data.setUnits

        .. automethod:: getLm
        .. automethod:: plotOrbit
        .. automethod:: plotSummary
        .. automethod:: plotSpectrogram
        .. automethod:: setUnits
        """

        def setUnits(self, per=None):
            """Set units of energy and flux/fluence

            If keyword 'per' is set to None, this method reports the units currently set.
            To set energy in MeV and flux/fluence in 'per MeV', set 'per=MeV'. Valid options are
            'eV', 'keV', 'Mev' and 'GeV'.

            Other Parameters
            ----------------
            per : string (optional)
                Energy units for both energy and flux/fluence
            """

            curr = self['Energy'].attrs['UNITS']
            particle_var = self.attrs['varname']
            if not per:
                print('Energy units: {0}'.format(curr))
                print('{0} units: {1}'.format(particle_var, self[particle_var].attrs['UNITS']))
                return
            unitlist = ['eV', 'keV', 'MeV', 'GeV']
            faclist = [1, 1e3, 1e6, 1e9]
            if per not in unitlist:
                raise ValueError("Units of {0} are not supported: Valid options are {1}".format(per, unitlist))
            else:
                if per == curr: return
                unitidx_from = unitlist.index(curr)
                unitidx_to = unitlist.index(per)
                self[particle_var] /= faclist[unitidx_from]  # convert to eV
                self[particle_var] *= faclist[unitidx_to]  # convert to target units
                self[particle_var].attrs['UNITS'] = self[particle_var].attrs['UNITS'].replace(curr, per)
                self['Energy'] *= faclist[unitidx_from]  # convert to eV
                self['Energy'] /= faclist[unitidx_to]  # convert to target units
                self['Energy'].attrs['UNITS'] = self['Energy'].attrs['UNITS'].replace(curr, per)

        def getLm(self, alpha=[90], model='T89'):
            """Calculate McIlwain L for the imported AE9/AP9 run and add to object
            """
            import spacepy.irbempy as ib
            ticks = spt.Ticktock(self['Epoch'])
            loci = spc.Coords(self['Coords'], self['Coords'].attrs['COORD_SYS'], 'car')
            loci.ticks = ticks
            retvals = ib.get_Lm(ticks, loci, alpha, extMag=model)
            self['Lm'] = dm.dmarray(retvals['Lm'].squeeze(), attrs={ 'MODEL': model })

        def _makeOrbitAxis(self, ser1, ser2, ax_target):
            """Helper function. Not intended for direct use.

            cx is the set of coordinates for plotting, ax_target is an axes object"""
            l1 = ax_target.plot(ser1, ser2)
            if np.abs((ax_target.get_xlim()[1] - ax_target.get_xlim()[0])) < 1:
                refpt = np.abs(np.max(ax_target.get_xlim()))
                ax_target.set_xlim([-1.25 * refpt, 1.25 * refpt])
                ax_target.set_ylim(ax_target.get_xlim())
                l1[0].set_marker('o')
                c1 = splot.plt.Circle([0, 0], radius=1.0, fc='none', ec='k')
                ax_target.add_artist(c1)
            else:
                splot.dual_half_circle(ax=ax_target)
            return ax_target

        def plotSummary(self, timerange=None, coord_sys=None, fig_target=None, spec=False, orbit_params=(False, True),
                        **kwargs):
            """Generate summary plot of AE9/AP9/SPM data loaded

            spec : if True, plot spectrogram instead of flux/fluence lineplot, requires 'ecol' keyword
            """
            if timerange:
                if isinstance(timerange, spt.Ticktock):
                    t1 = timerange.UTC[0]
                    t2 = timerange.UTC[-1]
                elif isinstance(timerange[0], dt.datetime):
                    t1 = timerange[0]
                    t2 = timerange[-1]
                else:
                    raise TypeError('Incorrect data type provided for timerange')
                # now select subset
                i_use = tb.tOverlapHalf([t1, t2], self['Epoch'])
                t_use = self['Epoch'][i_use]
                c_use = self['Coords'][i_use]
                f_use = self[self.attrs['varname']][i_use, ...]
            else:
                t_use = self['Epoch']
                c_use = self['Coords']
                f_use = self[self.attrs['varname']]

            if coord_sys and (coord_sys.upper() != c_use.attrs['COORD_SYS']):
                # TODO: We assume cartesian, make flexible so can take spherical
                cx = spc.Coords(c_use, c_use.attrs['COORD_SYS'], 'car')
                cx.ticks = spt.Ticktock(t_use)
                cx = cx.convert(coord_sys.upper(), 'car').data
            else:
                coord_sys = c_use.attrs['COORD_SYS']
                cx = c_use
            sys_subs = r'$_{' + coord_sys + r'}$'

            splot.style('spacepy')
            if orbit_params[0]:
                landscape = orbit_params[1]
                locs = [121, 122] if landscape else [211, 212]
            else:
                locs = [223, 224]
                landscape = True
            if fig_target:
                fig, ax1 = splot.set_target(fig_target, loc=locs[0])
            else:
                fig, ax1 = splot.set_target(fig_target, figsize=(8, 8), loc=locs[0])
            fig, ax2 = splot.set_target(fig, loc=locs[1])

            ax1 = self._makeOrbitAxis(cx[:, 0], cx[:, 1], ax1)
            ax2 = self._makeOrbitAxis(cx[:, 0], cx[:, 2], ax2)
            if landscape: ax1.set_xlabel('X{0} [{1}]'.format(sys_subs, self['Coords'].attrs['UNITS']))
            ax2.set_xlabel('X{0} [{1}]'.format(sys_subs, self['Coords'].attrs['UNITS']))
            ax1.set_ylabel('Y{0} [{1}]'.format(sys_subs, self['Coords'].attrs['UNITS']))
            ax2.set_ylabel('Z{0} [{1}]'.format(sys_subs, self['Coords'].attrs['UNITS']))
            ax1xl, ax1yl, ax2xl, ax2yl = ax1.get_xlim(), ax1.get_ylim(), ax2.get_xlim(), ax2.get_ylim()
            maxabslim = np.max(np.abs([ax1xl, ax1yl, ax2xl, ax2yl]))
            if np.abs((ax1.get_xlim()[1] - ax1.get_xlim()[0])) < 1:
                refpt = maxabslim
                ax1.set_xlim([-1.25 * refpt, 1.25 * refpt])
                ax1.set_ylim(ax1.get_xlim())
                refpt = ax2.get_xlim()[0]
                ax2.set_xlim([-1.25 * refpt, 1.25 * refpt])
                ax2.set_ylim(ax1.get_xlim())
            else:
                ax1.set_xlim([-maxabslim, maxabslim])
                ax1.set_ylim([-maxabslim, maxabslim])
                ax2.set_xlim([-maxabslim, maxabslim])
            ax2.set_ylim(ax2.get_xlim())
            ax1.invert_yaxis()
            ax1.invert_xaxis()
            ax2.invert_xaxis()
            ax1.set_aspect('equal')
            ax2.set_aspect('equal')

            if not orbit_params[0]:
                ax3 = splot.plt.subplot2grid((2, 2), (0, 0), colspan=2)
                if not spec:
                    l3 = ax3.semilogy(t_use, f_use)
                    ylab = '{0} ['.format(self.attrs['varname']) + re.sub('(\^[\d|-]*)+', _grp2mathmode,
                                                                          self[self.attrs['varname']].attrs['UNITS']) + ']'
                    ax3.set_ylabel(ylab)
                    for ll, nn in zip(l3, self['Energy']):
                        ll.set_label('{0} {1}'.format(nn, self['Energy'].attrs['UNITS']))
                    ncol = len(self['Energy']) // 2 if len(self['Energy']) <= 6 else 3
                    leg = ax3.legend(loc='center', bbox_to_anchor=(0.5, 1), ncol=ncol,
                                     frameon=True, framealpha=0.5)
                    lims3 = ax3.get_ylim()
                    newupper = 10 ** (np.log10(lims3[0]) + (np.log10(lims3[1] / lims3[0]) * 1.125))
                    ax3.set_ylim([lims3[0], newupper])
                    splot.applySmartTimeTicks(ax3, t_use)
                    fig.tight_layout()
                    pos3 = ax3.get_position()
                    ax3.set_position([pos3.x0, pos3.y0, pos3.width, pos3.height * 0.8])
                    # fig.suptitle('{model_type}\n'.format(**self.attrs) +
                    #             '{0} - {1}'.format(t_use[0].isoformat()[:19], t_use[-1].isoformat()[:19]))
                else:
                    if timerange: raise NotImplementedError('Time range selection not yet implemented for spectrograms')
                    pos3 = ax3.get_position()
                    ax3.set_position([pos3.x0, pos3.y0, pos3.width, pos3.height * 0.9])
                    ecol = kwargs['ecol'] if 'ecol' in kwargs else 0
                    ax3 = self.plotSpectrogram(target=ax3, ecol=ecol)
                    splot.plt.subplots_adjust(wspace=0.3)
                    pos1 = ax1.get_position()
                    ax1.set_position([pos3.x0, pos1.y0, pos1.width, pos1.height])

            splot.revert_style()
            return fig

        def plotOrbit(self, timerange=None, coord_sys=None, landscape=True, fig_target=None):
            """
            Plot X-Y and X-Z projections of satellite orbit in requested coordinate system
            """
            fig = self.plotSummary(timerange=timerange, coord_sys=coord_sys,
                                   fig_target=fig_target, orbit_params=(True, landscape))
            return fig

        def plotSpectrogram(self, ecol=0, pvars=None, **kwargs):
            '''
            Plot a spectrogram of the flux along the requested orbit, as a function of Lm and time

            Other Parameters
            ----------------
            zlim : list
                2-element list with upper and lower bounds for color scale
            colorbar_label : string
                text to appear next to colorbar (default is 'Flux' plus the units)
            ylabel : string
                text to label y-axis (default is 'Lm' plus the field model name)
            title : string
                text to appear above spectrogram (default is climatology model name, data type and energy)
            pvars : list
                list of plotting variable names in order [Epoch-like (X axis), Flux-like (Z axis), Energy (Index var for Flux-like)]
            ylim : list
                2-element list with upper and lower bounds for y axis
            '''
            import spacepy.plot as splot
            if 'Lm' not in self:
                self.getLm()
            sd = dm.SpaceData()

            if pvars is None:
                varname = self.attrs['varname']
                enname = 'Energy'
                epvar = 'Epoch'
            else:
                epvar = pvars[0]
                varname = pvars[1]
                enname = pvars[2]
            if len(self['Lm']) != len(self[epvar]): #Lm needs interpolating to new timebase
                import matplotlib.dates as mpd
                tmptimetarg, tmptimesrc = mpd.date2num(self[epvar]), mpd.date2num(self['Epoch'])
                sd['Lm'] = dm.dmarray(np.interp(tmptimetarg, tmptimesrc, self['Lm'], left=np.nan, right=np.nan))
            else:
                sd['Lm'] = self['Lm']
            #filter any bad Lm
            goodidx = sd['Lm']>1
            sd['Lm'] = sd['Lm'][goodidx]
            if 'ylim' in kwargs:
                Lm_lim = kwargs['ylim']
                del kwargs['ylim']
            else:
                Lm_lim = [2.0,8.0]
            #TODO: allow user-definition of bins in time and Lm
            sd['Epoch'] = dm.dmcopy(self[epvar])[goodidx] #TODO: assumes 1 pitch angle, generalize
            try:
                sd['1D_dataset'] = self[varname][goodidx,ecol] #TODO: assumes 1 pitch angle, generalize
            except IndexError: #1-D
                sd['1D_dataset'] = self[varname][goodidx]
            bins = [[],[]]
            if 'tbins' not in kwargs:
                bins[0] = spt.tickrange(self[epvar][0], self[epvar][-1], 3/24.).UTC
            else:
                bins[0] = kwargs['tbins']
                del kwargs['tbins']
            if 'Lbins' not in kwargs:
                bins[1] = np.arange(Lm_lim[0], Lm_lim[1], 1./4.)
            else:
                bins[1] = kwargs['Lbins']
                del kwargs['Lbins']
            spec = splot.spectrogram(sd, variables=['Epoch', 'Lm', '1D_dataset'], ylim=Lm_lim, bins=bins)

            if 'zlim' not in kwargs:
                zmax = 10 ** (int(np.log10(max(sd['1D_dataset']))) + 1)
                idx = np.logical_and(sd['1D_dataset'] > 0, sd['Lm'] > Lm_lim[0])
                idx = np.logical_and(idx, sd['Lm'] <= Lm_lim[1])
                zmin = 10 ** int(np.log10(min(sd['1D_dataset'][idx])))
                kwargs['zlim'] = [zmin, zmax]
            if 'colorbar_label' not in kwargs:
                flux_units = self[varname].attrs['UNITS']
                kwargs['colorbar_label'] = '{0} ['.format(varname) + re.sub('(\^[\d|-]*)+', _grp2mathmode, flux_units) + ']'
            if 'ylabel' not in kwargs:
                kwargs['ylabel'] = 'L$_M$' + ' ' + '[{0}]'.format(self['Lm'].attrs['MODEL'])
            if 'title' not in kwargs:
                kwargs['title'] = '{model_type} {varname}: '.format(**self.attrs) + \
                                  '{0:5.2f} {1}'.format(self[enname][ecol], self[enname].attrs['UNITS'])
            reset_shrink = splot.mpl.mathtext.SHRINK_FACTOR
            splot.mpl.mathtext.SHRINK_FACTOR = 0.85
            splot.mpl.mathtext.GROW_FACTOR = 1 / 0.85
            ax = spec.plot(cmap='plasma', **kwargs)
            splot.mpl.mathtext.SHRINK_FACTOR = reset_shrink
            splot.mpl.mathtext.GROW_FACTOR = 1 / reset_shrink
            return ax

    def readFile(fname, comments='#'):
        """
        read a model generated file into a datamodel.SpaceData object

        Parameters
        ==========
        fname : str
            filename of the file

        Other Parameters
        ================
        comments : str (optional)
            String that is the comments in the data file

        Returns
        =======
        out : :mod:`~spacepy.datamodel.SpaceData`
            Data contained in the  file

        Examples
        ========
        >>> from spacepy import ae9ap9
        >>> ae9ap9.readFile('ephem_sat.dat').tree(verbose=1)
        +
        |____Epoch (spacepy.datamodel.dmarray (121,))
        |____Coords (spacepy.datamodel.dmarray (121, 3))
        |____MJD (spacepy.datamodel.dmarray (121,))
        |____posComp (spacepy.datamodel.dmarray (3,))
        """
        if not os.path.isfile(fname):
            raise (ValueError("File {0} not found".format(fname)))
        # get the header information first
        header = parseHeader(fname)
        # and read in all the data
        data = np.loadtxt(fname, delimiter=header['delimiter'], comments=comments)
        # now parse the data
        ans = Ae9Data()
        # parse the datetime if it is there (it is always first)
        if 'datetime' in header['columns'][0]:
            if header['time_format'] == 'eDOY':  # have to massage the data first
                year = data[:, 0].astype(int)
                frac = data[:, 1]
                time = spt.Ticktock([dt.datetime(y, 1, 1) + relativedelta.relativedelta(days=v)
                                     for y, v in zip(year, frac)], 'UTC')
                ans[header['time_format']] = dm.dmarray(data[:, 0:2])
                ans[header['time_format']].attrs['VAR_TYPE'] = 'support_data'
                ans[header['time_format']].attrs['LABL_PTR_1'] = 'timeComp'
                ans['timeComp'] = dm.dmarray(['Year', 'eDOY'])
                ans['timeComp'].attrs['VAR_TYPE'] = 'metadata'
                del header['columns'][0:2]
                data = data[:, 2:]
            elif header['time_format'] == 'UTC':  # have to massage the data first
                tm = data[:, 0:6].astype(int)
                time = spt.Ticktock([dt.datetime(*v) for v in tm], 'UTC')
                ans[header['time_format']] = dm.dmarray(data[:, 0:6])
                ans[header['time_format']].attrs['VAR_TYPE'] = 'support_data'
                ans[header['time_format']].attrs['LABL_PTR_1'] = 'timeComp'
                ans['timeComp'] = dm.dmarray(['Year', 'Month', 'Day', 'Hour', 'Minute', 'Second'])
                ans['timeComp'].attrs['VAR_TYPE'] = 'metadata'
                del header['columns'][0:6]
                data = data[:, 6:]
            else:
                time = spt.Ticktock(data[:, 0], header['time_format'])
                ans[header['time_format']] = dm.dmarray(data[:, 0])
                ans[header['time_format']].attrs['VAR_TYPE'] = 'support_data'
                del header['columns'][0]
                data = data[:, 1:]

            ans['Epoch'] = dm.dmarray(time.UTC)
            ans['Epoch'].attrs['VAR_TYPE'] = 'support_data'
        # parse the position, it is always next
        if 'posx' in header['columns'][0]:
            varname = 'Coords'
            pos = dm.dmarray(data[:, 0:3])
            ans[varname] = pos
            ans[varname].attrs['UNITS'] = header['coord_system'][1]
            ans[varname].attrs['FIELDNAM'] = varname
            ans[varname].attrs['LABLAXIS'] = 'Position'
            ans[varname].attrs['DEPEND_0'] = 'Epoch'
            ans[varname].attrs['LABL_PTR_1'] = 'posComp'
            ans[varname].attrs['VAR_TYPE'] = 'data'
            ans[varname].attrs['COORD_SYS'] = header['coord_system'][0]
            ans['posComp'] = dm.dmarray(['X', 'Y', 'Z'])
            ans['posComp'].attrs['VAR_TYPE'] = 'metadata'
            del header['columns'][0:3]
            data = data[:, 3:]

        # if there are PA they are now
        if header['columns'] and ('pitch' in header['columns'][0]):
            raise (NotImplementedError("Sorry have not implemented pitch angle resolved files yet"))
        ## This is an issue and the files have repeadedlines at the same time for each pitch
        ##  angle. They would all need to be read in them combines into a multi-d array

        # now parse fluence, flux, or doserate
        # do all the rest of the column headers match?
        col_arr = np.asarray(header['columns'])
        if header['columns'] and (col_arr == col_arr[0]).all():
            varname = col_arr[0][0].title()
            ans[varname] = dm.dmarray(data)
            ans[varname].attrs['UNITS'] = str(col_arr[0][1])
            ans[varname].attrs['FIELDNAM'] = varname
            ans[varname].attrs['LABLAXIS'] = varname
            ans[varname].attrs['VAR_TYPE'] = 'data'
            ans[varname].attrs['SCALETYP'] = 'log'
            header['varname'] = varname
            ans[varname].attrs[
                'Description'] = '{flux_direction} {varname} based on {flux_type} from the {model_type} model'.format(
                **header)
            ans[varname].attrs['DEPEND_0'] = 'Epoch'
            if ans[varname].shape[1] == header['energy'][0].shape[0]:
                ans[varname].attrs['DEPEND_1'] = 'Energy'
            if 'percentile' in header:
                ans[varname].attrs['TITLE'] = '{0} percentile'.format(header['percentile'])

        # create the energy variable
        if 'energy' in header:
            ans['Energy'] = dm.dmarray(header['energy'][0])
            ans['Energy'].attrs['UNITS'] = header['energy'][1]
            ans['Energy'].attrs['FIELDNAM'] = 'Energy'
            ans['Energy'].attrs['LABLAXIS'] = 'Energy'
            ans['Energy'].attrs['VAR_TYPE'] = 'data'
            ans['Energy'].attrs['SCALETYP'] = 'log'
            ans['Energy'].attrs['VAR_TYPE'] = 'support_data'

        skip_header = ['columns', 'coord_system', 'delimiter', 'energy', 'flux_direction']
        for k in header:
            if k not in skip_header:
                ans.attrs[k] = header[k]
        return ans

    def _parseInfo(header):
        """
        given a header parse and return the common information in all headers

        # Time format:              Modified Julian Date
        # Coordinate system:        GEI (Geocentric Equatorial Inertial) Cartesian in Earth radii
        # Data Delimiter:           comma
        """
        # get the time format, it is in all the files
        ans = { }
        for ind, val in enumerate(header):
            if "Time format:" in val:
                if "Modified Julian Date" in val:
                    ans['time_format'] = 'MJD'
                elif "Year, day_of_year.frac" in val:
                    ans['time_format'] = 'eDOY'
                elif "Year, Month, Day, Hour, Minute, Seconds" in val:
                    ans['time_format'] = 'UTC'
                else:
                    raise (NotImplementedError("Sorry can't read that time format yet: {0}".format(val)))
            elif "Coordinate system:" in val:
                coord_sys = val.split(':')[1].strip().split()[0]
                if "in Earth radii" in val:
                    ans['coord_system'] = (coord_sys, 'Re')
                else:
                    ans['coord_system'] = (coord_sys, '')
            elif "Data Delimiter:" in val:
                if "comma" in val:
                    ans['delimiter'] = ','
                else:
                    raise (NotImplementedError("Sorry can't read that delimiter yet"))
            elif "Ae9Ap9 Software Version:" in val:
                match = re.match(r'^Ae9Ap9 Software Version:.*(\d\.\d\d\.\d\d\d)$', val)
                version = match.group(1).strip()
                ans['software_version'] = version
            elif "Model type:" in val:
                match = re.match(r'^Model type:(.*)$', val)
                mtype = match.group(1).strip()
                ans['model_type'] = mtype
            elif "CRRESELE Software Version:" in val:
                match = re.match(r'^CRRESELE Software Version:.*(\d\.\d\d\.\d\d\d)$', val)
                version = match.group(1).strip()
                ans['software_version'] = version
                ans['model_type'] = 'CRRESELE/PRO'
            elif "Flux type:" in val:
                match = re.match(r'^Flux type:(.*)$', val)
                ftype = match.group(1).strip()
                ans['flux_type'] = ftype
            elif "Flux direction:" in val:
                if '=' not in val:
                    match = re.match(r'^Flux direction:(.*)$', val)
                    fdir = match.group(1).strip()
                else:
                    pa = val.split('=')[1].strip()
                    fdir = np.asarray(pa.split())
                ans['flux_direction'] = fdir
            elif "Energy levels" in val:
                match = re.match(r'^Energy levels.*\((.*)\):(.*)$', val)
                ans['energy'] = (np.asarray(match.group(2).strip().split()).astype(float), match.group(1).strip())
            # Get the orbital element propagator, two versions based on AE9 model changes
            # new format data file
            elif "Propagator" in val:  # New format
                match = re.search(r'Propagator:\ (.*)$', val)
                ans['propagator'] = match.group(1).strip()
            elif "generated from specified elements" in val:  # In both old and new
                match = re.search(r'^generated from specified elements.*:\ (.*)$', val)
                if match:  # But old has propagator on this line; process and warn
                    warnings.warn(
                        "Support for orbit files from AE9AP9 model <1.5 is deprecated; please update to model 1.5 or later.",
                        DeprecationWarning)
                    ans['propagator'] = match.group(1).strip()
        return ans

    def _readHeader(fname):
        """
        read only the header from a Ae9Ap9 file
        """
        dat = []
        if fname.endswith('.gz'):
            try:
                fp = gzip.open(fname, 'rt')
            except ValueError: #Python 2.7 (Windows) compatibility
                fp = gzip.open(fname)
        else:
            fp = open(fname, 'rt')
        while True:
            tmp = fp.readline()
            if tmp[0] not in ('#', 35):
                break
            dat.append(tmp.strip())
        dat = [v[1:].strip() for v in dat]
        fp.close()
        return dat

    def parseHeader(fname):
        """
        given an AE9AP9 output test file parse the header and return the information in a
        dictionary

            .. versionchanged:: 0.3.0

            The underlying AE9AP9 model changed the ephem file format and this
            reader was updated to match. Reading the old format will issue
            DeprecationWarning.

        Parameters
        ==========
        fname : str
            filename of the file

        Returns
        =======
        out : dict
            Dictionary of the header information in the file
        """
        ans = { }
        ans['filename'] = os.path.abspath(fname)

        header = _readHeader(fname)
        # get the information
        ans.update(_parseInfo(header))

        columns = header[-1].split(ans['delimiter'])
        columns_unq = _unique_elements_order(columns)
        len_columns_unq = len(columns)

        # pull the units off the column headings
        cols = []
        while columns_unq:
            tmp = columns_unq.pop(0)
            umatch = re.match(r'(.*)\((.*)\).*', tmp)
            if umatch:
                cols.append((umatch.group(1), umatch.group(2)))
            else:
                cols.append(tmp)
        # now go through and see if any non-tuple entry has the same name as a tuple entry
        to_remove = []
        for ii, c1 in enumerate(cols):
            if not isinstance(c1, tuple):
                for jj, c2 in enumerate(cols[ii:]):
                    if c2[0] == c1:
                        to_remove.append(ii)
        if to_remove:
            for v in to_remove:
                del cols[v]
        # and now replicate the last entry enough times to make it the right length
        cols.extend([cols[-1]] * (len_columns_unq - len(cols)))

        ans['columns'] = cols

        if 'pctile' in os.path.basename(fname):  # capture the percentile
            match = re.match(r'.*_(\d\d).txt', os.path.basename(fname))
            ans['percentile'] = float(match.group(1))
        return ans

    def _unique_elements_order(seq):
        """
        see http://stackoverflow.com/questions/480214/how-do-you-remove-duplicates-from-a-list-in-python-whilst-preserving-order
        """
        seen = set()
        seen_add = seen.add
        return [x for x in seq if not (x in seen or seen_add(x))]

    def combinePercentiles(files, timeframe='all', verbose=True):
        """
        combine files at different percentiles into one file with the spectra at different
        percentiles for easy plotting and analysis

        NOTE: Requires pandas for any timeframe other than 'all'

        Parameters
        ==========
        files : str
            Filenames of the percentile files

        Other Parameters
        ================
        timeframe : str
            Timeframe to average the input spectra over (either 'all' or a pandas understoop resample() time
        verbose : boolean
            Print out information while reading the files

        Returns
        =======
        out : :mod:`~spacepy.datamodel.SpaceData`
            Combined output spectra of the file
        """
        if not files:
            raise (ValueError("Must input files"))
        data = { }
        for fname in files:
            if verbose: print("Reading: {0}".format(fname))
            tmp = readFile(fname)
            if 'percentile' not in tmp.attrs:
                raise (ValueError("File {0} does not have a percentile key".format(fname)))
            data[tmp.attrs['percentile']] = tmp

        varnames = ['Flux', 'Fluence', 'Dose', None]
        for varname in varnames:
            if varname in data[list(data.keys())[0]]:
                break  # this leaves it where it was found
        if varname is None:
            raise (ValueError("Did not understand file type could not find one of {0}".format(varnames)))

        # now that they are all read in, we collapse them down to different time bases
        if timeframe == 'all':
            # this is just a full collapse
            # find the var to collapse
            d_comp = np.asarray([])
            ps = sorted(data.keys())[::-1]
            for p in ps:
                if data[p][varname].attrs['DEPEND_0'] == 'Epoch':
                    d_comp = np.concatenate((d_comp, data[p][varname].mean(axis=0)))  # full time collapse
                else:
                    d_comp = np.concatenate((d_comp, data[p][varname]))  # full time collapse
            d_comp = d_comp.reshape(len(ps), -1).T
        else:
            try:
                import pandas as pd
            except ImportError:
                raise (ImportError("panads is required for timeframe other than 'all'"))
            raise (NotImplementedError("Timeframes other than 'all' not yet implemented"))
            # make a dataframe
            # resample to timeframe
            # join everything
            # meet up at same place for output
        # now that we have all the data prep it to a spacedata
        ans = dm.SpaceData()
        ans[varname] = dm.dmarray(d_comp)
        if timeframe == 'all':
            ans[varname].attrs['DEPEND_0'] = 'Energy'
            ans[varname].attrs['DEPEND_1'] = 'Percentile'
            ans[varname].attrs['DISPLAY_TYPE'] = 'time_series'
        for k in data[p][varname].attrs:
            if 'DEPEND' in k:
                continue
            else:
                ans[varname].attrs[k] = data[p][varname].attrs[k]
        ans['Energy'] = data[p]['Energy']
        ans['Percentile'] = dm.dmarray(ps)
        ans['Percentile'].attrs['FIELDNAM'] = 'Percentile'
        ans['Percentile'].attrs['LABLAXIS'] = 'Percentile'
        ans['Percentile'].attrs['VALIDMIN'] = 0.0
        ans['Percentile'].attrs['VALIDMAX'] = 100.0
        ans['Percentile'].attrs['SCALETYP'] = 'support_data'
        return ans

    def _getData(fnames):
        """
        helper routine for dm.toCDF and dm.toHDF5 and dm.toJSONheadedASCII since the data
        prep is all the same
        """
        if hasattr(fnames, 'strip'):  # it is a string
            return readFile(fnames)
        else:
            return combinePercentiles(fnames)

    def _grp2mathmode(matchobj):
        """Helper function to latex-ify the units for plotting"""
        unitstr = matchobj.group()
        unitstr = unitstr[0] + '{' + unitstr[1:]
        return '$' + unitstr + '}$'