How to use the pynbody.units.Unit function in pynbody

def SolarNeighborhood(r1=units.Unit("5 kpc"), r2=units.Unit("10 kpc"), height=units.Unit("2 kpc"), cen=(0, 0, 0)):
    """
    Convenience function that returns a filter which selects particles
    in a disc between radii `r1` and `r2` and thickness `height`.
    """

    x = Disc(max(r1, r2), height, cen) & ~Disc(min(r1, r2), height, cen)
    x._descriptor = "Solar Neighborhood"
    return x

"""

    if center:
        angmom.faceon(sim)

    if min:
        min_r = min
    else:
        min_r = sim['rxy'].min()
    if max:
        max_r = max
    else:
        max_r = sim['rxy'].max()

    if isinstance(pretime, str):
        pretime = units.Unit(pretime)

    diskstars = sim.star[filt.Disc(max_r, disk_height)]
    youngstars = np.where(diskstars['tform'].in_units("Myr") >
                          sim.properties['time'].in_units(
                              "Myr", **sim.conversion_context())
                          - pretime.in_units('Myr'))[0]

    pro = profile.Profile(diskstars[youngstars], type=bin_spacing,
                          nbins=nbins, min=min_r, max=max_r)

    r = pro['rbins'].in_units(r_units)
    fourierprof = pro['fourier']
    a2 = fourierprof['amp'][2]

    if clear:
        p.clf()

width = _units.Unit(width)
        width = width.in_units(sim['pos'].units, **sim.conversion_context())

    width = float(width)

    pixel_size = width / vector_resolution
    X, Y = np.meshgrid(np.arange(-width / 2 + pixel_size/2, width / 2 + pixel_size/2, pixel_size ),
                       np.arange(-width / 2 + pixel_size/2, width / 2 + pixel_size/2, pixel_size))

    im = image(sim, width=width, **kwargs)

    if mode == 'quiver':
        if scale is None:
            Q = p.quiver(X, Y, vx, vy, color=vector_color, edgecolor=edgecolor)
        else:
            Q = p.quiver(X, Y, vx, vy, scale=_units.Unit(scale).in_units(
                sim['vel'].units), color=vector_color, edgecolor=edgecolor)
        if quiverkey:
        	qk = p.quiverkey(Q, key_x, key_y, key_unit.in_units(sim['vel'].units, **sim.conversion_context()),
                    r"$\mathbf{" + key_unit.latex() + "}$", labelcolor=key_color, color=key_color, fontproperties={'size': 24})
        if  quiverkey_bg_color is not None:
            qk.text.set_backgroundcolor(quiverkey_bg_color)
    elif mode == 'stream' :
        Q = p.streamplot(X, Y, vx, vy, color=vector_color, density=density)

	# RS - if a axis object is passed in, use the right limit call
    if subplot:
        p.set_xlim(-width/2, width/2)
        p.set_ylim(-width/2, width/2)
    else:
        p.xlim(-width/2, width/2)
        p.ylim(-width/2, width/2)

simulation is in comoving units and you specify a different
    redshift z. Specifically, the physical density for the redshift
    you specify is calulated, but expressed as a comoving density *at
    the redshift of the snapshot*. This is intentional behaviour."""

    if z is None:
        z = f.properties['z']

    if unit is None:
        try:
            unit = f.dm["mass"].units / f.dm["pos"].units ** 3
        except units.UnitsException:
            unit = units.NoUnit()

    if hasattr(unit, "_no_unit"):
        unit = units.Unit("Msol kpc^-3 a^-3")

    omM = f.properties['omegaM0']
    omL = f.properties['omegaL0']
    h0 = f.properties['h']
    a = 1.0 / (1.0 + z)

    H_z = _a_dot(a, h0, omM, omL) / a
    H_z = units.Unit("100 km s^-1 Mpc^-1") * H_z

    rho_crit = (3 * H_z ** 2) / (8 * math.pi * units.G)

    return rho_crit.ratio(unit, **f.conversion_context())

def _get_units_from_description(self, description, expectedCgsConversionFactor=None):
        arr_units = units.Unit('1.0')
        conversion = 1.0
        for unitname in self._hdf_unitvar.keys():
            power = 1.
            if unitname in description:
                sstart = description.find(unitname)
                if sstart > 0:
                    if description[sstart - 1] == "/":
                        power *= -1.
                if len(description) > sstart + len(unitname):
                    # Just check we're not at the end of the line
                    if description[sstart + len(unitname)] == '^':
                        ## Has an index, check if this is negative
                        if description[sstart + len(unitname) + 1] == "-":
                            power *= -1.
                            power *= float(
                                description[sstart + len(unitname) + 2:-1].split()[0])  ## Search for the power

self.max = units.Unit(kwargs['max']).ratio(x.units,
                                                               **sim.conversion_context())
                else:
                    self.max = kwargs['max']
            else:
                self.max = np.max(x)
            if kwargs.has_key('bins'):
                self.nbins = len(kwargs['bins']) - 1
            elif kwargs.has_key('nbins'):
                self.nbins = kwargs['nbins']
            else:
                self.nbins = 100

            if kwargs.has_key('min'):
                if isinstance(kwargs['min'], str):
                    self.min = units.Unit(kwargs['min']).ratio(x.units,
                                                               **sim.conversion_context())
                else:
                    self.min = kwargs['min']
            else:
                self.min = np.min(x[x > 0])

            if kwargs.has_key('bins'):
                self._properties['bin_edges'] = kwargs['bins']
                self.min = kwargs['bins'].min()
                self.max = kwargs['bins'].max()
            elif type == 'log':
                self._properties['bin_edges'] = np.logspace(
                    np.log10(self.min), np.log10(self.max), num=self.nbins + 1)
            elif type == 'lin':
                self._properties['bin_edges'] = np.linspace(
                    self.min, self.max, num=self.nbins + 1)

@Profile.profile_property
def sb(self,band='v') :
    # At 10 pc (distance for absolute magnitudes), 1 arcsec is 10 AU=1/2.06e4 pc
    # In [5]: (np.tan(np.pi/180/3600)*10.0)**2
    # Out[5]: 2.3504430539466191e-09
    # 1 square arcsecond is thus 2.35e-9 pc^2
    sqarcsec_in_bin = self._binsize.in_units('pc^2') / 2.3504430539466191e-09
    bin_luminosity = 10.0**(-0.4*self['magnitudes,'+band])
    #import pdb; pdb.set_trace()
    surfb = -2.5*np.log10(bin_luminosity / sqarcsec_in_bin)
    surfb = array.SimArray(surfb, units.Unit('1'))
    surfb.sim = self.sim
    return surfb

it is used in place of the redshift in output f."""

    if z is None:
        z = f.properties['z']
    a = 1. / (1. + z)
    omegam0 = f.properties['omegaM0']
    omegal0 = f.properties['omegaL0']

    b, X = _lingrowthfac(z, omegam0, omegal0, return_norm=True)
    I = _lingrowthintegrand(a, omegam0)

    term1 = -(1.5 * omegam0 * a ** -3) * b / \
        math.sqrt(1. - omegam0 + omegam0 * a ** -3)
    term2 = (2.5 * omegam0) * hzoverh0(a, omegam0) ** 2 * a * I / X

    res = units.h * (term1 + term2) * 100. * units.Unit("km s^-1 Mpc^-1")

    return res.in_units(unit, **f.conversion_context())