Python cosmology.Planck13類代碼示例

def mk_mock_coords(radeczfile, outfile, simul_cosmo):

    if simul_cosmo == "Planck":
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5

    rad = np.arange(1.0, 67.0, 5.0)

    radecz = h5_arr(radeczfile, "radecz")

    cart = np.zeros(radecz.shape)

    for i, rdz in enumerate(radecz):

        ra = Angle(rdz[0], u.deg)
        dec = Angle(rdz[1], u.deg)

        losd = cosmo.comoving_distance(rdz[2])
        dis = Distance(losd, u.Mpc)

        coord = ICRSCoordinates(ra, dec, distance=dis)

        cart[i, :] = np.array([coord.x, coord.y, coord.z])

    arr2h5(cart, outfile, "coords", mode='w')

def mk_coords(radecfile, outfile, cosmology):

    # Set the cosmology with h free
    if cosmology == "Planck":
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif cosmology == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5

    f_in = h5.File(radecfile)
    radecz = f_in["radecz"]

    f_out = h5.File(outfile)
    cart = f_out.create_dataset("cart_pts", shape=(radecz.shape[0], 3),
                                dtype='float64')

    for i in range(radecz.shape[0]):
        ra = Angle(radecz[i, 0], u.deg)
        dec = Angle(radecz[i, 1], u.deg)

        losd = cosmo.comoving_distance(radecz[i, 2])
        dis = Distance(losd)

        coord = ICRSCoordinates(ra, dec, distance=dis)

        cart[i, :] = np.array([coord.x, coord.y, coord.z])

    f_in.close()
    f_out.close()

def get_comv(cosmology):

    if cosmology == "Planck":
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif cosmology == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5

    return cosmo.comoving_distance

def r200(z,M_200): ## unit: arcsec

	critical_density = Planck13.critical_density(z).si.value # unit: kg/m^3
	m2arcs = np.degrees(1.0/Planck13.angular_diameter_distance(z).si.value)*3600.
	critical_density  = critical_density*Planck13.H0/100./m2arcs**3/2e30 # unit: Msun/arcsec^3/h

	r200_cube = 200*critical_density/M_200/(4./3.*np.pi)
	r200 = r200_cube**(1./3.)

	return r200

def get_inv_efunc(cosmology):

    if cosmology == "Planck":
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif cosmology == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5

    return cosmo.inv_efunc

def mk_mock_srch(radecfile, nzdictfile, Nsph, simul_cosmo):

    if simul_cosmo == "Planck":
        # First make h free
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5
    comv = cosmo.comoving_distance

    radecarr = h5_arr(radecfile, "good_pts")
    nzdict = json.load(open(nzdictfile))

    Nrands = radecarr.shape[0]
    Narrs = Nsph / Nrands
    remain = Nsph % Nrands

    radecz = np.zeros((Nsph, 3))

    for i in range(Narrs):

        start = Nrands * i
        stop = Nrands * (i + 1)
        radecz[start:stop, :2] = radecarr[:, :]

    endchunk = Nrands * (Narrs)
    radecz[endchunk:, :2] = radecarr[:remain, :]

    rad = np.arange(1.0, 67.0, 5.0)
    zlo = nzdict["zlo"]
    zhi = nzdict["zhi"]

    radeczlist = len(rad) * [radecz]

    for r_i, r in enumerate(rad):

        dis_near = Distance(comv(zlo) + r, u.Mpc)
        dis_far = Distance(comv(zhi) - r, u.Mpc)

        z_a = dis_near.compute_z(cosmology=cosmo)

        z_b = dis_far.compute_z(cosmology=cosmo)

        randz = (z_a ** 3 + \
                 (z_b ** 3 - z_a ** 3) * np.random.rand(Nsph)) ** (1. / 3.)

        radeczlist[r_i][:, 2] = randz[:]

        arr2h5(radeczlist[r_i], "{0}/{1}/mocks/mock_srch_pts.hdf5".format(os.path.dirname(radecfile), simul_cosmo), "radecz_{0}".format(str(r_i * 5 + 1)))

