Python FlatLambdaCDM.comoving_distance方法代码示例

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_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 comoving_distance [as 别名]
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)

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
    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)

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
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

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
 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

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
def fill_cone(ralim,declim,zlim,xyz=None,b=None,repmethod='rotation',
              npts=80000,nvoids=2000,nstep=100,rmin=None,cosmo=None):
    '''
    Fill the volume of an observation light cone with small cubes of 3D points.

    The cube is stacked in 3D along the desired volume with as many cubes as 
    needed. This can generate artificial strucure so each cube can be : 

    (1) Copied as is (``repmethod='copy'``). This will likely introduce artificial patterns.
    
    (2) Randomly rotated and mirrored (``repmethod='rotation'``)
    
    (3) Generated randomly, so each box is uniquely random (``repmethod='fullrandom'``). Naturally this is slower
    
    In case (3), the routine ``filament_box()`` is used to generate boxes
    with filamentary structure

    Parameters
    ----------
    ralim : list (ra0, ra1)
        RA limits [deg]
    declim : list (dec0, dec1)
        DEC limits [deg]
    zlim : list (z0, z1)
        Redshift limits
    xyz : array of shape (:,3)
        The building block of points to be stacked. This is **not used** when
        ``repmethod=fullrandom``, as then each box is generated uniquely
    b : float
        Box size in Mpc. Default 150
    repmethod : string
        * 'copy' : the same box is repeated
        * 'rotation' : the same box is repeated but rotated and mirrored. Default
        * 'fullrandom' : generate unique random boxes each time
    npts : integer
        Number of points in the initial unit sphere. Default 80000. Used only 
        if ``repmethod='fullrandom'``
    nvoids : integer
        Number of voids in the initial unit sphere. Default 2000. Used only 
        if ``repmethod='fullrandom'``
    nsteps : integer
        Number of steps that points move away from void. In general, less steps
        means softer and weaker filaments. Default 100. Used only if 
        ``repmethod='fullrandom'``
    rmin : float
        After box is created, remove 1 member of each pair closer than rmin [Mpc]. 
        Used only if ``repmethod='fullrandom'``
    cosmo : astropy cosmology object
        If not given, defaults to a flat LambdaCDM with H_0=100 and Om=0.3
        
    Notes
    -----
    The box is always unit size, which is then scaled to any meaningfull size
    by ``b``. The units depend on your intepretation. A choice of ``nvoids=2000`` 
    and ``b=150`` gets a box of 300 Mpc with voids and filaments of size similar 
    to real ones. In general :
    
    * For a fixed ``b``, the larger ``nsteps`` the stronger the features
    
    * For a fixed ``b``, The larger ``nvoids`` the smaller the structure
    
    * Smaller structure means more points can be included to trace small scales
      without getting unrealistic features.
    
    It is hard to estimate a priori how many points will fall inside the
    observation cone. Just give a try and change ``npts`` accordingly.

    Returns
    ----------
    kone : array of shape (n,7)
        Array of 7 columns [ra,dec,redshift,comdis(redshift),x,y,z]
    '''
    print '-------------------------------------------------------------'
    print 'Filling light cone with cubes'

    if cosmo is None: cosmo = FlatLambdaCDM(H0=100, Om0=0.3)

    # Get (ra/dec/z) into units of radians and comov_Mpc  ---------------------
    torad = np.pi/180
    ra0, ra1  = ralim[0]*torad, ralim[1]*torad
    dec0,dec1 = declim[0]*torad, declim[1]*torad
    d0        = cosmo.comoving_distance(zlim[0]).value
    d1        = cosmo.comoving_distance(zlim[1]).value

    d_max = np.floor(1 + d1/b)*b  # maximum comov distance in box lengths
    xyza = np.zeros([1,3]) # accumulation array with dummy zero element

    # Loop for a 3d grid of cubes coverting the cone length  ------------------
    nb = 0
    for x0 in np.arange(-1*d_max,d_max,b):
        for y0 in np.arange(-1*d_max,d_max,b):
            for z0 in np.arange(-1*d_max,d_max,b):
                x1 = x0 + b
                y1 = y0 + b
                z1 = z0 + b
                # For each cube find corners (0 and 1) in ra/dec/dist/x/y/z space
                if (x0>=0) and (y0>=0):
                    tra0 = np.arctan2(y0,x1)
                    tra1 = np.arctan2(y1,x0)
                    if (z0>=0):
