%matplotlib inline
from matplotlib.pylab import *
rcParams['figure.figsize'] = (10,6)
rcParams['font.size'] = 18

import pynbody
s = pynbody.load('/Users/rokstar/Downloads/testdata/g15784.lr.01024.gz')

s

len(s)

len(s.stars)

len(s.g), len(s.gas), len(s.dark)

s.properties

s.properties['time'].in_units('Gyr')

s.keys()

s.loadable_keys()

s['pos'].in_units('kpc')

s['mass'], s['vel'], s.g['HeI'], s.s['posform']

s.physical_units()

s['vt'],s['vr']

s.derivable_keys()[0:10]

h = s.halos()

pynbody.analysis.halo.center(h[1],mode='ssc')

h1 = h[1]
h1

len(h1.g), len(h1.s)

h[1]['mass'].sum().in_units('1e12 Msol')

im = pynbody.plot.image(s.d,width='50 Mpc', units='Msol kpc^-2', cmap=cm.Greys)

pynbody.plot.image(h[1].g,width='200 kpc', av_z=True, cmap=cm.Blues)

pynbody.analysis.angmom.faceon(h[1],cen=(0,0,0))

s.rotate_x(90)
pynbody.plot.image(h[1].g,width='50 kpc',cmap=cm.Blues)

pynbody.plot.image(h1.g,qty='temp',width='40 kpc', cmap=cm.YlOrRd)
s.rotate_x(-90)

p = pynbody.analysis.profile.Profile(h1,max='200 kpc',min='.1 kpc',type='log',nbins=100)
ps = pynbody.analysis.profile.Profile(h1.s,max='50 kpc',nbins=50)
pg = pynbody.analysis.profile.Profile(h1.g,max='50 kpc',nbins=50)
plot(p['rbins'], p['density'], label = 'all')
plot(ps['rbins'], ps['density'], label = 'stars')
plot(pg['rbins'], pg['density'], label = 'gas')
semilogy()
xlim(0,50)
legend()
xlabel('$R$ [kpc]')
ylabel('$\Sigma$ [M$_{\odot}$ kpc$^{-2}$]')

plot(p['rbins'], p['v_circ'], label = 'all')
plot(ps['rbins'], ps['v_circ'], label = 'stars')
plot(pg['rbins'], pg['v_circ'], label = 'gas')
xlim(0,50)
xlabel('$R$ [kpc]')
ylabel('$V_{circ}$')
legend()

plot(ps['rbins'],ps['age'])
xlabel('$R$ [kpc]')
ylabel('Age [Gyr]')

plot(ps['rbins'],ps['d_v_circ'])
xlabel('$R$ [kpc]')
ylabel('d$v_{c}$/d$R$')

plot(ps['rbins'],ps['vr_rms'])
xlabel('R [kpc]')
ylabel('$v_{rms}$')

agefilt = pynbody.filt.BandPass('age', '1 Gyr','2 Gyr')
agefilt

sphere = pynbody.filt.Sphere('50 kpc')
h1[sphere].s

h1.s[sphere]['mass'].sum().in_units('1e12 Msol')

h1.s[sphere&agefilt]['mass'].sum().in_units('1e12 Msol')


s['mass'] + s['x']

s['mass']*s['x']