def mock_vpf(mock_cart_coords, spheresfile, simul_cosmo, rad):

    gals = h5_arr(mock_cart_coords, "coords")

    print gals

    name = mock_cart_coords.split("/")[-1].split(".")[0]

    if simul_cosmo == "Planck":
        # First make h free
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5
    comv = cosmo.comoving_distance

    gal_baum = cKDTree(gals)

    spheres = h5_arr(spheresfile, "radecz_{0}".format(str(int(rad))))

    print spheres

    for i, sphere in enumerate(spheres):

        rang = Angle(sphere[0], u.deg)
        decang = Angle(sphere[1], u.deg)

        dis = Distance(comv(sphere[2]), u.Mpc)

        coord = ICRSCoordinates(rang, decang, distance=dis)

        sph_cen = np.array([coord.x, coord.y, coord.z])

        nn = gal_baum.query(sph_cen)

        print "rad: ", rad, ", sphere: ", i

        f = open("{0}/vpf_out/ascii/{1}_{2}.dat".format(os.path.dirname(spheresfile), name, str(int(rad))), 'a')

        if not nn[0] < rad:
            f.write("1\n")
        else:
            f.write("0\n")

        f.close()

    def __init__(self, redshifts, cosmology='Planck13', cm='DuttonMaccio'):
        if type(redshifts) != np.ndarray:
            redshifts = np.array(redshifts)
        if redshifts.ndim != 1:
            raise ValueError("Input redshift array must have 1 dimension.")
        if np.sum(redshifts < 0.) > 0:
            raise ValueError("Redshifts cannot be negative.")

        if cosmology == 'Planck13':
            from astropy.cosmology import Planck13 as cosmo
        elif cosmology == 'WMAP9':
            from astropy.cosmology import WMAP9 as cosmo
        elif cosmology == 'WMAP7':
            from astropy.cosmology import WMAP7 as cosmo
        elif cosmology == 'WMAP5':
            from astropy.cosmology import WMAP5 as cosmo
        else:
            raise ValueError('Input cosmology must be one of: \
                              Planck13, WMAP9, WMAP7, WMAP5.')
        self._cosmo = cosmo

        if cm == 'DuttonMaccio':
            self._cm = 'DuttonMaccio'
        elif cm == 'Prada':
            self._cm = 'Prada'
        elif cm == 'Duffy':
            self._cm = 'Duffy'
        else:
            raise ValueError('Input concentration-mass relation must be \
                              one of: DuttonMaccio, Prada, Duffy.')

        self.describe = "Ensemble of galaxy clusters and their properties."
        self.number = redshifts.shape[0]
        self._z = redshifts
        self._rho_crit = cosmo.critical_density(self._z)
        self._massrich_norm = 2.7 * (10**13) * units.Msun
        self._massrich_slope = 1.4
        self._df = pd.DataFrame(self._z, columns=['z'])
        self._Dang_l = cosmo.angular_diameter_distance(self._z)
        self._m200 = None
        self._n200 = None
        self._r200 = None
        self._rs = None
        self._c200 = None
        self._deltac = None

def mock_vpf(mock_cart_coords, spheresfile, simul_cosmo):

    gals = h5_arr(mock_cart_coords, "coords")
    name = mock_cart_coords.split("/")[-1].split(".")[0]

    if simul_cosmo == "Planck":
        # First make h free
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5
    comv = cosmo.comoving_distance

    gal_baum = cKDTree(gals)

    rad = np.arange(1.0, 67.0, 5.0)

    for r_i, r in enumerate(rad):

        spheres = h5_arr(spheresfile, "radecz_{0}".format(str(r_i * 5 + 1)))
        voids = np.zeros(spheres.shape[0])

        for i, sphere in enumerate(spheres):

            rang = Angle(sphere[0], u.deg)
            decang = Angle(sphere[1], u.deg)

            dis = Distance(comv(sphere[2]), u.Mpc)

            coord = ICRSCoordinates(rang, decang, distance=dis)

            sph_cen = np.array([coord.x.value, coord.y.value, coord.z.value])

            nn = gal_baum.query(sph_cen)

            print "rad: ", r, ", sphere: ", i

            if not nn[0] < r:

                voids[i] = 1

        arr2h5(voids,
                "{0}/vpf_out/{1}.hdf5".format(os.path.dirname(spheresfile), name),
                "voids_{0}".format(str(r_i * 5 + 1)))

def bolometric(model, luminosity=True):
    phase = model.source._phase
    flux = np.sum(model.source.flux(phase, model.source._wave) * dwave, axis=1)
    if luminosity:
        z = model.get('z')
        dl = Planck13.luminosity_distance(z).to('cm').value
        L = 4 * np.pi * flux * dl * dl
        return L
    return flux

def z_to_time(z,cosmology=None):
    """
    gives the time, in years, since the big bang given an astropy cosmology
    and a redshift.  if no cosmology passed in, assumes Planck13
    """
    if cosmology is None:
        from astropy.cosmology import Planck13 as cosmology
    from yt import YTArray
    t = cosmology.age(z)
    return YTArray(t*t.unit.in_units('yr'),'yr')

