Python Planck13.__init__方法代碼示例

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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')

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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()

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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)))

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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()

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
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)))

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
import numpy as np
from astropy.cosmology import Planck13, WMAP5
from scipy.interpolate import interp2d
from scipy.spatial import cKDTree


simul_cosmo = "WMAP"


def spherical_cap(h):
    return 0.75 * (h ** 2) * (1 - h / 3)


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

Nsph = 10000000
rad = np.arange(5.0, 66.0, 5.0)

As = np.arange(0.0, 1.0, 0.05)
Bs = np.arange(0.0, 1.0, 0.05)

splarr = np.loadtxt("test_dat/edge_splarr.dat")
A, B = np.meshgrid(As, Bs)

inty = interp2d(A[0, :], B[:, 0], splarr)

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
def vpf(dat_dir, Nsph, simul_cosmo, rad):

    # Grab the data coordinates
    gals = h5_arr("./dat/out/{0}/{1}/gals_cart_coords.hdf5".
                      format(dat_dir, simul_cosmo), "cart_pts")

    # Get details about the redshift interval being considered
    nbar_dict = json.load(open("./dat/out/{0}/{1}/nbar_zrange.json".
                                   format(dat_dir, simul_cosmo)))
    zlo = nbar_dict["zlo"]
    zhi = nbar_dict["zhi"]

    # Get the search points
    good_pts = h5_arr("./dat/out/{0}/srch_radec.hdf5".format(dat_dir), "good_pts")
    bad_pts = h5_arr("./dat/out/{0}/veto.hdf5".format(dat_dir),
                     "bad_pts")

    # Set angular radius of effective area around bad points
    bad_r = np.arccos(1.0 - (np.pi * 9.8544099e-05) / (2 * 180 ** 2))
    bad_r_deg = np.rad2deg(bad_r)

    # Set the cosmology with h free
    # Here the cosmology is based on WMAP (for first MultiDark simulation)
    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

    # Build the trees

    # galaxy tree
    gal_baum = cKDTree(gals)

    # tree of bad points (angular coordinates on unit sphere)
    bad_xyz = radec2xyz(bad_pts)
    veto_baum = cKDTree(bad_xyz)

    # Initialise final output arrays
#    rad = np.arange(1.0, 67.0, 5.0)  doing it one radius at a time
#    P_0 = np.zeros(rad.shape)

    # No. of spheres and norm
