Python cosmology.FlatLambdaCDM类代码示例

def LogLikelihood(pos, obs, sigmas, par):
    """
        Input:
            Array with the position in the parameters space
            Array with observed values
            Array with sigmas
            Dictionary of indexes of parameters in the arrays
        Return:
            Log likelihood
    """
    mb_obs, x1_obs, c_obs, z_obs = obs
    sigma_mb, sigma_x1, sigma_c = sigmas

    cosmo = FlatLambdaCDM(H0=70, Om0=pos[par["Omega_m"]])
    mu = cosmo.distmod(z_obs).value

    x1 = pos[-2 * len(obs[0]) : -len(obs[0])]
    c = pos[-len(obs[0]) :]

    mb_true = pos[par["MB"]] - pos[par["alpha"]] * x1 + pos[par["beta"]] * c + mu  # +  pos[ par['sigma_int'] ]

    #   I'm not sure this is correct. We need to check how to deal with the two
    #   scatters (sigma_int and sigma_mb)
    likelihood_m_obs = LogGaussian(mb_obs, mb_true, sigma_mb + pos[par["sigma_int"]])
    likelihood_x1 = LogGaussian(x1_obs, x1, sigma_x1)
    likelihood_c = LogGaussian(c_obs, c, sigma_c)

    return np.sum([likelihood_m_obs, likelihood_x1, likelihood_c])

def f(data=data4,ra=ra,dec=dec,freq=freq,label='ICs_10'):
    h=67.77
    cosmo = FlatLambdaCDM(H0=h, Om0=0.307)
    z=1420.4/(freq/10**6)-1
    dc = np.array(cosmo.comoving_distance(z))*h/100        #z: Mpc/h
    N=50
    data=data[N]*10**-3
    freq=freq[N]
    dc=dc[N]
    Hx=(dec[1]-dec[0])/180.*np.pi*dc
    Hy=-(ra[1]-ra[0])/180.*np.pi*dc
    #print Hy
    Lx=Hx*len(dec)
    Ly=Hy*len(ra)
    print np.pi*2/Lx
    freq_x=np.fft.fftfreq(len(dec),1./len(dec))
    freq_y=np.fft.fftfreq(len(ra),1./len(ra))
    deltak=np.fft.fft2(data)
    Pk=np.abs(deltak)**2*(Lx*Ly/len(ra)**2/len(dec)**2)
    k=np.sqrt((2*np.pi/Lx)**2*freq_x[None,:]**2+(2*np.pi/Ly)**2*freq_y[:,None]**2)
    ################################################################################
    bin=10
    edges=np.linspace(k.min(),k.max(),bin+1,endpoint=True)
    n=np.histogram(k,edges)[0]
    print n
    k_bin=np.histogram(k,edges,weights=k)[0]
    print data.shape
    print k.shape
    print Pk.shape
    pk=np.histogram(k,edges,weights=Pk)[0]
    plt.semilogy(k_bin/n,pk/n,'.--',label=label)

def overplot_ruler(ax, z, pixsize=0.396, rlength_arcsec=10., nx=64, ny=64):
    """
    Params
    ------
    ax:
    kpc_per_arcsec
    pixsize=0.396
        in arcsec
    rulerlength_arcsec=10.
        in arcsec
    """
    # import
    from astropy.cosmology import FlatLambdaCDM
    cosmo = FlatLambdaCDM(H0=70, Om0=0.3)

    # set up
    xmid = 0.5*nx+5 # x ending point of bar
    y = ny-10. # y loc of bar

    # conversion
    kpc_per_arcsec = 1./cosmo.arcsec_per_kpc_proper(z)
    rlength_pix = rlength_arcsec/pixsize
    rlength_kpc = (rlength_arcsec*kpc_per_arcsec).value

    #===== plotting
    ax.plot([xmid-rlength_pix, xmid], [y, y], color='white', lw=2)
    ax.plot([xmid-rlength_pix, xmid-rlength_pix], [y+1., y-1.], color='white', lw=2)
    ax.plot([xmid, xmid], [y+1., y-1.], color='white', lw=2)
    ax.text(xmid+4., y+1, '10" ('+'%.0f'%rlength_kpc+' kpc)', color='white', fontsize=12)

def test_approx_DL():
    for z in np.linspace(0.01, 4, num=10):
        z = 2.
        v1 = approx_DL()(z)
        cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=None)
        v2 = cosmo.luminosity_distance(z).value
        assert abs(v1/v2 - 1) < 0.01

