Python FlatLambdaCDM.critical_density方法代码示例

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import critical_density [as 别名]
	def __init__(self):
		ombh2 = 0.022
		omb0 = ombh2/param_dict['h']**2.
		#param_dict['omega_m'] = 0.3089
		from astropy.cosmology import FlatLambdaCDM
		cosmo = FlatLambdaCDM(H0 = param_dict['h']*100., Om0 = param_dict['omega_m'], Ob0 = omb0)

		self.rho_bar = cosmo.critical_density(0.).to('M_sun/Mpc3').value * cosmo.Om0
		self.kmin = 1e-4
		self.kmax = 40.
		self.h = cosmo.h
		self.nk = 1000
		self.ombh2 = 0.022#cosmo.Ob0 * cosmo.h**2.
		self.omch2 = cosmo.Om0 * cosmo.h**2.
		self.ns = 0.965

		pars = camb.CAMBparams()
		pars.set_cosmology(H0=self.h * 100, ombh2=self.ombh2, omch2=self.omch2)
		pars.set_dark_energy() #re-set defaults
		pars.InitPower.set_params(ns=self.ns)
		self.cambmatterpower = camb.get_matter_power_interpolator(pars, hubble_units = 0, kmax = self.kmax, k_hunit = 0, nonlinear = 0)

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import critical_density [as 别名]
def make_cosmological_sources(exp_time, fov, sky_center, cat_center=None,
                              absorb_model="wabs", nH=0.05, area=40000.0, 
                              output_sources=None, prng=None):
    r"""
    Make an X-ray source made up of contributions from
    galaxy clusters, galaxy groups, and galaxies. 

    Parameters
    ----------
    exp_time : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The exposure time of the observation in seconds.
    fov : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The field of view in arcminutes.
    sky_center : array-like
        The center RA, Dec of the field of view in degrees.
    cat_center : array-like
        The center of the field in the coordinates of the
        halo catalog, which range from -5.0 to 5.0 in 
        degrees in both directions. If None is given, a 
        center will be randomly chosen.
    absorb_model : string, optional
        The absorption model to use, "wabs" or "tbabs". Default: "wabs"
    nH : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional
        The hydrogen column in units of 10**22 atoms/cm**2. 
        Default: 0.05
    area : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional
        The effective area in cm**2. It must be large enough 
        so that a sufficiently large sample is drawn for the 
        ARF. Default: 40000.
    output_sources : string, optional
        If set to a filename, output the properties of the sources
        within the field of view to a file. Default: None
    prng : :class:`~numpy.random.RandomState` object, integer, or None
        A pseudo-random number generator. Typically will only 
        be specified if you have a reason to generate the same 
        set of random numbers, such as for a test. Default is None, 
        which sets the seed based on the system time.
    """
    exp_time = parse_value(exp_time, "s")
    fov = parse_value(fov, "arcmin")
    if nH is not None:
        nH = parse_value(nH, "1.0e22*cm**-2")
    area = parse_value(area, "cm**2")
    prng = parse_prng(prng)
    cosmo = FlatLambdaCDM(H0=100.0*h0, Om0=omega_m)
    agen = ApecGenerator(0.1, 10.0, 10000, broadening=False)

    mylog.info("Creating photons from cosmological sources.")

    mylog.info("Loading halo data from catalog: %s" % halos_cat_file)
    halo_data = h5py.File(halos_cat_file, "r")

    scale = cosmo.kpc_proper_per_arcmin(halo_data["redshift"]).to("Mpc/arcmin")

    # 600. arcmin = 10 degrees (total FOV of catalog = 100 deg^2)
    fov_cat = 10.0*60.0
    w = construct_wcs(*sky_center)

    cat_min = -0.5*fov_cat
    cat_max = 0.5*fov_cat

    if cat_center is None:
        xc, yc = prng.uniform(low=cat_min+0.5*fov, high=cat_max-0.5*fov, size=2)
    else:
        xc, yc = cat_center
        xc *= 60.0
        yc *= 60.0
        xc, yc = np.clip([xc, yc], cat_min+0.5*fov, cat_max-0.5*fov)

    mylog.info("Coordinates of the FOV within the catalog are (%g, %g) deg." %
               (xc/60.0, yc/60.0))

    xlo = (xc-1.1*0.5*fov)*scale.value*h0
    xhi = (xc+1.1*0.5*fov)*scale.value*h0
    ylo = (yc-1.1*0.5*fov)*scale.value*h0
    yhi = (yc+1.1*0.5*fov)*scale.value*h0

    mylog.info("Selecting halos in the FOV.")

    fov_idxs = (halo_data["x"] >= xlo) & (halo_data["x"] <= xhi)
    fov_idxs = (halo_data["y"] >= ylo) & (halo_data["y"] <= yhi) & fov_idxs

    n_halos = fov_idxs.sum()

    mylog.info("Number of halos in the field of view: %d" % n_halos)

    # Now select the specific halos which are in the FOV
    z = halo_data["redshift"][fov_idxs].astype("float64")
    m = halo_data["M500c"][fov_idxs].astype("float64")/h0
    s = scale[fov_idxs].to("Mpc/arcsec").value
    ra0, dec0 = w.wcs_pix2world(halo_data["x"][fov_idxs]/(h0*s)-xc*60.0,
                                halo_data["y"][fov_idxs]/(h0*s)-yc*60.0, 1)

    # Close the halo catalog file
    halo_data.close()

    # Some cosmological stuff
    rho_crit = cosmo.critical_density(z).to("Msun/Mpc**3").value

    # halo temperature and k-corrected flux
#.........这里部分代码省略.........