How to use the pynbody.array.SimArray function in pynbody

@Profile.profile_property
def magnitudes(self, band='v'):
    """
    Calculate magnitudes in each bin
    """
    from . import luminosity

    magnitudes = np.zeros(self.nbins)
    for i in range(self.nbins):
        magnitudes[i] = luminosity.halo_mag(
            self.sim[self.binind[i]], band=band)
    magnitudes = array.SimArray(magnitudes, units.Unit('1'))
    magnitudes.sim = self.sim
    return magnitudes

@Profile.profile_property
def mass(self):
    """
    Calculate mass in each bin
    """
    
    if pynbody.config['verbose'] : print 'Profile: mass()'
    mass = array.SimArray(np.zeros(self.nbins), self.sim['mass'].units)

    with self.sim.immediate_mode : 
        pmass = self.sim['mass'].view(np.ndarray)

    for i in range(self.nbins):
        mass[i] = (pmass[self.binind[i]]).sum()

    mass.sim = self.sim

    return mass

for r in rxy_points:
        # Do four samples like Tipsy does
        r_acc_r = []
        for pos in [(r, 0, 0), (0, r, 0), (-r, 0, 0), (0, -r, 0)]:
            acc = accel(f, pos, eps)
            r_acc_r.append(np.dot(acc, pos))

        vel2 = np.mean(r_acc_r)
        if vel2 > 0:
            vel = math.sqrt(vel2)
        else:
            vel = 0

        vels.append(vel)

    x = array.SimArray(vels, units=u_out)
    x.sim = f.ancestor
    return x

@units.takes_arg_in_units((0, "Mpc h^-1"))
def correlation(r, powspec=PowerSpectrumCAMB):

    if hasattr(r, '__len__'):
        ax = pynbody.array.SimArray([correlation(ri,  powspec) for ri in r])
        ax.units = powspec.Pk.units * powspec.k.units ** 3
        return ax

    # Because sin kr becomes so highly oscilliatory, normal
    # quadrature is slow/inaccurate for this problem. The following
    # is the best way I could come up with to overcome that.
    #
    # For small kr, sin kr/kr is represented as a Taylor expansion and
    # each segment of the power spectrum is integrated over, summing
    # over the Taylor series to convergence.
    #
    # When the convergence of this starts to fail, each segment of the
    # power spectrum is still represented by a power law, but the
    # exact integral boils down to a normal incomplete gamma function
    # extended into the complex plane.
    #

def variance(M, f_filter=TophatFilter, powspec=PowerSpectrumCAMB, arg_is_R=False):
    if hasattr(M, '__len__'):
        ax = pynbody.array.SimArray(
            [variance(Mi, f_filter, powspec, arg_is_R) for Mi in M])
        # hopefully dimensionless
        ax.units = powspec.Pk.units * powspec.k.units ** 3
        return ax

    if arg_is_R:
        R = M
    else:
        R = f_filter.M_to_R(M)

    integrand = lambda k: k ** 2 * powspec(k) * f_filter.Wk(k * R) ** 2
    integrand_ln_k = lambda k: np.exp(k) * integrand(np.exp(k))
    v = scipy.integrate.romberg(integrand_ln_k, math.log(powspec.min_k), math.log(
        1. / R) + 3, divmax=10, rtol=1.e-4) / (2 * math.pi ** 2)

    return v

def __setitem__(self, name, item):
        """Set the contents of an array in this snapshot"""
        if self.is_derived_array(name) and not self.auto_propagate_off:
            raise RuntimeError("Derived array is not writable")

        if isinstance(name, tuple) or isinstance(name, list):
            index = name[1]
            name = name[0]
        else:
            index = None

        self._assert_not_family_array(name)

        if isinstance(item, array.SimArray):
            ax = item
        else:
            ax = np.asanyarray(item).view(array.SimArray)

        if name not in self.keys():
            # Array needs to be created. We do this through the
            # private _create_array method, so that if we are operating
            # within a particle-specific subview we automatically create
            # a particle-specific array
            try:
                ndim = len(ax[0])
            except TypeError:
                ndim = 1
            except IndexError:
                ndim = ax.shape[-1] if len(ax.shape) > 1 else 1

