How to use the pynbody.units function in pynbody

*verbose* (default=False): if True, prints out the diagnostics at
     each iteration. Useful to determine whether the centering is
     zeroing in on the wrong part of the simulation.

    """
    import os

    if r is None:
        # use rough estimate for a maximum radius
        # results will be insensitive to the exact value chosen
        r = (sim["x"].max() - sim["x"].min()) / 2

    elif isinstance(r, str) or issubclass(r.__class__, units.UnitBase):
        if isinstance(r, str):
            r = units.Unit(r)
        r = r.in_units(sim['pos'].units, **sim.conversion_context())

    mass = np.asarray(sim['mass'], dtype='double')
    pos = np.asarray(sim['pos'], dtype='double')

    com = _com.shrink_sphere_center(pos, mass, min_particles, shrink_factor, r)
    logger.info("Final SSC=%s", com)

    return array.SimArray(com, sim['pos'].units)

gv = gpro['rotation_curve_spherical'].in_units(v_units)
            dv = dpro['rotation_curve_spherical'].in_units(v_units)
            sv = spro['rotation_curve_spherical'].in_units(v_units)
        else:
            gv = gpro['v_circ'].in_units(v_units)
            dv = dpro['v_circ'].in_units(v_units)
            sv = spro['v_circ'].in_units(v_units)
        p.plot(r, gv, "--", label="gas")
        p.plot(r, dv, label="dark")
        p.plot(r, sv, linestyle="dotted", label="star")
    else:
        p.plot(r, v, **kwargs)

    if yrange:
        p.axis(
            [min_r, units.Unit(max_r).in_units(r.units), yrange[0], yrange[1]])

    if not axes:
        p.xlabel("r / $" + r.units.latex() + "$", fontsize='large')
        p.ylabel("v$_c / " + v.units.latex() + '$', fontsize='large')

    if legend:
        p.legend(loc=0)

    if (filename):
        logger.info("Saving %s", filename)
        p.savefig(filename)

    return r, v

all of `velocity`, `distance`, `mass` and `temperature`
        units. The units can be given as strings or as pynbody `Unit`
        objects.

        If any of the units are not specified and a previous
        `file_units_system` does not exist, the defaults are used.
        """
        from .. import config_parser
        import ConfigParser

        # if the units system doesn't exist (if this is a new snapshot), create
        # one
        if len(self._file_units_system) < 3:
            warnings.warn("Previous unit system incomplete -- using defaults")
            self._file_units_system = [
                units.Unit(x) for x in ('G', '1 kpc', '1e10 Msol')]

        else:
            # we want to change the base units -- so convert to original
            # units first and then set all arrays to new unit system
            self.original_units()


        # if any are missing, work them out from what we already have:

        if velocity is None:
            velocity = self.infer_original_units('km s^-1')

        if distance is None:
            distance = self.infer_original_units('kpc')

        if mass is None:

def __get_dtype_dims_and_units(self, fam, translated_name):
        if fam is None:
            fam = self.families()[0]

        inferred_units = units.NoUnit()
        representative_dset = None
        representative_hdf = None
        # not all arrays are present in all hdfs so need to loop
        # until we find one
        for hdf0 in self._hdf_files:
            try:
                representative_dset = self._get_hdf_dataset(hdf0[
                                                  self._family_to_group_map[fam][0]], translated_name)
                representative_hdf = hdf0
                if hasattr(representative_dset, "attrs"):
                    inferred_units = self._get_units_from_hdf_attr(representative_dset.attrs)

                if len(representative_dset)!=0:
                    # suitable for figuring out everything we need to know about this array
                    break
            except KeyError:

