How to use the pynbody.analysis.angmom function in pynbody

fifmyrf = filt.LowPass('age','15 Myr')
twokpcf = filt.Sphere('2 kpc')
i=1
if (len(sys.argv) > 2):
    photiords = np.genfromtxt(sys.argv[2],dtype='i8')
    frac = np.float(len(np.where(np.in1d(photiords,h[i]['iord']))[0]))/len(photiords)
    print 'i: %d frac: %.2f'%(i,frac)
    while(((frac) < 0.5) & (i<100)): 
        i=i+1
        frac = np.float(len(np.where(np.in1d(photiords,h[i]['iord']))[0]))/len(photiords)
        print 'i: %d frac: %.2f'%(i,frac)
else:
    while len(h[i].star) <2: i=i+1

if (i==100): sys.exit()
pynbody.analysis.angmom.faceon(h[i])
try:
    s.g['tconv'] = 4.0 * s.g['smooth'] / np.sqrt(s.g['u']+s.g['uNoncool'])
    s.g['tconv'].units = s.s['tform'].units
except:
    pass
s.physical_units()
#diskf = filt.Disc(np.max(h[i].s[diskf]['r']),np.max(h[i].s[diskf]['z']))
notdiskf = filt.Not(diskf)
notcooldenf = filt.And(filt.LowPass('rho','0.1 m_p cm^-3'),
                       filt.HighPass('temp',3e4))
hot=filt.HighPass('temp',3e4)
cold=filt.LowPass('temp',3e4)
Jtot = np.sqrt(((np.multiply(h[i]['j'].transpose(),h[i]['mass']).sum(axis=1))**2).sum())
W = np.sum(h[i]['phi']*h[i]['mass'])
K = np.sum(h[i]['ke']*h[i]['mass'])
absE = np.fabs(W+K)

def fourier_profile(sim, center=True, disk_height='2 kpc', nbins=50,
                    pretime='2 Gyr', r_units='kpc', bin_spacing='equaln',
                    clear=True, min=False, max=False, filename=None, **kwargs):
    """
    Centre on potential minimum, align so that the disk is in the
    x-y plane, then plot the amplitude of the 2nd fourier mode as a 
    function of radius.

    **needs description of the keyword arguments**

    """

    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(

import pynbody
import matplotlib.pylab as plt

# 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])

# convert all units to something reasonable (kpc, Msol etc)
s.physical_units()

# create a profile object for the stars (by default this is a 2D profile)
p = pynbody.analysis.profile.Profile(h[1].s,min=.01,max=50)

# make the figure and sub plots
f, axs = plt.subplots(1,2,figsize=(14,6))

# make the plot
axs[0].plot(p['rbins'],p['density'], 'k')
axs[0].semilogy()
axs[0].set_xlabel('R [kpc]')
axs[0].set_ylabel(r'$\Sigma_{\star}$ [M$_{\odot}$ kpc$^{-2}$]')

bin_spacing='equaln', clear=True, quick=False,
                   filename=None, min=False, max=False, yrange=False,
                   legend=False, parts=False, axes=False, **kwargs):
    """

    Centre on potential minimum, align so that the disk is in the
    x-y plane, then use the potential in that plane to generate and
    plot a rotation curve.

    **needs documentation/description of the keyword arguments**

    """
    import pylab as p

    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()

    pro = profile.Profile(sim, type=bin_spacing, nbins=nbins,
                          min=min_r, max=max_r)

    r = pro['rbins'].in_units(r_units)
    if quick:
        v = pro['rotation_curve_spherical'].in_units(v_units)

*align* : whether to align the snapshot based on the angular momentum in the central region (default = True)

    *halo* : halo to center on (default = 0)
    """

    s = pynbody.load(output)
    hop_center(s, halo)

    st = s[pynbody.filt.Sphere('100 kpc')]

    cen = pynbody.analysis.halo.center(
        st, retcen=True, mode='ssc', verbose=True)

    if align:
        pynbody.analysis.angmom.faceon(
            st.s, disk_size='10 kpc', cen=cen, mode='ssc')
    else:
        s['pos'] -= cen

    s['pos'].convert_units('kpc')
    s['vel'].convert_units('km s^-1')

    return s

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]')

def half_light_r(sim, band='v',center=True):
    '''Calculate half light radius

    Calculates entire luminosity of simulation, finds half that, sorts
    stars by distance from halo center, and finds out inside which radius
    the half luminosity is reached.
    '''
    import pynbody
    import pynbody.filt as f
    sim = sim.s
    half_l = halo_lum(sim, band=band) * 0.5
    if center: pynbody.analysis.angmom.faceon(sim)
    indr = np.argsort(sim['r'])
    sun_abs_mag = {'u':5.56,'b':5.45,'v':4.8,'r':4.46,'i':4.1,'j':3.66,
                   'h':3.32,'k':3.28}[band]
    lumcumsum = np.cumsum(10.0 ** ((sun_abs_mag - 
                                    sim[indr][band + '_mag']) / 2.5))

    return sim[indr][np.argwhere(lumcumsum < half_l).max()]['r']

import pynbody
import matplotlib.pylab as plt

# load the snapshot and set to physical units
s = pynbody.load('testdata/g15784.lr.01024.gz')
s.physical_units()

# load the halos
h = s.halos()

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

#create an image using the default bands (i, v, u)
pynbody.plot.stars.render(s,width='20 kpc')

import pynbody
import pynbody.plot.sph as sph
import matplotlib.pylab as plt

# load the snapshot and set to physical units
s = pynbody.load('testdata/g15784.lr.01024.gz')
s.physical_units()

# load the halos
h = s.halos()

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

# create the subplots
f, axs = plt.subplots(1,2,figsize=(14,6))

#create a simple slice showing the gas temperature, with velocity vectors overlaid
sph.velocity_image(h[1].g, vector_color="cyan", qty="temp",width=50,cmap="YlOrRd", 
                   denoise=True,approximate_fast=False, subplot=axs[0], show_cbar = False)

#you can also make a stream visualization instead of a quiver plot
pynbody.analysis.angmom.faceon(h[1])
s['pos'].convert_units('Mpc')
sph.velocity_image(s.g, width='3 Mpc', cmap = "Greys_r", mode='stream', units='Msol kpc^-2', 
                   density = 2.0, vector_resolution=100, vmin=1e-1,subplot=axs[1], 
                   show_cbar=False, vector_color='black')

import pynbody
import pynbody.plot.sph as sph
import matplotlib.pylab as plt

# load the snapshot and set to physical units
s = pynbody.load('testdata/g15784.lr.01024.gz')
s.physical_units()

# load the halos
h = s.halos()

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

#create a simple slice showing the gas temperature
sph.image(h[1].g,qty="temp",width=50,cmap="YlOrRd", denoise=True,approximate_fast=False)