#     Nsph_arr = Nsph * np.array(4 * [0.01] + 4 * [0.1] + 4 * [1.0])
#     norm = 1. / Nsph_arr
#    norm = 1. / Nsph

    rand_i = 0

    for r_i, r in enumerate(rad):

        # start the count of successful voids
        count = 0

        # Custom zrange for sphere size
        dis_near = Distance(comv(zlo).value + r, u.Mpc)
        dis_far = Distance(comv(zhi).value - r, u.Mpc)

        z_a = dis_near.compute_z(cosmology=cosmo)

        z_b = dis_far.compute_z(cosmology=cosmo)

        for i in range(Nsph):  # _arr[r_i]):

            # compensate for finite length of mask file
            rand_i = rand_i % 999999

            radec = good_pts[rand_i, :]

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

            randz = (z_a ** 3 + \
                     (z_b ** 3 - z_a ** 3) * np.random.rand(1)[0]) ** (1. / 3.)
            dis = Distance(comv(randz), 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:

                # add instance to probability count
                count += 1

                # record quality of sphere using spline values for intersection
                # with bad points

                # Get radius of circular projection of sphere
                R = np.arcsin(r / np.sqrt(np.sum(sph_cen[:] ** 2)))

                # Get coordinates of circle centre on unit sphere
                crc_cen = radec2xyz(radec)[0]
#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
def process_nbar(nbarfile, nz_dict_file, cosmology, radeczfile=None):
    """
    Parameters
    ---------

    nbarfile : str
        the path to and name of the corrected nbar file
    nz_dict_file : str
        path to and name of the json file with the nbar dict
    cosmology : str, "WMAP" or "Planck"
        the cosmology to compute shell volumes with
    radeczfile : str, "data" or "mock"
        the data or mock file to process
    """

    # magic number for width around maximum
    Q = 0.65
    # magic number for shell vol computation
    Nfrac = (6769.0358 * np.pi) / 129600

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

    nbar_corr = np.loadtxt(nbarfile)
    nz_dict = {"tophat height for zrange": Q}

    # Cut out the first bit of crap (works for CMASS, dunno about LOWZ)
    ind03 = np.abs(nbar_corr[:, 0] - 0.3).argmin()

    nbar_corr = nbar_corr[ind03:, :]

    zcen = nbar_corr[:, 0]
    z_near = nbar_corr[:, 1]
    z_far = nbar_corr[:, 2]
    corr_gal_counts = nbar_corr[:, 6]

    nbar = []
    shell_vols = []

    for i in range(len(zcen)):

        shell_vols.append(Nfrac * calc_shell_vol(comv, z_near[i], z_far[i], zcen[i]))
        nbar.append(corr_gal_counts[i] / shell_vols[i])

    nbar = np.array(nbar)

    # Find nbar peak and index
    max_nbar = np.max(nbar)
    max_i = int(np.where(nbar == max_nbar)[0])

    nz_dict["max_nbar_corr"] = max_nbar
    nz_dict["nbar_corr_tophat"] = Q * max_nbar
    nz_dict["z_nbar_max"] = zcen[max_i]

    # get the interval edge indices
    L = np.abs(nbar[:max_i] - max_nbar * Q).argmin()
    R = max_i + np.abs(nbar[max_i:] - max_nbar * Q).argmin()

    nbar = nbar[L:R + 1]
    shell_vols = shell_vols[L:R + 1]

    nz_dict["zlo"] = zcen[L]
    nz_dict["zhi"] = zcen[R]

    nz_dict["avg_nbar_corr"] = np.average(nbar)
    nz_dict["total_shell_vol"] = np.sum(shell_vols)

    if radeczfile:

        radecz = h5_arr(radeczfile, "radecz")

        # Make the redshift cut in the nbar array with right cosmology
        nbar_corr = nbar_corr[(nz_dict["zlo"] <= nbar_corr[:, 0]) * \
                            (nbar_corr[:, 0] <= nz_dict["zhi"])]

        # Get binning those observed galaxies
        zbinedges = np.append(nbar_corr[0, 1], nbar_corr[:, 2])

        # Find the counts per bin
        H = np.histogram(radecz[:, 2], bins=zbinedges)

        # The number to downsample to in each bin
        # (multiply bin number by the relative fraction determined from
        #  corrected distribution of nbar)
        num_down = np.rint((nz_dict["nbar_corr_tophat"] / nbar[:]) * H[0])
        num_down = num_down.astype(int)

        # make a mask for the final array for analysis within the redshift limits
        finmask = np.array(radecz.shape[0] * [False])

        for i, nd in enumerate(num_down):
            """Turn on the right amount of galaxies in each bin."""
            zbin_ids = np.where(((zbinedges[i] < radecz[:, 2]) * (radecz[:, 2] <= zbinedges[i + 1])) == True)

            if zbin_ids[0].shape[0] == 0:
#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import __init__ [as 別名]
def box_completeness(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

    rad = np.arange(5.0, 66.0, 5.0)

    As = np.arange(0.0, 1.0, 0.05)
    Bs = np.arange(0.0, 1.0, 0.05)

    splarr = np.loadtxt("test_dat/edge_splarr.dat")
    A, B = np.meshgrid(As, Bs)

    inty = interp2d(A[0, :], B[:, 0], splarr)

    # this number from survey I think
    # should be 2% of "sky" area
    # had 138621 for some reason, now used 2% to get
    bad_pts = 1000 * np.random.rand(34834737093, 2)
    bad_r = 0.0004275  # determined from quick calculation for now
    bad_A = np.pi * bad_r ** 2

    badbaum = cKDTree(bad_pts)

    for r_i, r in enumerate(rad):

        spheres = 1000 * np.random.rand(Nsph, 2)

        bound_bool = (spheres[:, 0] < r) + (spheres[:, 1] < r) + \
                     ((1000 - spheres[:, 0]) < r) + \
                     ((1000 - spheres[:, 1]) < r)
        bad_inds = np.where(bound_bool == True)
        badsphs = spheres[bound_bool]

        pickle_bool = ((badsphs[:, 0] ** 2 + badsphs[:, 1] ** 2) < r) + \
               (((1000 - badsphs[:, 0]) ** 2 + (1000 - badsphs[:, 0]) ** 2)
                   < r)
        pickle_inds = bad_inds[pickle_bool]

        for i, sph in enumerate(spheres):

            badvol = 0.

            pierce_pts = badbaum.query_ball_point(sph, r)

            for pt in pierce_pts:
                # retrieve coordinates of points within sphere
                pt_coord = bad_pts[pt]
                # calculate fractional projected distance from centre
                dis = np.sqrt((sph[0] - pt_coord[0]) ** 2 + \
                              (sph[1] - pt_coord[1]) ** 2) / r
                # calculate length pierced through sphere
                l = 2 * np.sqrt(1 - dis ** 2)

                badvol += l * bad_A

            # check if sphere at boundary
            if i in bad_inds:

                if i in pickle_inds:

                    badvol += inty(sph[0] / r, sph[1] / r)

                else:
                    if sph[0] < r:
                        badvol += spherical_cap(1 - sph[0] / r)
                    elif 1000 - sph[0] < r:
                        badvol += spherical_cap(1 - (1000 - sph[0]) / r)

                    if sph[1] < r:
                        badvol += spherical_cap(1 - sph[1] / r)
                    elif 1000 - sph[1] < r:
                        badvol += spherical_cap(1 - (1000 - sph[1]) / r)

            f = open("test_dat/simul_badvol.dat", 'a')
            f.write("{0}\n".format(badvol))
            f.close()