*mag_range*: float, float (default: None)
		 If provided, the brightest and faintest surface brightnesses in the range,
		 in mag arcsec^-2. Takes precedence over dynamic_range.

	   *with_dust*: bool, (default: False) 
		 If True, the image is rendered with a simple dust screening model
		 based on Calzetti's law.

	   *z_range*: float, (default: 50.0)
		 If with_dust is True this parameter specifies the z range
		 over which the column density will be calculated.
		 The default value is 50 kpc.

	'''

	if isinstance(width, str) or issubclass(width.__class__, _units.UnitBase):
		if isinstance(width, str):
			width = _units.Unit(width)
		width = width.in_units(sim['pos'].units, **sim.conversion_context())

	if starsize is not None:
		smf = filt.HighPass('smooth', str(starsize) + ' kpc')
		sim.s[smf]['smooth'] = array.SimArray(starsize, 'kpc', sim=sim)

	r = image(sim.s, qty=r_band + '_lum_den', width=width, log=False,
			  units="pc^-2", clear=False, noplot=True, resolution=resolution) * r_scale
	g = image(sim.s, qty=g_band + '_lum_den', width=width, log=False,
			  units="pc^-2", clear=False, noplot=True, resolution=resolution) * g_scale
	b = image(sim.s, qty=b_band + '_lum_den', width=width, log=False,
			  units="pc^-2", clear=False, noplot=True, resolution=resolution) * b_scale

	# convert all channels to mag arcsec^-2

try:
            f = open(self.filename + "." + array_name + ".pynbody-meta", 'r')
        except IOError:
            return self._default_units_for(array_name), None, None

        res = {}
        for l in f:

            X = l.split(":")

            if len(X) == 2:
                res[X[0].strip()] = X[1].strip()

        try:
            u = units.Unit(res['units'])
        except:
            u = None
        try:
            fams = [family.get_family(x.strip())
                    for x in res['families'].split(" ")]
        except:
            fams = None

        try:
            dtype = np.dtype(res['dtype'])
        except:
            dtype = None

        return u, fams, dtype

If not None, sets the maximum size of stars in the image

       *ret_im*: bool (default: False)
         if True, the NxNx3 image array is returned

       *plot*: bool (default: True)
         if True, the image is plotted

       *axes*: matplotlib axes object (deault: None)
         if not None, the axes object to plot to

       *dynamic_range*: float (default: 2.0)
         The number of dex in luminosity over which the image brightness ranges
    '''

    if isinstance(width, str) or issubclass(width.__class__, _units.UnitBase):
        if isinstance(width, str):
            width = _units.Unit(width)
        width = width.in_units(sim['pos'].units, **sim.conversion_context())

    if starsize is not None:
        smf = filt.HighPass('smooth', str(starsize) + ' kpc')
        sim.s[smf]['smooth'] = array.SimArray(starsize, 'kpc', sim=sim)

    r = image(sim.s, qty=r_band + '_lum_den', width=width, log=False,
              av_z=True, clear=False, noplot=True) * r_scale
    g = image(sim.s, qty=g_band + '_lum_den', width=width, log=False,
              av_z=True, clear=False, noplot=True) * g_scale
    b = image(sim.s, qty=b_band + '_lum_den', width=width, log=False,
              av_z=True, clear=False, noplot=True) * b_scale

    #r,g,b = nw_scale_rgb(r,g,b)

@RamsesSnap.decorator
def translate_info(sim):

    cosmo = 'aexp' in sim._info

    if sim._info['H0']>10:
        sim.properties['a'] = sim._info['aexp']
        sim.properties['omegaM0'] = sim._info['omega_m']
        sim.properties['omegaL0'] = sim._info['omega_l']
        sim.properties['h'] = sim._info['H0'] / 100.

    # N.B. these conversion factors provided by ramses already have the
    # correction from comoving to physical units
    d_unit = sim._info['unit_d'] * units.Unit("g cm^-3")
    t_unit = sim._info['unit_t'] * units.Unit("s")
    l_unit = sim._info['unit_l'] * units.Unit("cm")

    sim.properties['boxsize'] = sim._info['boxlen'] * l_unit

    if sim._info['omega_k'] == sim._info['omega_l'] == sim._info['omega_b'] == 0.0:
        sim.properties['time'] = sim._info['time'] * t_unit
    else:
        sim.properties['time'] = analysis.cosmology.age(
            sim) * units.Unit('Gyr')

    sim._file_units_system = [d_unit, t_unit, l_unit]

def _get_units_from_hdf_attr(self, hdfattrs) :
        """Return the units based on HDF attributes VarDescription"""



        VarDescription = str(hdfattrs['VarDescription'])
        CGSConversionFactor = float(hdfattrs['CGSConversionFactor'])
        aexp = hdfattrs['aexp-scale-exponent']
        hexp = hdfattrs['h-scale-exponent']

        arr_units = self._get_units_from_description(VarDescription, CGSConversionFactor)

        if not np.allclose(aexp, 0.0):
            arr_units *= (units.a)**util.fractions.Fraction.from_float(float(aexp)).limit_denominator()
        if not np.allclose(hexp, 0.0):    
            arr_units *= (units.h)**util.fractions.Fraction.from_float(float(hexp)).limit_denominator()
  
        return arr_units

def _load_array(self, array_name, fam=None):
        if fam is None:
            self._attempt_load_all_families(array_name)
            return

        f = self._open_file_for_array(fam, array_name)

        _, nbod, ndim, dtype = self._load_header(f)
        if array_name not in self.keys():
            self._create_family_array(array_name, fam, ndim=ndim,dtype=dtype)
        r = self[fam][array_name]
        if units.has_units(r):
            r.convert_units(self._default_units_for(array_name))
        else:
            r.set_default_units()

        disk_dtype = dtype.newbyteorder('>')

        buf_shape = _max_buf
        if ndim > 1:
            buf_shape = (_max_buf, ndim)

        b = np.empty(buf_shape)

        # skip over min and max values (see issue #211)
        np.fromfile(f, dtype=disk_dtype, count=2 * ndim)

        for readlen, buf_index, mem_index in self._load_control.iterate(fam, fam):