def cartesian_to_keplerian(pos, vel, kep):
    """Given observables XYZ position and velocity, compute orbital elements.  Position must be in
       au and velocity in km/s.  Returns astropy unit objects for all orbital elements.
        Inputs:
            pos (3xN arr) [au]: position in xyz coords in au, with 
                        x = pos[0], y = pos[1], z = pos[2] for each of N orbits
                        +x = +Dec, +y = +RA, +z = towards observer
            vel (3xN arr) [km/s]: velocity in xyz plane.
        Returns:
            sma (1xN arr flt) [au]: semi-major axis in au, must be an astropy units object
            ecc (1xN arr flt) [unitless]: eccentricity
            inc (1xN arr flt) [deg]: inclination
            argp (1xN arr flt) [deg]: argument of periastron
            lon (1xN arr flt) [deg]: longitude of ascending node
            meananom (1xN arr flt) [radians]: mean anomaly 
        Written by Logan Pearce, 2019, inspired by Sarah Blunt
    """
    import numpy as np
    import astropy.units as u
    rcrossv = np.cross(pos, vel) * u.au * (u.km / u.s)
    h = np.sqrt(rcrossv[0] ** 2 + rcrossv[1] ** 2 + rcrossv[2] ** 2)
    n = rcrossv / h
    inc = np.arccos(n[2])
    lon = np.arctan2(n[0], -n[1])
    r = np.sqrt(pos[0] ** 2 + pos[1] ** 2 + pos[2] ** 2)
    v = np.sqrt(vel[0] ** 2 + vel[1] ** 2 + vel[2] ** 2)
    sma = 1 / (2.0 / r - v ** 2 / kep)
    rdotv = pos[0] * vel[0] + pos[1] * vel[1] + pos[2] * vel[2]
    rdot = rdotv / r
    parameter = h ** 2 / kep
    ecosf = parameter / r - 1
    esinf = h * rdot / kep.to(u.m ** 3 / u.s ** 2)
    ecc = np.sqrt(ecosf ** 2 + esinf ** 2)
    f = np.arctan2(esinf, ecosf)
    f = f.value % (2.0 * np.pi)
    E = 2.0 * np.arctan(np.sqrt((1 - ecc.value) / (1 + ecc.value)) * np.tan
        (f / 2.0))
    M = E - ecc * np.sin(E)
    rcosu = pos[0] * np.cos(lon) + pos[1] * np.sin(lon)
    rsinu = (-pos[0] * np.sin(lon) + pos[1] * np.cos(lon)) / np.cos(inc)
    uangle = np.arctan2(rsinu, rcosu)
    argp = uangle.value - f
    return sma.to(u.au), ecc, np.degrees(inc), np.degrees(argp
        ) % 360.0 * u.deg, np.degrees(lon.value) % 360.0 * u.deg, M