#.........这里部分代码省略.........

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]

#.........这里部分代码省略.........
        uu = np.random.randint(4,nnotunif+1,ncen-1)
        tt = np.append(uu,[4,nnotunif+1])
        tt.sort()
        Cmem = np.diff(tt)  # n_members by some proportionality to extension

    # Generate clusters   -----------------------------------------------------
    print 'Generating clusters...'
    for i in range(ncen):
        print '  |--',i,Cra[i],Cdec[i],Cred[i],Crad[i],Cmem[i]
        # Real space center for the 3D gaussian (center of cluster)
        mm = np.array([Cx[i],Cy[i],Cz[i]])
        # Real space std for the 3D gaussian (sort of size of cluster)
        if rand_elong : 
            efac = np.random.random(3)*3.
        else :
            efac = np.ones(3)
        stdx = Crad[i]*efac[0]
        stdy = Crad[i]*efac[1]
        stdz = Crad[i]*efac[2]
        cv = np.array([[stdx,0.,0.],[0.,stdy,0.],[0.,0.,stdz]])
        # Number of members of the cluster
        #nmem = np.int(1.0*nnotunif/ncen)  #same nr for all clusters
        nmem = Cmem[i]
        # Finally generate the cluster
        if Cmem[i] > 0:
            m=np.random.multivariate_normal(mm,cv,nmem)
            if i == 0 : 
                abc = m
            else:
                abc = np.concatenate((abc,m),axis=0)
    x1,y1,z1 = abc[:,0],abc[:,1],abc[:,2]
            
    # Convert clusters from (x,y,z) to (ra,dec,redshift)  ---------------------
    ttt = cartesian_to_spherical(x1,y1,z1)
    comd1,dec1,ra1 = ttt[0].value,ttt[1].value*180./np.pi,ttt[2].value*180./np.pi
    tmpz = np.linspace(0.,zlim[1]+7.,500)
    tmpd = cosmo.comoving_distance(tmpz).value
    reds1 = np.interp(comd1,tmpd,tmpz,left=-1., right=-1.)
    if (reds1<0).any(): raise Exception('Problem with redshifts!')

    # Cut clusters to desired ra,dec,redshift window   ------------------------
    idx,=np.where( (ra1>ralim[0]) & (ra1<ralim[1]) & (dec1>declim[0]) & (dec1<declim[1]) & (reds1>zlim[0]) & (reds1<zlim[1]) )
    ra1,dec1,reds1=ra1[idx],dec1[idx],reds1[idx]
    ncpop=len(ra1)  # this is the effective nr of clustered objects in the cone
    
    # Add uniform background points to complete n objects  --------------------
    nback = n-ncen-ncpop
    ra2,dec2 = uniform_sky(ralim,declim,n=nback)
    reds2    = (zlim[1]-zlim[0])*np.random.random(nback) + zlim[0]
    x2,y2,z2 = rdz2xyz(ra2,dec2,reds2,cosmo)
    #print 'n,nunif,nnotunif,ncen,ncpop,nback',n,nunif,nnotunif,ncen,ncpop,nback
    print 'Results'
    print '  |--Nr of clusters              :',ncen
    print '  |--Nr of cluster members       :',ncpop
    print '  |--Nr of non-cluster objects   :',nback
    print '  |--Total nr of objects in cone :',n

    # Join everything and do plot  --------------------------------------------
    ra   = np.concatenate([Cra,ra1,ra2])
    dec  = np.concatenate([Cdec,dec1,dec2])
    reds = np.concatenate([Cred,reds1,reds2])
    distC = cosmo.comoving_distance(Cred).value
    dist1 = cosmo.comoving_distance(reds1).value
    dist2 = cosmo.comoving_distance(reds2).value
    dist   = np.concatenate([distC,dist1,dist2])
    x = np.concatenate([Cx,x1,x2])
    y = np.concatenate([Cy,y1,y2])
    z = np.concatenate([Cz,z1,z2])

    if colorize:
        c2,c1,cc = 'k','b','r'
        spts, scen = 0.2, 20
    else:
        c2,c1,cc = 'k','k','k'
        spts, scen = 0.2, 0.2

    if doplot:
        #x,y=raz2xy(ra,z)
        #plt.scatter(x,y,s=0.4)  #c=z,cmap='jet'
        #plt.scatter(ra,dec,s=spts,c=z)
        x2,y2=raz2xy(ra2,reds2)
        x1,y1=raz2xy(ra1,reds1)
        xC,yC=raz2xy(Cra,Cred)
        plt.scatter(x2,y2,s=spts,color=c2)
        plt.scatter(x1,y1,s=spts,color=c1)
        plt.scatter(xC,yC,s=scen,color=cc)
        
        plt.figure()
        plt.scatter(ra2,dec2,s=spts,color=c2)
        plt.scatter(ra1,dec1,s=spts,color=c1)
        plt.scatter(Cra,Cdec,s=scen,color=cc)
        plt.axis('equal')

    # Choose ouput format  ----------------------------------------------------    
    kone = np.asarray(zip(ra,dec,reds,dist,x,y,z))
    if oformat == 'table' :
        cols = ['ra','dec','z','comd','px','py','pz']
        kone = Table(data=kone, names=cols)

    return kone

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
############ verifier que les X_, Y_, Z_ sont bien en Mpc/h!!!!!