def v_mean(self):
    """SPH-smoothed mean velocity"""
    import sph

    sph.build_tree(self)

    nsmooth = config['sph']['smooth-particles']

    logger.info(
        'Calculating mean velocity with %d nearest neighbours' % nsmooth)

    sm = array.SimArray(np.empty_like(self['vel']),
                        self['vel'].units)

    self.kdtree.set_array_ref('rho',self['rho'])
    self.kdtree.set_array_ref('smooth',self['smooth'])
    self.kdtree.set_array_ref('mass',self['mass'])
    self.kdtree.set_array_ref('qty',self['vel'])
    self.kdtree.set_array_ref('qty_sm',sm)

    start = time.time()
    self.kdtree.populate('qty_mean',nsmooth)
    end = time.time()

    logger.info('Mean velocity done in %5.3g s' % (end - start))

    return sm

def test_functionals_nfw():

    r = SimArray(np.linspace(10, 500, 100), units="kpc")

    rho_from_static = NFW1.profile_functional_static(r, rhos, rs)
    rho_from_instance = NFW1.profile_functional(r)

    truth = SimArray([  2.50000000e+07,   1.07460773e+07,   5.62154816e+06,
            3.31384573e+06,   2.11880566e+06,   1.43727440e+06,
            1.01995927e+06,   7.50047528e+05,   5.67701535e+05,
            4.40058190e+05,   3.48027380e+05,   2.79994777e+05,
            2.28614491e+05,   1.89084016e+05,   1.58172652e+05,
            1.33652241e+05,   1.13951988e+05,   9.79427606e+04,
            8.47986748e+04,   7.39061323e+04,   6.48027682e+04,
            5.71356745e+04,   5.06323136e+04,   4.50799429e+04,
            4.03108479e+04,   3.61916013e+04,   3.26151535e+04,
            2.94949429e+04,   2.67604620e+04,   2.43538877e+04,
            2.22274964e+04,   2.03416640e+04,   1.86633064e+04,
            1.71646535e+04,   1.58222797e+04,   1.46163307e+04,
            1.35299042e+04,   1.25485503e+04,   1.16598657e+04,
            1.08531640e+04,   1.01192041e+04,   9.44996741e+03,
            8.83847317e+03,   8.27862511e+03,   7.76508333e+03,
            7.29315708e+03,   6.85871449e+03,   6.45810642e+03,
            6.08810188e+03,   5.74583322e+03,   5.42874931e+03,

loadblock = lambda count: np.fromstring(
                    f.read(count * 4), dtype=dtype, count=count).byteswap()
                # data = np.fromstring(f.read(3*len(self)*4),dtype).byteswap()
            else:
                loadblock = lambda count: np.fromstring(
                    f.read(count * 4), dtype=dtype, count=count)
                # data = np.fromstring(f.read(3*len(self)*4),dtype)

        ndim = 1

        self.ancestor._tipsy_arrays_binary = binary

        all_fam = [family.dm, family.gas, family.star]
        if fam is None:
            fam = all_fam
            r = np.empty(len(self), dtype=dtype).view(array.SimArray)
        else:
            r = np.empty(len(self[fam]), dtype=dtype).view(array.SimArray)

        for readlen, buf_index, mem_index in self._load_control.iterate(all_fam, fam):
            buf = loadblock(readlen)
            if mem_index is not None:
                r[mem_index] = buf[buf_index]


        if units is not None:
            r.units = units

        return r

def __init__(self, group_id, *args) :
        super(SubFindFOFGroup,self).__init__(group_id, *args)

        self._subhalo_catalogue = SubFindHDFSubhaloCatalogue(group_id, self._halo_catalogue)

        self._descriptor = "fof_group_"+str(group_id)

        # load properties
        for key in self._halo_catalogue._fof_properties.keys() :
            self.properties[key] = array.SimArray(self._halo_catalogue._fof_properties[key][group_id],
                                            self._halo_catalogue._fof_properties[key].units)
            self.properties[key].sim = self.base