Set kernel = biff



In [ ]:

from __future__ import division, print_function

# Third-party
import astropy.units as u
from astropy.constants import G as _G
import h5py
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline
pl.style.use("apw-notebook") # comment this out to run at home
importgala.potential as gp
fromgala.units import galactic

import biff
G = _G.decompose(galactic).value




In [ ]:

def _plummer_density(x, y, z, M, r_s):
"""
.. warning::

THIS IS AN INTERNAL FUNCTION -- USE plummer_density() INSTEAD.

"""
r2 = x*x + y*y + z*z
return (3*M / (4*np.pi*r_s**3)) * (1 + r2/r_s**2)**(-5/2.)




In [ ]:

true_M = 1/G
true_r_s = 1.

x = np.logspace(-2,1,512)
xyz = np.zeros((len(x),3))
xyz[:,0] = x

pot = gp.PlummerPotential(m=true_M, b=true_r_s, units=galactic)
true_pot = pot.value(xyz.T).value
true_dens = pot.density(xyz.T).value




In [ ]:

nmax = 16
lmax = 0

Snlm = np.zeros((nmax+1,lmax+1,lmax+1))
Serr = np.zeros((nmax+1,lmax+1,lmax+1))
Tnlm = np.zeros((nmax+1,lmax+1,lmax+1))
Terr = np.zeros((nmax+1,lmax+1,lmax+1))

nlms = []
for n in range(nmax+1):
for l in range(lmax+1):
for m in range(l+1):
nlms.append([n,l,m])

for nlm in nlms:
n,l,m = nlm
print(n,l,m)
(S,S_err),(T,T_err) = biff.compute_coeffs(_plummer_density, nlm=nlm,
M=true_M, r_s=true_r_s, args=(true_M,true_r_s),
epsrel=1E-9)
Snlm[n,l,m] = S
Serr[n,l,m] = S_err
Tnlm[n,l,m] = T
Terr[n,l,m] = T_err

#     nmax = f['nlm'].attrs['nmax']
#     lmax = f['nlm'].attrs['lmax']
#     nlm = np.array(f['nlm'])
#     _Anlm = np.array(f['Anlm'])
#     Anlm_err = np.array(f['Anlm_err'])

#     ix = np.abs(_Anlm) > np.abs(Anlm_err)
#     nlm = nlm[ix]
#     _Anlm = _Anlm[ix]

# Anlm = np.zeros((nmax+1, lmax+1, lmax+1))
# for (n,l,m),A in zip(nlm, _Anlm):
#     Anlm[n,l,m] = A




In [ ]:

pl.semilogy(np.array(nlms)[::2,0], np.abs(Snlm.flat/Snlm[0,0,0])[::2])
pl.xlim(0,10)
pl.ylim(1E-6, 1)




In [ ]:

bfe_dens = biff.density(xyz, Snlm, Tnlm, nmax, lmax, true_M, true_r_s)
bfe_pot = biff.potential(xyz, Snlm, Tnlm, nmax, lmax, G, true_M, true_r_s)




In [ ]:




In [ ]:

fig,axes = pl.subplots(3, 1, figsize=(8,10), sharex=True)

axes[0].set_ylabel(r'$\rho$')
axes[0].loglog(x, true_dens, marker=None)
axes[0].loglog(x, bfe_dens, marker=None)

axes[1].set_ylabel(r'$\Phi$')
axes[1].semilogx(x, true_pot, marker=None)
axes[1].semilogx(x, bfe_pot, marker=None)

axes[2].set_ylabel(r'${\rm d}\Phi/{\rm d}x$')

fig.tight_layout()




In [ ]:

fig,axes = pl.subplots(3, 1, figsize=(8,10), sharey=True)

axes[0].set_ylabel(r'$\rho$')
axes[0].loglog(x, np.abs((true_dens-bfe_dens)/true_dens), marker=None)

axes[1].set_ylabel(r'$\Phi$')
axes[1].loglog(x, np.abs((true_pot-bfe_pot)/true_pot), marker=None)

axes[2].set_ylabel(r'${\rm d}\Phi/{\rm d}x$')

axes[0].set_ylim(1E-12, 1E0)
fig.tight_layout()




In [ ]: