In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [97]:

    

def merger_time_scale(m_ratio=1.0, hubble=1.0, eta=0.5, r_c=0.1, r_vir=1.0):
    """
    Merger Time-scale for a halo merger. Reference: Boylan-Kolchin (2007)
    
    Inputs: 
        m_ratio: mass ratio (<1.0) of the merger. The original paper only treated m_ratio<0.1
        hubble: hubble constant in units of Gyr^-1
        eta: circularity parameter j/j_c(E), where j_c(E) is the 
            specific angular momentum of a circular orbit of energy E.
            eta is related to excentricity as eta = (1-e^2)^1/2
        r_c: circular radius of an orbit with energy E (in Mpc)
        r_vir: virual radius (in Mpc)
    """
    
    t_df = 0.0216/hubble
    t_df = t_df * m_ratio**(1.3)/(log(1.0+m_ratio))
    t_df = t_df * exp(1.9*eta) * (r_c/r_vir)
    return t_df


    

In [98]:

    

G_GRAV = 4.38E-15 #Longitud en Mpc, masa en Msun y tiempo en Gyr.
OVERDENSITY = 360
KM_PER_SEC_TO_MPC_PER_GYR = 0.00105 # converts km/s into Mpc/Gyr 
OMEGA_BARYONS = 0.05
OMEGA_DM = 0.27
HUBBLE = 70.0
HUBBLE_T = HUBBLE * KM_PER_SEC_TO_MPC_PER_GYR
RHO_U = 7.5E10 # average density of the universe at z=0.0 in Msun/Mpc^3


def func_rhalo(mhalo=1.0):
    """
    Input: 
        mhalo: Virial mass in units of Msun
    Output: 
        The virial radius in units of kpc
    """
    part_a= mhalo/(OVERDENSITY * RHO_U * (4.0/3.0)*pi)
    part_b = (part_a**(1.0/3.0))

    return part_b*1000.0

def func_sigma_halo(mhalo=1.0):
    """
        Inputs:
            mhalo: halo mass in Msun
        Outputs:
            sigma_halo: velocity dispersion inside the halo in km/s
    """
    rhalo = func_rhalo(mhalo=mhalo)/1000.0
    sigma_SIS = sqrt(G_GRAV*mhalo/rhalo)
    sigma_SIS = sigma_SIS/KM_PER_SEC_TO_MPC_PER_GYR
    return sigma_SIS


    

In [99]:

    

func_sigma_halo(mhalo=3E12)


    

    
Out[99]:





199.90802873488232

In [100]:

    

merger_time_scale(m_ratio=0.1, hubble=HUBBLE_T, eta=0.0, r_c=0.1, r_vir=1.0)


    

    
Out[100]:





0.015453509491164806

In [111]:

    

halo_mass  = array([1.0E12, 2.0E12, 3.0E12])
sat_mass = 1.0E11
n_points = size(halo_mass)
time_r_A = zeros(n_points)
time_r_B = zeros(n_points)
radius_r = zeros(n_points)


    

In [120]:

    

for i in range(n_points):
    radius_r[i] = func_rhalo(mhalo=halo_mass[i])
    time_r_A[i] = merger_time_scale(m_ratio=sat_mass/halo_mass[i], hubble=HUBBLE_T,\
                                    eta=1.0, r_c=1.0 * radius_r[i], r_vir=radius_r[i])
    time_r_B[i] = merger_time_scale(m_ratio=sat_mass/halo_mass[i], hubble=HUBBLE_T,\
                                    eta=0.6, r_c=1.0 * radius_r[i], r_vir=radius_r[i])


    

In [121]:

    

plot(halo_mass, time_r_A)
plot(halo_mass, time_r_B)


    

    
Out[121]:





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






    

In [ ]: