def main():
    p = ArgumentParser(description='Makes unit flux eV^-1 as input to GLOW or Transcar to create ionospheric eigenprofiles')
    p.add_argument('-i', '--inputgridfn', help='original Zettergren input flux grid to base off of', default='zettflux.csv')
    p.add_argument('-o', '--outfn', help='hdf5 file to write with ionospheric response (eigenprofiles)')
    p.add_argument('-t', '--simtime', help='yyyy-mm-ddTHH:MM:SSZ time of sim', nargs='+', default=['1999-12-21T00:00:00Z'])
    p.add_argument('-c', '--latlon', help='geodetic latitude/longitude (deg)', type=float, nargs=2, default=[65, -148.])
#    p.add_argument('-m', '--makeplot', help='show to show plots, png to save pngs of plots', nargs='+', default=['show'])
    p.add_argument('-M', '--model', help='specify auroral model (glow,rees,transcar)', default='glow')
    p.add_argument('-z', '--zlim', help='minimum,maximum altitude [km] to plot', nargs=2, default=(None, None), type=float)
    p.add_argument('--isotropic', help='(rees model only) isotropic or non-isotropic pitch angle', action='store_true')
    p.add_argument('--vlim', help='plotting limits on energy dep and production plots', nargs=2, type=float, default=(1e-7, 1e1))

    p = p.parse_args()

    if not p.outfn:
        print('you have not specified an output file with -o options, so I will only plot and not save result')

#    makeplot = p.makeplot

    if len(p.simtime) == 1:
        T = [parse(p.simtime[0])]
    elif len(p.simtime) == 2:
        T = list(rrule.rrule(rrule.HOURLY,
                             dtstart=parse(p.simtime[0]),
                             until=parse(p.simtime[1])))
# %% input unit flux
    Egrid = loadregress(Path(p.inputgridfn).expanduser())
    Ebins = makebin(Egrid)[:3]
    EKpcolor, EK, diffnumflux = ekpcolor(Ebins)
# %% ionospheric response
    """ three output eigenprofiles
    1) ver (optical emissions) 4-D array: time x energy x altitude x wavelength
    2) prates (production) 4-D array:     time x energy x altitude x reaction
    3) lrates (loss) 4-D array:           time x energy x altitude x reaction
    """
    model = p.model.lower()
    glat, glon = p.latlon

    if model == 'glow':
        ver, photIon, isr, phitop, zceta, sza, prates, lrates, tezs, sion = makeeigen(EK, diffnumflux, T, p.latlon,
                                                                                      p.makeplot, p.outfn, p.zlim)

        writeeigen(p.outfn, EKpcolor, T, ver.z_km, diffnumflux, ver, prates, lrates, tezs, p.latlon)
        # %% plots
        # input
        doplot(p.inputgridfn, Ebins)

        # output
        sim = namedtuple('sim', ['reacreq', 'opticalfilter'])
        sim.reacreq = sim.opticalfilter = ''

        for t in ver:  # TODO for each time
            # VER eigenprofiles, summed over wavelength
            ploteigver(EKpcolor, ver.z_km, ver.sum('wavelength_nm'),
                       (None,)*6, sim,
                       '{} Vol. Emis. Rate '.format(t))
            # volume production rate, summed over reaction
            plotprodloss(prates.loc[:, 'final', ...].sum('reaction'),
                         lrates.loc[:, 'final', ...].sum('reaction'),
                         t, glat, glon, p.zlim)
            # energy deposition
            plotenerdep(tezs, t, glat, glon, p.zlim)

    elif model == 'rees':
        assert len(T) == 1, 'only one time with rees for now.'
        z = glowalt()
        q = reesiono(T, z, Ebins.loc[:, 'low'], glat, glon, p.isotropic)

        writeeigen(p.outfn, Ebins, T, z, prates=q, tezs=None, latlon=(glat, glon))

        plotA(q, 'Volume Production Rate {}  {} {}'.format(T, glat, glon), p.vlim)
    elif model == 'transcar':
        raise NotImplementedError('Transcar by request')
    else:
        raise NotImplementedError('I am not yet able to handle your model {}'.format(model))
# %% plots

    show()