How to use the pynbody.analysis.profile.Profile function in pynbody

def j_phi(self): 
    """
    Angle that the angular momentum vector of the bin makes with the x-axis in the xy plane.
    """
    j_phi = np.zeros(self.nbins)

    for i in range(self.nbins) : 
        subs = self.sim[self.binind[i]]
        jx = (subs['j'][:,0]*subs['mass']).sum()/self['mass'][i]
        jy = (subs['j'][:,1]*subs['mass']).sum()/self['mass'][i]
        j_phi[i] = np.arctan(jy,jx)
    
    return j_phi
    

class InclinedProfile(Profile) :
    """
    
    Creates a profile object to be used with a snapshot inclined by
    some known angle to the xy plane.

    In addition to the SimSnap object, it also requires the angle to
    initialize.

    **Example:**

    >>> s = pynbody.load('sim')
    >>> pynbody.analysis.angmom.faceon(s)
    >>> s.rotate_x(60)
    >>> p = pynbody.profile.InclinedProfile(s, 60)

    """

**needs a description of all keywords**

    By default, sbprof will use the formation mass of the star.
    In tipsy, this will be taken from the starlog file.

    '''

    if center:
        logger.info("Centering...")
        angmom.faceon(sim)

    logger.info("Selecting disk stars")
    diskstars = sim.star[filt.Disc(rmax, diskheight)]
    logger.info("Creating profile")
    ps = profile.Profile(diskstars, type=binning)
    logger.info("Plotting")
    r = ps['rbins'].in_units('kpc')

    if axes:
        plt = axes
    else:
        import matplotlib.pyplot as plt
    if clear:
        plt.clf()

    plt.plot(r, ps['sb,' + band], linewidth=2, **kwargs)
    if axes:
        plt.set_ylim(max(ps['sb,' + band]), min(ps['sb,' + band]))
    else:
        plt.ylim(max(ps['sb,' + band]), min(ps['sb,' + band]))
    if fit_exp:

mtwokpcstar=np.sum([h[i][twokpcf].s['mass'].in_units('Msol')]),
mhalogas=np.sum([h[i][notdiskf].g['mass'].in_units('Msol')]),
# 3D density profile
rhoprof = profile.Profile(h[i].dm,ndim=3,type='log')
gashaloprof = profile.Profile(h[i].g,ndim=3,type='log')
spheregasprof = profile.Profile(s.g,ndim=3,type='equaln',max=2*rvir,bins=200)
spherecoldgasprof = profile.Profile(s.g[cold],ndim=3,type='equaln',max=2*rvir,bins=200)
spherehotgasprof = profile.Profile(s.g[hot],ndim=3,type='equaln',max=2*rvir,bins=200)
spherestarprof = profile.Profile(s.s,ndim=3,type='equaln',max=2*rvir,bins=200)
spheredarkprof = profile.Profile(s.dm,ndim=3,type='equaln',max=2*rvir,bins=200)
# Rotation curve
rcpro = profile.Profile(h[i], type='equaln', nbins = 200, max = '40 kpc')
# surface brightness profile
diskstars = h[i].star[diskf]
sbprof = profile.Profile(diskstars,type='equaln',bins=200)
surfcoldgasprof = profile.Profile(h[i].g[diskf][cold],type='equaln',bins=200)
surfhotgasprof = profile.Profile(h[i].g[diskf][hot],type='equaln',bins=200)
surfgasprof = profile.Profile(h[i].g[diskf],type='equaln',bins=200)
# Kinematic decomposition
#decompprof = pynbody.analysis.decomp(h[i])
#dec = h[i].star['decomp']
sfh,sfhtimes = find_sfh(h[i],bins=100)
hrsfh,hrsfhtimes = find_sfh(h[i],bins=600)
sfr = np.sum(h[i].star[fifmyrf]['mass'].in_units('Msol')) / 1.5e7

# where did halo come from?  SN or virial shocks
# look at all particles inside 2.5 r_vir
twopfivef = pynbody.filt.Sphere(2.5*rvir)
stfrvir = s[twopfivef]
halogas = stfrvir[notdiskf].g[notcooldenf]
mtfhalogas = np.sum(halogas['mass'].in_units('Msol'))
iposcoolontime = halogas['coolontime'].in_units('yr')>0

import pynbody
import matplotlib.pylab as plt
from os import environ
# load the snapshot and set to physical units
s = pynbody.load('testdata/g15784.lr.01024.gz')

# load the halos
h = s.halos()

# center on the largest halo and align the disk
pynbody.analysis.angmom.faceon(h[1])

# create a profile object for the stars
s.physical_units()
p = pynbody.analysis.profile.Profile(h[1],min=.01,max=250,type='log',ndim=3)
pg = pynbody.analysis.profile.Profile(h[1].g,min=.01,max=250,type='log',ndim=3)
ps = pynbody.analysis.profile.Profile(h[1].s,min=.01,max=250,type='log',ndim=3)
pd = pynbody.analysis.profile.Profile(h[1].d,min=.01,max=250,type='log',ndim=3)

# make the plot
plt.plot(p['rbins'],p['v_circ'],label='total')
plt.plot(pg['rbins'],pg['v_circ'],label='gas')
plt.plot(ps['rbins'],ps['v_circ'],label='stars')
plt.plot(pd['rbins'],pd['v_circ'],label='dm')

plt.xlabel('R [kpc]')
plt.ylabel(r'$v_c$ [km/s]')
plt.legend()

import matplotlib.pylab as plt
from os import environ
# load the snapshot and set to physical units
s = pynbody.load('testdata/g15784.lr.01024.gz')

# load the halos
h = s.halos()

# center on the largest halo and align the disk
pynbody.analysis.angmom.faceon(h[1])

# create a profile object for the stars
s.physical_units()
p = pynbody.analysis.profile.Profile(h[1],min=.01,max=250,type='log',ndim=3)
pg = pynbody.analysis.profile.Profile(h[1].g,min=.01,max=250,type='log',ndim=3)
ps = pynbody.analysis.profile.Profile(h[1].s,min=.01,max=250,type='log',ndim=3)
pd = pynbody.analysis.profile.Profile(h[1].d,min=.01,max=250,type='log',ndim=3)

# make the plot
plt.plot(p['rbins'],p['v_circ'],label='total')
plt.plot(pg['rbins'],pg['v_circ'],label='gas')
plt.plot(ps['rbins'],ps['v_circ'],label='stars')
plt.plot(pd['rbins'],pd['v_circ'],label='dm')

plt.xlabel('R [kpc]')
plt.ylabel(r'$v_c$ [km/s]')
plt.legend()

h['te']-=te_max


    logger.info("Making disk rotation curve...")
    
    # Now make a rotation curve for the disk. We'll take everything
    # inside a vertical height of eps*3

    d = h[filt.Disc('1 Mpc', h['eps'].min()*3)]
    
    try :
        
        # attempt to load rotation curve off disk
        r, v_c = np.loadtxt(h.ancestor.filename+".rot."+str(h.properties["halo_id"]), skiprows=1, unpack=True)
        
        pro_d = profile.Profile(d, nbins=100, type='log')
        r_me = pro_d["rbins"].in_units("kpc")
        r_x = np.concatenate(([0], r, [r.max()*2]))
        v_c = np.concatenate(([v_c[0]], v_c, [v_c[-1]]))
        v_c = interp.interp1d(r_x, v_c, bounds_error=False)(r_me)

        logger.info(" - found existing rotation curve on disk, using that")
            
        v_c = v_c.view(array.SimArray)
        v_c.units = "km s^-1"
        v_c.sim = d

        v_c.convert_units(h['vel'].units)
        
        pro_d._profiles['v_circ'] = v_c
        pro_d.v_circ_loaded = True

@Profile.profile_property
def pattern_frequency(pro):
    """Estimate the pattern speed from the m=2 Fourier mode"""
    return pro['fourier']['dphi_dt'][2,:]/2

formkey = 'tform'

    try:
        diskstars['tform']
    except KeyError:
        formkey = 'timeform'

    youngstars = np.where(diskstars[formkey].in_units("Myr") >
                          sim.properties['time'].in_units(
                              "Myr", **sim.conversion_context())
                          - pretime.in_units('Myr'))[0]

    # calculate surface densities
    if radial:
        ps = profile.Profile(diskstars[youngstars], nbins=bins)
        pg = profile.Profile(diskgas, nbins=bins)
    else:
        # make bins 2 kpc
        nbins = rmax * 2 / binsize
        pg, x, y = np.histogram2d(diskgas['x'], diskgas['y'], bins=nbins,
                                  weights=diskgas['mass'],
                                  range=[(-rmax, rmax), (-rmax, rmax)])
        ps, x, y = np.histogram2d(diskstars[youngstars]['x'],
                                  diskstars[youngstars]['y'],
                                  weights=diskstars['mass'],
                                  bins=nbins, range=[(-rmax, rmax), (-rmax, rmax)])

    if clear:
        plt.clf()

    plt.loglog(pg['density'].in_units('Msol pc^-2'),