### X_cone, Y_cone and Z_cone
print 'deriving points coordinates...'
cells_points_Mpc = (cells.points-center) * info["unit_length"].express(C.Mpc)

### conversion boite proper ou comoving?



### luminosity distance
print 'deriving comoving distances...'

## The redshift corresponds to the center point of the simulation box
X_cone_center = cosmo.comoving_distance(z_center).value ## Mpc



#AMR_position = (xcell,ycell,zcell)
print 'deriving AMR coordinates...'

### passer en matriciel
X_cone = np.array([np.dot(los,a) for a in cells_points_Mpc]) + np.array(X_cone_center)
Y_cone = np.array([np.dot(RA,a) for a in cells_points_Mpc])
Z_cone =np.array([np.dot(DEC,a) for a in cells_points_Mpc])

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
pool.close()

mag_app = mag_app[0]

if Make_LF:
	nfilters = 0
	for j in np.arange(nFil):
		id_filter_i = np.where((id_filter == j) & 
					  (~np.isinf(mag_app)))
		gals_filt = len(id_filter_i[0])
		if gals_filt == 0:
			continue
		
		z_min = p10p90W[j,0]/L_line - 1
		z_max = p10p90W[j,1]/L_line - 1
		d1 = cosmo.comoving_distance(z_max)
		d0 = cosmo.comoving_distance(z_min)
		Vol = 4./3*np.pi*(d1.value**3 - d0.value**3)/8.0

		mag_sel = np.zeros(len(id_filter_i[0]))
		for i in range(len(id_filter_i[0])): mag_sel[i] = mag_app[id_filter_i[0][i]]
#		mag_sel = mag_app[0][id_filter_i
#		allmag_lf[j] = get_lf(mag_app[0][id_filter_i],
		allmag_lf[j] = get_lf(list(mag_sel),
		Vol = Vol, binsize=lfbin,minsample = lfmin,maxsample=lfmax)
		print 'LF computed for filter ',j
		nfilters += 1