def main():
    fil = pysysp.StarSpectrum()
    fil.setwavelength(np.arange(1000, 30000, 0.4))
    fil.setflux( 1e-15 * (fil.wavelength/10000.0)**(-1.5))
    i = pysysp.BandPass('/Users/jselsing/github/pysysp/pysysp/filters/ugriz/i.dat')


    dl = (cosmo.luminosity_distance(1.2)) * 1e5
    Mz0 = -5 * np.log10(dl.value) + fil.apmag(band=i, mag='AB')
    print('Mz0 = '+str(Mz0))


    i.wavelength /= (1 + 2)
    # fil.flux = fil.flux * (1 + 2)   

    i.update_bandpass()

    dl = (cosmo.luminosity_distance(1.2)) * 1e5
    Mz2 = -5 * np.log10(dl.value) + fil.apmag(band=i, mag='AB') #- 2.5 * np.log10(1 + 2)
    print('Mz2 = '+str(Mz2))




    # i.wavelength *= (1 + 2)

    # i.update_bandpass()

    # fil.wavelength = fil.wavelength * (1 + 2)
    # fil.flux = fil.flux / (1 + 2)
    # dl = (cosmo.luminosity_distance(1.2)) * 1e5
    # Mz2 = -5 * np.log10(dl.value) + fil.apmag(band=i, mag='AB')
    # print('Mz2 = '+str(Mz2))
    print('Mz0 - Mz2 = '+str(Mz0 - Mz2))

def slice_z_cosmo(slice, z, radius):
    """Slices in redshift considering a distance to cut"""
    zmax = z_at_value(cosmo.comoving_distance, cosmo.comoving_distance(z) + radius)
    zmin = z_at_value(cosmo.comoving_distance, cosmo.comoving_distance(z) - radius)

    print zmin - zmax

    slice = slice[np.where(slice["zb_1"] < zmax)]
    slice = slice[np.where(slice["zb_1"] > zmin)]
    return slice

def producing_lensed_images(xi1,xi2,source_cat,lens_cat):

    xlc1 = lens_cat[0]
    xlc2 = lens_cat[1]
    ql   = lens_cat[2]
    rlc  = 0.0
    rle  = re_sv(lens_cat[3],lens_cat[5],source_cat[7])       #Einstein radius of lens Arc Seconds.

    phl = lens_cat[4]
    #----------------------------------------------------------------------
    ai1,ai2,mua = lens_equation_sie(xi1,xi2,xlc1,xlc2,ql,rlc,rle,phl)

    yi1 = xi1-ai1*0.0
    yi2 = xi2-ai2*0.0

    ysc1 = source_cat[0]
    ysc2 = source_cat[1]
    mag_tot = source_cat[2]
    Reff_arc = np.sqrt(source_cat[3]*source_cat[4])
    qs   = np.sqrt(source_cat[3]/source_cat[4])
    phs  = source_cat[5]
    ndex = source_cat[6]
    zs = source_cat[7]

    #g_limage = sersic_SB_2D_tot(yi1,yi2,ysc1,ysc2,mag_tot,Reff_arc,qs,phs,zs)
    Da_s = p13.angular_diameter_distance(2.0).value
    Reff = Reff_arc*Da_s/apr*1e6 # pc
    Ieff = sersic_mag_tot_to_Ieff(mag_tot,Reff,ndex,zs)
    g_limage = sersic_2d(yi1,yi2,ysc1,ysc2,Ieff,Reff_arc,qs,phs,ndex)

    return g_limage

def find_rvirial(dd, ds, center, start_rad = 0, delta_rad_coarse = 20, delta_rad_fine = 1, rad_units = 'kpc'):
    vir_check = 0
    r0 = ds.arr(start_rad, rad_units)
    critical_density = cosmo.critical_density(ds.current_redshift).value   #is in g/cm^3
    max_ndens_arr=center

    while True:
        r0_prev = r0
        r0 = r0_prev + ds.arr(delta_rad_coarse, rad_units)
        v_sphere = ds.sphere(max_ndens_arr, r0)
        dark_mass 	= v_sphere[('DarkMatter', 'Mass')].in_units('Msun').sum()	
        rho_internal = dark_mass.in_units('g')/((r0.in_units('cm'))**3.*(pi*4/3.))
        print rho_internal, r0, dark_mass
        if rho_internal < 200*ds.arr(critical_density,'g')/ds.arr(1.,'cm')**3.:
            #now run fine test
            print('\tNow running fine search on the virial radius')
            r0 = r0_prev
            while True:
                r0 += ds.arr(delta_rad_fine, rad_units)
                v_sphere = ds.sphere(max_ndens_arr, r0)
                dark_mass 	= v_sphere[('DarkMatter', 'Mass')].in_units('Msun').sum()	
                rho_internal = dark_mass.in_units('g')/((r0.in_units('cm'))**3.*(pi*4/3.))
                print rho_internal, r0
                if rho_internal < 200*ds.arr(critical_density,'g')/ds.arr(1.,'cm')**3.:
                    rvir = r0
                    print dark_mass.in_units('Msun'), rvir.in_units('kpc')
                    return rvir