Python FlatLambdaCDM.angular_diameter_distance方法代码示例

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
def DistanceFraction(H0, Om, z_gals, z_clust):
    Cos = FlatLambdaCDM(H0=H0, Om0=Om)
    Ds = Cos.angular_diameter_distance(z_gals)
    Dl = Cos.angular_diameter_distance(z_clust)
    
    #Calculate Angular Diameter Distance between objects and Cluster
    DM1 = Cos.comoving_distance(z_clust)
    DM2 = Cos.comoving_distance(z_gals)
    Dls = (DM2 - DM1)/(1 + z_gals)
    return NP.array(Ds/(Dls*Dl))

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
def find_nearest_in_Mpc(sam, ref, cosmo=None):
    if cosmo is None:
        from astropy.cosmology import FlatLambdaCDM
        cosmo = FlatLambdaCDM(70, 0.3)
    ang_diam_dists = cosmo.angular_diameter_distance(ref['z'])
    ref_idx = []
    sep_Mpc = []
    for ra, dec in sam['ra', 'dec']:
        seps = angsep(ref['ra'], ref['dec'], ra, dec, sepunits='radian')
        seps *= ang_diam_dists.value
        sep_Mpc.append(seps.min())
        ref_idx.append(seps.argmin())
    ref_idx = np.array(ref_idx)
    sep_Mpc = np.array(sep_Mpc)
    return ref_idx, sep_Mpc

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
class Source(object):
    
    def __init__(self, Zs):  
    
        self.Zs = Zs
        self.compute_distances()
        self.setup_grid(NX=100,NY=100,pixscale=0.1)
    
        return    
# ----------------------------------------------------------------------
    def compute_distances(self):
        self.cosmological = FlatLambdaCDM(H0=71.0, Om0=0.2669)
        self.Ds = self.cosmological.angular_diameter_distance(self.Zs)
# ----------------------------------------------------------------------   
    def read_source_from(self, fitsfile):
        '''
        Read an image from a fitsfile, and setup its grid.  Here we need
        to properly read in the wcs coordinate information so that the 
        grid is spaced correctly. This means we have to extract the pixel 
        scale from the FITS header.
        '''
        if fitsfile is None:
            raise Exception("You need to provide an image.\n")
            
        hdulist = fits.open(fitsfile)
        self.hdr = hdulist[0].header
        self.intensity = hdulist[0].data
        hdulist.close()
              
        if self.hdr['NAXIS'] == 2:
            if self.intensity.shape == (self.hdr['NAXIS1'],self.hdr['NAXIS2']):
                self.NX,self.NY = self.intensity.shape
            elif self.intensity.shape ==(self.hdr['NAXIS2'],self.hdr['NAXIS1']):
                self.NY,self.NX = self.intensity.shape
            else:
                raise Exception("Your image is formatted incorrectly.\n")    
        else:
            assert len(self.intensity.shape) == 3
            if self.intensity.shape == (self.hdr['NAXIS'],self.hdr['NAXIS1'],self.hdr['NAXIS2']):
                self.Naxes,self.NX,self.NY = self.intensity.shape
            elif self.intensity.shape ==(self.hdr['NAXIS'],self.hdr['NAXIS2'],self.hdr['NAXIS1']):
                self.Naxes,self.NY,self.NX = self.intensity.shape
            else:
                raise Exception("Your image is formatted incorrectly.\n")
        self.set_pixscale()
        
        # Set up a new pixel grid to go with this new kappa map:
        self.setup_grid()
        
        return
# ----------------------------------------------------------------------
    
    def setup_grid(self, NX=None, NY=None, pixscale=None):
        '''
        Make two arrays, x and y, that define the extent of the maps
        - pixscale is the size of a pixel, in arcsec.
        - 
        '''
        if NX is not None: 
            self.NX = NX
        if NY is not None: 
            self.NY = NY
        if pixscale is not None: 
            self.pixscale = pixscale        
        
        xgrid = np.arange(-self.NX/2.0,(self.NX)/2.0,1.0)*self.pixscale+self.pixscale
        ygrid = np.arange(-self.NY/2.0,(self.NY)/2.0,1.0)*self.pixscale+self.pixscale
        
        self.beta_x, self.beta_y = np.meshgrid(xgrid,ygrid)        
        
        return

# ----------------------------------------------------------------------
    
    def set_pixscale(self):
        
        # Modern FITS files:
        if 'CD1_1' in self.hdr.keys():            
            determinant = self.hdr['CD1_1']*self.hdr['CD2_2'] \
                          - self.hdr['CD1_2']*self.hdr['CD2_1']
            self.pixscale = 3600.0*np.sqrt(np.abs(determinant))

        # Older FITS files:
        elif 'CDELT1' in self.hdr.keys():
            self.pixscale = 3600.0*np.sqrt(np.abs(self.hdr['CDELT1']*self.hdr['CDELT2']))

        # Simple FITS files with no WCS information (bad):
        else:
            self.pixscale = 1.0
            
        return        

# ----------------------------------------------------------------------
    def build_from_clumps(self,size=2.0,clump_size = 0.1,axis_ratio=1.0, orientation=0.0,center=[0,0], Nclumps=50, n = 1 , error =10**-8,singlesource=False,seeds=[1,2,3],Flux=1.0):
        #raise Exception("cannot build source from clumps yet. \n")
        '''
        Build source from gaussian clumps centered about specified position.
        Accepted parameters are as follows:
        - Size of the source, in kpc (approximately the half light radius)
        - Size of individual clumps.  RMS determines spread in clump size.
#.........这里部分代码省略.........

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
    def __init__(self, zmin=0.5, zmax=7.0, Om0=0.315, H0=67.7, phib0=-3.02,
                 gamma_sfmf=0.4, ninterp=1000):
        """ Initializer.

        Parameters
        ----------
        zmin : float
          Minimum redshift

        zmax: float
          Maximum redshift

        Om0: float
          Matter density parameter

        H0: float
          Hubble constant in km / s / Mpc

        phib0: float
          log 10 number density at SFMF break in comoving Mpc^-3

        gamma_sfmf: float
          Evolution of the number density of the SFMF at z > 1

        ninterp: int
          Number of interpolation samples to use

        """

        from scipy.interpolate import interp1d
        from scipy.integrate import trapz
        from astropy.cosmology import FlatLambdaCDM
        from astropy.units import Quantity

        self._zmin = float(zmin)
        self._zmax = float(zmax)
        if self._zmin == self._zmax:
            raise ValueError("No range between zmin and zmax")
        if self._zmin > self._zmax:
            self._zmin, self._zmax = self._zmax, self._zmin
        if self._zmin < 0.0:
            raise ValueError("zmin must be >= 0: {:f}".format(self._zmin))
        self._Om0 = float(Om0)
        if self._Om0 <= 0.0:
            raise ValueError("Om0 must be positive: {:f}".format(self._Om0))
        self._H0 = float(H0)
        if self._H0 <= 0.0:
            raise ValueError("H0 must be positive: {:f}".format(self._H0))
        self._ninterp = int(ninterp)
        if self._ninterp <= 0:
            raise ValueError("Ninterp must be > 0: {:d}".format(self._ninterp))
        self._phib0 = float(phib0)
        self._gamma_sfmf = float(gamma_sfmf)

        self._zvals = np.linspace(self._zmin, self._zmax, self._ninterp)

        # Cosmology bit
        c_over_H0 = 299792.458 / self._H0  # in Mpc
        cos = FlatLambdaCDM(H0=self._H0, Om0=self._Om0)

        # in comoving Mpc^3
        ang_diam_dist = cos.angular_diameter_distance(self._zvals)
        if isinstance(ang_diam_dist, Quantity):
            ang_diam_dist = ang_diam_dist.value
        self._dVdzdOmega = c_over_H0 * (1.0 + self._zvals)**2 * \
            ang_diam_dist**2 / np.sqrt((1.0 + self._zvals)**3 *
                                       self._Om0 + (1.0 - self._Om0))

        # Schecter evolution bit
        phi = self._phib0 * np.ones(self._ninterp, dtype=np.float64)
        wgt1 = np.nonzero(self._zvals > 1.0)[0]
        if len(wgt1) > 0:
            phi[wgt1] += self._gamma_sfmf * (1.0 - self._zvals[wgt1])

        # Combined
        comb = 10**phi * self._dVdzdOmega

        # Needed to understand normalization
        self._dVPhidzdOmega = trapz(comb, x=self._zvals)

        # Form inverse cumulative array needed to generate samples
        cumsum = comb.cumsum()
        cumsum -= cumsum[0]  # So that 0 corresponds to the bottom
        cumsum /= cumsum[-1]  # Normalization -> 0-1 is full range
        self._interpolant = interp1d(cumsum, self._zvals, kind='linear')

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
wig_post_daonh = wig_post_da / wig_post_h
wig_post_daonh_error = np.sqrt((wig_post_da_error/wig_post_da)**2 + (wig_post_h_error/wig_post_h)**2) * wig_post_daonh
wig_post_dadt = wig_post_h / (67.74 * (1 + wig_post_zs))
wig_post_dadt_error = (wig_post_h_error / wig_post_h) * wig_post_dadt

planck_om = 0.3089
planck_om_error = 0.0062
planck_h = 67.74
planck_h_error = 0.46


top_cosmology = FlatLambdaCDM(planck_h + planck_h_error, planck_om + planck_om_error)
bottom_cosmology = FlatLambdaCDM(planck_h - planck_h_error, planck_om - planck_om_error)
zs = np.linspace(0.2, 0.9, 100)
das = np.linspace(1000, 1500, 1000)
top_da = top_cosmology.angular_diameter_distance(zs).value
bottom_da = bottom_cosmology.angular_diameter_distance(zs).value
top_h = top_cosmology.H(zs).value
bottom_h = bottom_cosmology.H(zs).value
top_div = top_da / top_h
bottom_div = bottom_da / bottom_h
hconst = 67.74 * (1 + zs)
daconst = 3e5 / (67.74 * (1 + zs)) * np.log(1 + zs)
hup = top_h / (67.74 * (1 + zs))
hdown = bottom_h / (67.74 * (1 + zs))


def clamp(val, minimum=0, maximum=255):  # pragma: no cover
    if val < minimum:
        return minimum
    if val > maximum:

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
# hopefully that the arrays shape will be broadcast
N_halos = [2] * no
centroid_steps = [1. / np.sqrt(2) / 60. / 60.] * no
cent_new = map(get_new_centroids, N_halos, centroid_old, centroid_steps)
cent_new = np.array(cent_new)
[NW_ra, NW_dec, NW_z, SE_ra, SE_dec, SE_z] = cent_new.transpose()

# inputs of angular separation has to be in units of radians
ang_sep = \
    ang_util.angular_separation(NW_ra / 180. * np.pi,
                                NW_dec / 180. * np.pi,
                                SE_ra / 180. * np.pi,
                                SE_dec / 180. * np.pi)
cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
D_proj = cosmo.angular_diameter_distance(0.87) * ang_sep


###
# Debug input
###

def plot_dproj():
    plt.title('D-proj')
    plt.xlabel('Mpc')
    plt.hist(D_proj, bins=100, normed=True, histtype='step')
    plt.savefig(out_prefix + '_D_proj.png')
    plt.close()


def plot_hist_m_main():

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
import astropy.units as u
import astropy.constants as const
from astropy.cosmology import FlatLambdaCDM

G = const.G.to(u.Mpc**3/u.Gyr**2/(1.e10*u.solMass)).value
sig0 = (1.e-3 *u.g/u.cm**2).to(1.e10*u.solMass/u.Mpc**2).value
pram0 = (1.e-11 * u.erg/u.cm**3).to(1.e10*u.solMass/u.Mpc/u.Gyr**2).value

nH_ne = 0.852

H0 = 70.
cosmo = FlatLambdaCDM(H0=H0, Om0=0.3)

DA = lambda z: cosmo.angular_diameter_distance(z).to(u.cm).value
# Return proper distance [Mpc] for a given redshift and sep. in arcsec
properDist = lambda arcsec, z: arcsec*cosmo.kpc_proper_arcmin(z)/60./1000.

def EI(norm, z):
  return 10 ** (np.log10(norm) + np.log10(4*np.pi) + 2*np.log10( DA(z)*(1 + z) )  + 14)

def R_delta(delta, z, kT, beta_T=1.05):
  delta_z = (delta * cosmo.Om0) / (18*np.pi**2 * cosmo.Om(z)) 
  return 3.8 * np.sqrt(beta_T* (kT/10.) / (delta_z*(1+z)**3) ) / (H0/50.)

def calcLum(flux, z):
  return 4*np.pi*10**(2*np.log10(cosmo.luminosity_distance(z).to(u.cm).value) + np.log10(flux) )

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
class Calculator:
    """
    Calculate various quantities under a specific cosmology, as well
    as some values with respect to Chandra ACIS detector properties.
    """
    def __init__(self, H0=71.0, Om0=0.27, Ob0=0.046):
        self.H0 = H0  # [km/s/Mpc]
        self.Om0 = Om0
        self.Ode0 = 1.0 - Om0
        self.Ob0 = Ob0
        self._cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Ob0=Ob0)

    def evolution_factor(self, z):
        return self._cosmo.efunc(z)

    def luminosity_distance(self, z, unit="Mpc"):
        dist = self._cosmo.luminosity_distance(z)
        return dist.to(au.Unit(unit)).value

    def angular_diameter_distance(self, z, unit="Mpc"):
        dist = self._cosmo.angular_diameter_distance(z)
        return dist.to(au.Unit(unit)).value

    def kpc_per_arcsec(self, z):
        """
        Calculate the transversal length (unit: kpc) corresponding to
        1 arcsec at the *angular diameter distance* of z.
        """
        dist_kpc = self.angular_diameter_distance(z, unit="kpc")
        return dist_kpc * au.arcsec.to(au.rad)

    def kpc_per_pix(self, z):
        """
        Calculate the transversal length (unit: kpc) corresponding to
        1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance*
        of z.
        """
        pix = ACIS.pixel2arcsec * au.arcsec
        dist_kpc = self.angular_diameter_distance(z, unit="kpc")
        return dist_kpc * pix.to(au.rad).value

    def cm_per_pix(self, z):
        """
        Calculate the transversal length (unit: cm) corresponding to
        1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance*
        of z.
        """
        return self.kpc_per_pix(z) * au.kpc.to(au.cm)

    def norm_apec(self, z):
        """
        The normalization factor of the XSPEC APEC model assuming
        EM = 1 (i.e., n_e = n_H = 1 cm^-3, and V = 1 cm^3)

        norm = 1e-14 / (4*pi* (D_A * (1+z))^2) * int(n_e * n_H) dV
        unit: [cm^-5]

        This value will be used to calculate the cooling function values.

        References
        ----------
        * XSPEC: APEC model
          https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSmodelApec.html
        """
        da = self.angular_diameter_distance(z, unit="cm")
        norm = 1e-14 / (4*math.pi * (da * (1+z))**2)
        return norm

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import angular_diameter_distance [as 别名]
def _LS(lense, source, rand, weight = None, Boost = False):
    """
        Calculate Lensing Signal deltaSigma, Boost Factor, Corrected deltaSigma
        
        parameter
        ---------
        lense : lense catalog
        source : source catalog (shear catalog that contains e1, e2)
        rand : random catalog
        
        """
    
    import treecorr
    
    weight_lense, weight_source, weight_rand = weight
    
    z_s_min = 0.7
    z_s_max = 1.0
    
    num_lense_bin = 3
    num_source_bin = 3
    z_l_bins = np.linspace(0.45, 0.55, num_lense_bin)
    z_s_bins = np.linspace(z_s_min, z_s_max, num_source_bin)
    #matrix1, matrix2 = np.mgrid[0:z_l_bins.size, 0:z_s_bins.size]
    source = source[ (source['DESDM_ZP'] > z_s_min ) & (source['DESDM_ZP'] < z_s_max )]
    lense = lense[ (lense['DESDM_ZP'] > 0.45) & (lense['DESDM_ZP'] < 0.55) ]
    rand = rand[ (rand['Z'] > 0.45) & (rand['Z'] < 0.55) ]
    
    z_l_bincenter, lense_binned_cat,_ = divide_bins( lense, Tag = 'DESDM_ZP', min = 0.45, max = 0.55, bin_num = num_lense_bin)
    z_r_bincenter, rand_binned_cat,_ = divide_bins( rand, Tag = 'Z', min = 0.45, max = 0.55, bin_num = num_lense_bin)
    z_s_bincenter, source_binned_cat,_ = divide_bins( source, Tag = 'DESDM_ZP', min = z_s_min, max = z_s_max, bin_num = num_source_bin)
    
    # angular diameter
    from astropy.cosmology import FlatLambdaCDM
    cosmo = FlatLambdaCDM(H0=70,Om0=0.274)
    
    from astropy import constants as const
    from astropy import units as u
    
    h = 0.7
    c = const.c.to(u.megaparsec / u.s)
    G = const.G.to(u.megaparsec**3 / (u.M_sun * u.s**2))
    
    dA_l = cosmo.angular_diameter_distance(z_l_bins) * h
    dA_s = cosmo.angular_diameter_distance(z_s_bins) * h
    matrix1, matrix2 = np.mgrid[0:z_l_bins.size, 0:z_s_bins.size]
    dA_ls = cosmo.angular_diameter_distance_z1z2(z_l_bins[matrix1],z_s_bins[matrix2]) * h
    dA_l_matrix = dA_l[matrix1]
    dA_s_matrix = dA_s[matrix2]
    Sigma_critical = c**2 / (4 * np.pi * G) * dA_s_matrix/(dA_l_matrix * dA_ls)
    
    
    
    
    
    min_sep_1 = 0.001
    bin_size = 0.4
    nbins = 30
    
    theta = [min_sep_1 * np.exp(bin_size * i) for i in range(nbins) ]
    r_p_bins = theta * dA_l[0]
    
    
    if Boost is True:
        n_SL = []
        n_SR = []
        N_SL = []
        N_SR = []
        
        weight = 1./(Sigma_critical)**2
        weight = np.ones(Sigma_critical.shape)
        
        print " **  To do : add Boost codes Weight "
        
        for i, S in enumerate(source_binned_cat):
            source_cat = treecorr.Catalog(ra=S['RA'], dec=S['DEC'], ra_units='deg', dec_units='deg')
            for (j, dA), L, R in zip( enumerate(dA_l), lense_binned_cat, rand_binned_cat ):
                min_sep = theta[0] * dA_l[0]/dA
                lense_cat = treecorr.Catalog(ra=L['RA'], dec=L['DEC'], ra_units='deg', dec_units='deg')
                rand_cat = treecorr.Catalog(ra=R['RA'], dec=R['DEC'], ra_units='deg', dec_units='deg')
                SL = treecorr.NNCorrelation(min_sep=min_sep, bin_size=bin_size, nbins = nbins, sep_units='degree')
                SR = treecorr.NNCorrelation(min_sep=min_sep, bin_size=bin_size, nbins = nbins, sep_units='degree')
                SL.process(lense_cat, source_cat)
                SR.process(rand_cat, source_cat)
                w = 1./weight[j,i]**2
                n_SL.append(w * SL.npairs/np.sum(SL.npairs))
                n_SR.append(w * SR.npairs/np.sum(SR.npairs))
                N_SL.append(SL.npairs)
                N_SR.append(SR.npairs)
        
            N_SL_list = np.sum( N_SL, axis = 0)
            N_SR_list = np.sum( N_SR, axis = 0)
            n_SL_list = np.sum( n_SL, axis = 0)
            n_SR_list = np.sum( n_SR, axis = 0)
            
            BoostFactor = n_SL_list/n_SR_list
            
            nbar = 1e-4
            P = 1e+4
            w_FKP = 1./(nbar*P + 1)
#.........这里部分代码省略.........