if Make_Figure:
	pl.figure()

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
			if ff==1:
				#20180109: add offsets correction here
				MODEL_FFT = np.fft.fft2( STACKED_KAPPA_MODEL )

				fcen = 0.78 #page 14 of https://arxiv.org/pdf/1601.00621.pdf (Rykoff 2016)
				fmis = 1 - fcen
				D_a_sigma_c = 0.42/param_dict['h'] #0.42h-1 Mpc from 1707.09369
				CLUS_REDSHIFTS = CLUS_IDENTIFIER[:,2]
				z_med = np.median( CLUS_REDSHIFTS )

				from astropy.cosmology import FlatLambdaCDM
				ombh2 = 0.022
				omb0 = ombh2/param_dict['h']**2.
				cosmo = FlatLambdaCDM(H0 = param_dict['h']*100., Om0 = param_dict['omega_m'], Ob0 = omb0)

				D_a = (cosmo.comoving_distance(z_med)/(1.+z_med) ).value

				new_nx, new_ny = stacked_kappa_qe.shape
				new_mapparams = [new_nx, new_ny, dx, dy]
				lx, ly = sims.get_lxly(new_mapparams)
				L = np.sqrt( lx ** 2. + ly ** 2. )
				sigma_c = D_a_sigma_c/D_a
				HOW_MUCH_SMEARING = ( fcen + fmis * np.exp( -0.5 * sigma_c**2. * L**2. ) )
				STACKED_KAPPA_MODEL_MISCENTERED_FFT = HOW_MUCH_SMEARING * MODEL_FFT
				STACKED_KAPPA_MODEL_MISCENTERED = np.fft.ifft2( STACKED_KAPPA_MODEL_MISCENTERED_FFT ).real

				"""
				subplot(131);imshow(STACKED_KAPPA_MODEL, vmin = -0.02, vmax = 0.2);colorbar()
				subplot(132);imshow(STACKED_KAPPA_MODEL_MISCENTERED, vmin = -0.02, vmax = 0.2);colorbar()
				subplot(133);imshow(STACKED_KAPPA_MODEL - STACKED_KAPPA_MODEL_MISCENTERED);colorbar()
				show()

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
class LightCone():
    def __init__(self, box_size=50, ra_min=0, ra_max=np.pi / 2, dec_min=0, dec_max=np.pi / 2, min_z=0, max_z=1, seed=None, hubble=71.0, omega_m=0.27, omega_d=0.73, unique=False, subcone_id=0):
        box_size = box_size
        self.box_size = box_size
        self.ra_min = ra_min * units.rad
        self.ra_max = ra_max * units.rad
        self.dec_min = dec_min * units.rad
        self.dec_max = dec_max * units.rad
        self.min_z = min_z
        self.max_z = max_z
        if seed is None:
            seed = random.randint(0, 1000000)
        self.seed = seed
        self.unique = unique

        self.cosmology = FlatLambdaCDM(hubble, omega_m, omega_d)
        self.subcone_id = subcone_id

    @property
    def min_dist(self):
        retval = self.cosmology.comoving_distance(self.min_z)
        # return redshift_to_comoving_distance(self.min_z)
        return retval

    @property
    def max_dist(self):
        return self.cosmology.comoving_distance(self.max_z)

    @property
    def dec(self):
        return self.dec_max - self.dec_min

    @property
    def ra(self):return self.ra_max - self.ra_min

    @property
    def dec_offset(self):
        return - self.dec_min

    @property
    def ra_offset(self):
        if not self.unique:
            return 0
        box_size = units.Mpc * self.box_size
        comparison = self.max_dist - self.min_dist * np.cos(self.ra)
        if (comparison <= box_size):
            return -self.ra_min

        def angle_root(x):
            phi = self.ra + x
            multiplier = (np.cos(x) - np.sin(x) / np.tan(phi))
            result = box_size - self.max_dist * multiplier

            single = 1 * units.Mpc
            return float(result / single)

        try:
            result = scipy.optimize.ridder(angle_root, 0.5 * np.pi, 0.0)
        except ValueError:
            return None

        return result

    @property
    def origin(self):
        if not self.unique:
            return 0, 0, 0

        theta = self.ra_offset
        phi = theta + self.ra

        # calculate the cone RA height and declination height
        h = self.max_dist * np.sin(phi) - self.min_dist * np.sin(theta)
        h_dec = self.max_dist * np.sin(self.dec)

        # How many will fit in the domain?
        ny = floor(self.box_size / h)

        return (
            self.min_dist * np.cos(phi),
            (self.subcone_id % ny) * h - self.min_dist * np.sin(theta),
            (self.subcone_id / ny) * h_dec
        )

    def tiles(self):
        #  Locate the first box to use. Start by finding a close corner
        #  of the lightcone.
        min_ra = self.ra_min + self.ra_offset
        min_dec = self.dec_min + self.dec_offset

        coord_tmp = SkyCoord(ra=min_ra, dec=min_dec, distance=self.min_dist)
        coord_tmp = coord_tmp.cartesian.x, coord_tmp.cartesian.y, coord_tmp.cartesian.z
        coords_tmp = (sum(x) for x in zip(self.origin, coord_tmp))

        # Make the box corner line up.
        # First is the box number?
        # Coords should always be within a single box
        first = []
        coords = []
        for coord in coords_tmp:
#.........这里部分代码省略.........

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
# loop over each Obs. plot

font = {'family' : 'TeX Gyre Heros',
	        'weight' : 'normal',
	        'size'   : 6}

pl.rc('font', **font)

pl.figure(1)
for i in range(nLfs):

	data = oiidata('comparat',zrArr['zmean'][i])
	nbb = len(data['Lum'])

	d1 = cosmo.comoving_distance(zmaxlist[i])
	d0 = cosmo.comoving_distance(zminlist[i])
	Vol_z = 4./3*np.pi*(d1.value**3 - d0.value**3)/8.0

	izr = np.where((GalArr['redshift'] < zmaxlist[i]) &
	(GalArr['redshift'] >= zminlist[i]))

	loii_i = loii[0][izr]

	print 'computing mock LF for redshift range ',zrArr['zmean'][i]
	lf0 = get_lf(loii_i+40.0,Vol=Vol_z,binsize=0.1,
	minsample=42.0)


	pl.subplot(240 + (i+1),xlim=[40,45.9],ylim=[-10,-1])

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]

#.........这里部分代码省略.........
        print 'R_comp = ', newrname
        print 'D_targ = ', newtname
        print 'R_targ = ', newrname_targ
    else:
        pass
    # cut the edges off ra, dec and z and z3
    ra1=ra[(ra<250)&(ra>120)&(dec<-10)&(dec>-70)]
    dec1=dec[(ra<250)&(ra>120)&(dec<-10)&(dec>-70)]
    z3 = zx3[(ra<250)&(ra>120)&(dec<-10)&(dec>-70)]
    zx = za[(ra<250)&(ra>120)&(dec<-10)&(dec>-70)]
    ra_complete = ra1[z3>0]
    dec_complete = dec1[z3>0]
    redshift_complete = zx[z3>0]
    ra_targeted = ra1[z3==1]
    dec_targeted = dec1[z3==1]
    redshift_targeted = zx[z3==1]
    a.close()
    
    # Write two new Aardvark file, first up: Complete. 
    print 'Writing ',newfname
    newf = h5py.File(newfname, "w")
    h5ra = newf.create_dataset('ra',data=ra_complete)
    h5ra.attrs["Definition"] = "Right Ascension in degrees [0,360]"
    h5dec = newf.create_dataset('dec',data=dec_complete)
    h5dec.attrs['Definition'] = "Declination in degrees [-90,90]" 
    h5z = newf.create_dataset('z',data=redshift_complete)
    h5z.attrs['Definition']="Redshift"
    zx2 = newf['z'][...]
    
    # Determining Weights
    y0,binEdges0=np.histogram(zx2,bins=60)
    z = 0.5*(binEdges0[1:]+binEdges0[:-1])
    dz = binEdges0[2] - binEdges0[1]
    rc = cosmo.comoving_distance(z)
    rc = rc / u.Mpc
    om0=1
    H0=100
    c = 299792.458
    dV = rc**2 * c * 1.738 * dz / ( H0 * (om0 * (1.0 + z)**3 + 1.0 - om0 )**0.5 )
    Nz = y0 / dV
    Nz[Nz==0]=3*10**(-9)
    yfit = np.poly1d(np.polyfit(z,np.log10(Nz),9))
    nis = np.power(10.0,yfit(zx2))
    if (weights==1)|(weights==3):
        print 'Writing unity weights'
        wg = np.ones(len(zx2))
        h5z = newf.create_dataset('weight',data=wg)
        h5z.attrs['Definition']="Weight"
        
    elif (weights==2)|(weights==4):
        print 'Writing redshift-dependant weights'
        wg = 1.0 / ( 1.0 + 3000.0 * nis)
        h5z = newf.create_dataset('weight',data=wg)
        h5z.attrs['Definition']="Weight"        
    else:
        print 'Not writing any weights'
    if debug==True:
        print 'Weight distribution for complete_data'
        xp = np.linspace(0,1.6,100)
        idxran = np.random.random_integers(0,high=len(wg)-1,size=2000000)
        xran = redshift_complete[idxran]
        yran = wg[idxran]
        fig = plt.figure()
        plt.plot(xran,yran,'+')
        plt.xlabel("Redshift")
        plt.ylabel("Weight per galaxy")

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
fout = open(fn_out, 'w')
fout_rand = open(fn_rand_out, 'w')


wd = fits[1].where("Z > 0.43 && Z < 0.7 && WEIGHT_SYSTOT >0")
sel = np.arange(len(wd))
#np.random.shuffle(sel)
data = fits[1][wd]##[sel[:min(ndata, len(wd))]]

wr = fits_rand[1].where(" Z > 0.43 && Z < 0.7")
sel = np.arange(len(wr))
np.random.shuffle(sel)
data_R = fits_rand[1][wr]##[sel[:min(ndataR, len(wr))]]

xyz = cor.ra_dec_to_xyz(data['RA'], data['DEC'])* cosmo.comoving_distance(data['Z']).value
xyz_rand = cor.ra_dec_to_xyz(data_R['RA'], data_R['DEC'])* cosmo.comoving_distance(data_R['Z']).value

for i in xrange(len(wd)):
    fout.write("%.10f %.10f %.10f \n"%(xyz[0,i], xyz[1,i], xyz[2,i]))

for i in xrange(len(wr)):
    fout_rand.write("%.10f %.10f %.10f \n"%(xyz_rand[0,i], xyz_rand[1,i], xyz_rand[2,i]))

fout.flush()
fout_rand.flush()

##============================================================##
##                        Results                             ##
##============================================================##

# 需要导入模块: from astropy.cosmology import FlatLambdaCDM [as 别名]
# 或者: from astropy.cosmology.FlatLambdaCDM import comoving_distance [as 别名]
def mcone_filam(ralim,declim,zlim,npts=80000,nvoids=2000,nstep=100,rmin=None,
                b=150.,repmethod='rotation',cosmo=None,oformat='table'):
    '''
    Generate a mock light cone of random 3D points emulating filamentary structure
    
    To create filaments, the algorithm evolves a set of random points over a 
    sphere, pushing them away in small steps from a random set of void centers.
    In each step, the coordinates are contracted a bit towards the origin. Then,
    the sphere is chopped into a cube.
    
    This cube is then stacked in 3D along the desired observation cone with 
    as many cubes as needed. This can generate artificial strucure so each cube
    can be : 

    (1) Copied as is (``repmethod='copy'``). This will likely introduce artificial patterns.
    
    (2) Randomly rotated and mirrored (``repmethod='rotation'``)
    
    (3) Generated randomly, so each box is uniquely random (``repmethod='fullrandom'``). Naturally this is slower
            
    Parameters
    ----------
    ralim : list (ra0, ra1)
        RA limits [deg]
    declim : list (dec0, dec1)
        DEC limits [deg]
    zlim : list (z0, z1)
        Redshift limits
    npts : integer
        Number of points in the initial unit sphere. Default 20000
    nvoids : integer
        Number of voids in the initial unit sphere. Default 1000
    nsteps : integer
        Number of steps that points move away from void. In general, less steps
        means softer and weaker filaments. Default 100
    rmin : float
        After box is created, remove 1 member of each pair closer than rmin [Mpc]
    b : float
        Box size in Mpc. Default 300
    repmethod : string
        Repetition method for the boxes along the cone volume
        * 'copy' : the same box is repeated
        * 'rotation' : the same box is repeated but rotated and mirrored. Default
        * 'fullrandom' : generate unique random boxes each time
    cosmo : astropy cosmology object
        If not given, defaults to a flat LambdaCDM with H_0=100 and Om=0.3
    oformat : string
        Select output format as
        * 'table' : astropy table. Default
        * 'array' : numpy array
        
    Notes
    -----
    The box is always unit size, which is then scaled to any meaningfull size
    by ``b``. The units depend on your intepretation. A choice of ``nvoids=2000`` 
    and ``b=150`` gets a box of 150 Mpc with voids and filaments of size similar 
    to real ones. In general :
    
    * For a fixed ``b``, the larger ``nsteps`` the stronger the features
    
    * For a fixed ``b``, The larger ``nvoids`` the smaller the structure
    
    * Smaller structure means more points can be included to trace small scales
      without getting unrealistic features.
    
    It is hard to estimate a priori how many points will fall inside the
    observation cone. Just give a try and change ``npts`` accordingly.

    Returns
    ----------
    kone : astropy.table / array of shape (nobj_in_cone,7)
        Table or array of 7 columns [ra,dec,redshift,comdis(redshift),x,y,z]
    '''
    if cosmo is None: cosmo = FlatLambdaCDM(H0=100, Om0=0.3)

    # Get (ra/dec/z) into units of radians and comov_Mpc  ---------------------
    #torad = np.pi/180
    #ra0, ra1  = ralim[0]*torad, ralim[1]*torad
    #dec0,dec1 = declim[0]*torad, declim[1]*torad
    d0        = cosmo.comoving_distance(zlim[0]).value
    d1        = cosmo.comoving_distance(zlim[1]).value
        
    # Build unit cube with filaments, scaled by b  ----------------------------
    xyz = filament_box(npts=npts,nvoids=nvoids,nstep=nstep,rmin=rmin,b=b)
  
    # Fill observation cone with cubes  ---------------------------------------
    xyza = fill_cone(ralim,declim,zlim,xyz=xyz,b=b,repmethod=repmethod,
                     npts=npts,nvoids=nvoids,nstep=nstep,rmin=rmin,cosmo=cosmo)

    # Prune points outside the exact observation cone  ------------------------
    x,y,z = xyza[:,0],xyza[:,1],xyza[:,2]
    ttt   = cartesian_to_spherical(x,y,z)
    comd  = ttt[0].value
    dec   = ttt[1].value*180./np.pi
    ra    = ttt[2].value*180./np.pi
    
    idx, = np.where( (ra>ralim[0]) & (ra<ralim[1]) & (dec>declim[0]) & (dec<declim[1]) & (comd>d0) & (comd<d1))
    nobj = len(idx)
    print 'Results'
    print '  |-- Prunned objects outside intersecting cubes :',len(ra)-nobj
#.........这里部分代码省略.........