def absolute_magnitude_lim(z, app_mag_lim, cosmo=None):
    """
    give the absolute magnitude limit as a function of redshift for a flux-limited survey.
    
    Parameters
    ----------
    z: array_like
        redshift
    
    app_mag_lim: float
       apparent magnitude limit
    
    cosmo: cosmology object
    
    Returns
    -------
    M,z: np.array, np.array
        absolute magnitude in mag+5loh(h) units
    """
    if cosmo==None:
        from astropy.cosmology import FlatLambdaCDM
        cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
        print('Warning, no cosmology specified, using default:',cosmo)
    
    d_L = cosmo.luminosity_distance(z).value
    M = apparent_to_absolute_magnitude(app_mag_lim, d_L)
    
    return M-5.0*np.log10(cosmo.h)

def test_fit_tofz(): 
    z_table, t_table = UT.zt_table()
    
    cosmo = FlatLambdaCDM(H0=70, Om0=0.274)

    prettyplot()
    fig = plt.figure()
    sub = fig.add_subplot(111)

    for deg in range(2,10): 
        coeff = UT.fit_tofz(deg)
        if deg > 5: 
            print 'deg = ', deg, coeff
        tofz = np.poly1d(coeff)
        
        z_arr = np.arange(0., 2., 0.1)
        t_arr = cosmo.age(z_arr).value 

        sub.plot(z_arr, (tofz(z_arr) - t_arr)/t_arr, label='Degree='+str(deg))

    t_of_z = interpolate.interp1d(z_arr, t_arr, kind='cubic') 
    tint = t_of_z(z_table[1:20])#np.interp(t_table[:20], t_arr, z_arr)

    sub.scatter(z_table[1:20], (t_table[1:20] - tint)/tint, c='k', s=30) 
    sub.plot(np.arange(0., 2., 0.1), np.repeat(-0.025, len(np.arange(0., 2., 0.1))), c='k', ls='--', lw=3)
    sub.plot(np.arange(0., 2., 0.1), np.repeat(0.025, len(np.arange(0., 2., 0.1))), c='k', ls='--', lw=3)

    sub.set_ylim([-0.05, 0.05])
    sub.set_xlim([0., 2.])
    sub.legend(loc='upper left')
    plt.show()

    def Grid(self,Nx=64,Ny=64,Nz=256):
        print 'griding',Nx,Ny,Nz
        h=self.h
        ra=self.ra
        ra=ra-ra.mean()
        dec=self.dec
        dec=dec-dec.mean()
        self.Nx=Nx
        self.Ny=Ny
        self.Nz=Nz
        cosmo = FlatLambdaCDM(H0=h, Om0=0.307)
        redshift=1420./(self.freq/10**6)-1
        dc = np.array(cosmo.comoving_distance(redshift))*h/100        #z: Mpc/h
        pos_x=dc[:,None,None]*np.cos(np.pi/180.*dec[None,None,:])*np.cos(np.pi/180.*ra[None,:,None])
        pos_y=dc[:,None,None]*np.cos(np.pi/180.*dec[None,None,:])*np.sin(np.pi/180.*ra[None,:,None])
        pos_z=dc[:,None,None]*np.sin(np.pi/180.*dec[None,None,:])+np.zeros_like(ra[None,:,None])
#       print pos_x.max(),pos_x.min()
#       print pos_y.max(),pos_y.min()
#       print pos_z.max(),pos_z.min()

        #put redshift direction on z axis. and before this, dec>>z ra>>y redshift>>x
        #after transformation, ra>>x dec>>y redshift>>z
        self.DPM=np.c_[pos_y.reshape(-1),pos_z.reshape(-1),pos_x.reshape(-1)]
        self.bin_x=np.linspace(pos_y.min()-10**-5,pos_y.max()+10**-5,Nx+1)
        self.bin_y=np.linspace(pos_z.min()-10**-5,pos_z.max()+10**-5,Ny+1)
        self.bin_z=np.linspace(pos_x.min()-10**-5,pos_x.max()+10**-5,Nz+1)

def get_lum(rmid, zfinal_dict):
    flux = get_raw_lum(rmid)
    z = zfinal_dict[int(rmid)]
    cosmo = FlatLambdaCDM(H0=70.0, Om0=0.3)
    dl_MPC = cosmo.luminosity_distance(z)
    dl_cm = float(dl_MPC.value * 3.085677581 * (10.0 ** 24.0))
    lum = 4.0 * 3.1415926 * dl_cm * dl_cm * flux[0] * 10. ** (0. - 17.) * 5100.0
    return [lum, flux[1] * lum]

def angular_diameter(size,redshift): ## size must be in kpc

    cosmo = FlatLambdaCDM(H0=70, Om0=0.3)  ## defining the cosmology  => see the documentation and details of the parameters and their values
    lumdist=cosmo.luminosity_distance(redshift)  ## value in kpc

    zang=size*(1+redshift)**2/(1000.*lumdist)*180./np.pi*3600. ## arcsec

    return zang

def simple_flux_from_greybody(lambdavector, Trf = None, b = None, Lrf = None, zin = None, ngal = None):
	''' 
	Return flux densities at any wavelength of interest (in the range 1-10000 micron),
	assuming a galaxy (at given redshift) graybody spectral energy distribution (SED),
	with a power law replacing the Wien part of the spectrum to account for the
	variability of dust temperatures within the galaxy. The two different functional
	forms are stitched together by imposing that the two functions and their first
	derivatives coincide. The code contains the nitty-gritty details explicitly.
	Cosmology assumed: H0=70.5, Omega_M=0.274, Omega_L=0.726 (Hinshaw et al. 2009)

	Inputs:
	alphain = spectral index of the power law replacing the Wien part of the spectrum, to account for the variability of dust temperatures within a galaxy [default = 2; see Blain 1999 and Blain et al. 2003]
	betain = spectral index of the emissivity law for the graybody [default = 2; see Hildebrand 1985]
	Trf = rest-frame temperature [in K; default = 20K]
	Lrf = rest-frame FIR bolometric luminosity [in L_sun; default = 10^10]
	zin = galaxy redshift [default = 0.001]
	lambdavector = array of wavelengths of interest [in microns; default = (24, 70, 160, 250, 350, 500)];
	
	AUTHOR:
	Lorenzo Moncelsi [[email protected]]
	
	HISTORY:
	20June2012: created in IDL
	November2015: converted to Python
	'''

	nwv = len(lambdavector)
	nuvector = c * 1.e6 / lambdavector # Hz

	nsed = 1e4
	lambda_mod = loggen(1e3, 8.0, nsed) # microns
	nu_mod = c * 1.e6/lambda_mod # Hz

	#Lorenzo's version had: H0=70.5, Omega_M=0.274, Omega_L=0.726 (Hinshaw et al. 2009)
	cosmo = FlatLambdaCDM(H0 = 70.5 * u.km / u.s / u.Mpc, Om0 = 0.273)
	conversion = 4.0 * np.pi *(1.0E-13 * cosmo.luminosity_distance(zin) * 3.08568025E22)**2.0 / L_sun # 4 * pi * D_L^2    units are L_sun/(Jy x Hz)

	Lir = Lrf / conversion # Jy x Hz

	Ain = np.zeros(ngal) + 1.0e-36 #good starting parameter
	betain =  np.zeros(ngal) + b 
	alphain=  np.zeros(ngal) + 2.0

	fit_params = Parameters()
	fit_params.add('Ain', value= Ain)
	#fit_params.add('Tin', value= Trf/(1.+zin), vary = False)
	#fit_params.add('betain', value= b, vary = False)
	#fit_params.add('alphain', value= alphain, vary = False)

	#pdb.set_trace()
	#THE LM FIT IS HERE
	#Pfin = minimize(sedint, fit_params, args=(nu_mod,Lir.value,ngal))
	Pfin = minimize(sedint, fit_params, args=(nu_mod,Lir.value,ngal,Trf/(1.+zin),b,alphain))

	#pdb.set_trace()
	flux_mJy=sed(Pfin.params,nuvector,ngal,Trf/(1.+zin),b,alphain)

	return flux_mJy

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))

def calculate_rf_lens_magnitudes(lens_redshift, velocity_dispersion, filters_dict):
    """Calculates the reference-frame lens magnitudes in multiple filters

    Parameters
    ----------
    lens_redshift : float
        Redshift of the lens
    velocity_dispersion : float
        Velocity dispersion of the lens
    filters_dict : dict
        (See output of `get_sdss_filters` for details)
        Throughputs of various filters

    Returns
    -------
    dict
        Each key is one of string characters 'u', 'g', 'r', 'i', 'z'
        representing the filter
        Each value is the reference-frame apparent magnitude of the quasar
        in the 'key' filter, of type float
    """
    from stellarpop import tools
    from lenspop import population_functions, distances
    from astropy.cosmology import FlatLambdaCDM

    # Instantiate Distance
    distance = distances.Distance() #TODO: necessary?
    # Instantiate LensPopulation
    lenspop = population_functions.LensPopulation_()
    # Instantiate FlatLambdaCDM cosmology with reasonable parameters
    cosmology = FlatLambdaCDM(H0=70.0, Om0=0.3)

    lens_sed = tools.getSED('BC_Z=1.0_age=9.000gyr')
    velocity_dispersion = np.atleast_1d(velocity_dispersion)

    # Absolute --> apparent magnitude conversion in the R-band
    lens_abmag_r = tools.ABFilterMagnitude(filters_dict['r'], lens_sed, lens_redshift)
    distance_modulus = cosmology.distmod(lens_redshift).value
    lens_appmag_r = lens_abmag_r + distance_modulus

    # [Reference frame] Absolute --> apparent magnitude conversion in the R-band
    rf_lens_abmag_r, _ = lenspop.EarlyTypeRelations(velocity_dispersion)
    rf_lens_appmag = {}
    rf_lens_appmag['r'] = rf_lens_abmag_r + distance_modulus

    # Quantity which is added to ~magnitude to convert it into reference-frame ~magnitude
    offset_rf = rf_lens_abmag_r - lens_abmag_r

    # Converting absolute magnitude to reference-frame apparent magnitude
    for band in 'ugiz':
        rf_lens_appmag[band] = tools.ABFilterMagnitude(filters_dict[band], lens_sed, lens_redshift) + offset_rf + distance_modulus
    return rf_lens_appmag

def L_nu_from_magAB(magAB = None, z = None):
    import numpy as np
    
    """Convert absolute magnitude into luminosity (erg s^-1 Hz^-1)."""
    
    from astropy.cosmology import FlatLambdaCDM
    cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
    
    m_zero = 30.
    d = cosmo.comoving_distance( z ).value
    factor = 4. * np.pi * d**2.
    L_nu =  factor * 10.**((magAB - m_zero)/(-2.5))
    
    return L_nu

def test_fit_cosmology():
    """Test fitting cosmology on simulated data."""

    # Generate some fake data.
    cosmo = FlatLambdaCDM(H0=70., Om0=0.25)
    z = np.random.rand(200)
    mb = -19.3 + cosmo.distmod(z).value
    mberr = 0.2 * np.ones_like(z)

    # fit to fake data
    fitted_cosmo = fitting.fit_cosmology(z, mb, mberr)

    # check that fitted H0 value is same as input
    assert_allclose(cosmo.H0.value, fitted_cosmo.H0.value, rtol=1e-4)

 def Pos(self,freq=None,ra=None,dec=None):
     # freq ->> x axis; ra ->> y axis; dec ->> z axis
     self.shape=(freq.shape[0],ra.shape[0],dec.shape[0])
     cosmo = FlatLambdaCDM(H0=self.h, Om0=self.Om0)
     redshift=1420./(freq/10**6)-1
     dc = np.array(cosmo.comoving_distance(redshift))*self.h/100        #z: Mpc/h
     self.pos_z=dc[:,None,None]*np.sin(dec[None,None,:]/180.*np.pi)*np.ones(self.shape)
     self.pos_y=dc[:,None,None]*np.sin(ra[None,:,None]/180.*np.pi)*np.ones(self.shape)
     self.pos_x=dc[:,None,None]*np.ones(self.shape)
     DPM=np.c_[self.pos_x.reshape(-1),self.pos_y.reshape(-1),self.pos_z.reshape(-1)]
     bin_x=np.linspace(self.pos_x.min()-10**-5,self.pos_x.max()+10**-5,self.Nx+1)
     bin_y=np.linspace(self.pos_y.min()-10**-5,self.pos_y.max()+10**-5,self.Ny+1)
     bin_z=np.linspace(self.pos_z.min()-10**-5,self.pos_z.max()+10**-5,self.Nz+1)
     return DPM,bin_x,bin_y,bin_z