In [1]:

    

from __future__ import division, print_function

# Third-party
from astropy import log as logger
import emcee
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
from scipy.misc import logsumexp
import triangle
%matplotlib inline

# Custom
import streams.coordinates as stc
import streamteam.io as io
import streamteam.dynamics as sd
import streamteam.potential as sp
from streamteam.units import galactic


    

In [2]:

    

potential = sp.LeeSutoNFWPotential(v_h=0.5, r_h=20., 
                                   a=1., b=1., c=1., units=galactic)


    

In [3]:

    

def true_m_enc(potential, r):
    rh = potential.parameters['r_h']
    vh = potential.parameters['v_h']
    G = potential.c_instance.G
    rr = r/rh
    
    return rh*vh**2/G * (np.log(1+rr) - rr/(1+rr))


    

In [4]:

    

def est_m_enc(potential, xyz):
    G = potential.c_instance.G    
    
    r = np.linalg.norm(xyz, axis=-1)
    h = 0.01
    epsilon = h * xyz/r[:,None]
    dPhi_dr = (potential.value(xyz + epsilon) - potential.value(xyz - epsilon)) / (2*h)
    print(h)
    
    return np.abs(1/G * r**2 * dPhi_dr)

def est_m_enc_cy(potential, q):
    G = potential.c_instance.G
    epsilon = np.empty((1,3))
    mass = np.zeros(len(q))
    for i in range(len(q)):

        # Step-size for estimating radial gradient of the potential
        r = np.sqrt(q[i,0]*q[i,0] + q[i,1]*q[i,1] + q[i,2]*q[i,2])

        h = 0.01
        
        for j in range(3):
            epsilon[0,j] = h * q[i,j]/r + q[i,j]
        dPhi_dr = potential.value(epsilon)

        for j in range(3):
            epsilon[0,j] = h * q[i,j]/r - q[i,j]
        dPhi_dr = dPhi_dr - potential.value(epsilon)
        
        dPhi_dr /= (2*h)
        
        mass[i] = np.abs(r*r * dPhi_dr / G)
    
    return mass


    

In [5]:

    

rs = np.linspace(10.,200.,25)
xs = np.zeros((len(rs),3))
xs[:,0] = rs

mass = true_m_enc(potential, rs)
mass_appx = est_m_enc(potential, xs)
mass_cy = est_m_enc_cy(potential, xs)
mass_cls = potential.c_instance.mass_enclosed(xs)
plt.plot(rs, mass, marker=None)
plt.plot(rs, mass_appx, linestyle='none', marker='o')
plt.plot(rs, mass_cy, linestyle='none', marker='o')
plt.plot(rs, mass_cls, linestyle='none', marker='o')


    

    




0.01






    
Out[5]:





[<matplotlib.lines.Line2D at 0x10af855d0>]






    

In [5]: