Python Planck13.kpc_comoving_per_arcmin方法代碼示例

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import kpc_comoving_per_arcmin [as 別名]
def get_densities(slice, large_slice):
    """Measures the density for each galaxy in the slice
    using a larger slice, because of borders effects from z
    (RA Dec not taken in accounts) and also to include all masses
    and not only the current mass bin"""

    radius = .5*u.Mpc

    densities = np.array([])
    meanmasses = np.array([])
    medianmasses = np.array([])

    print datetime.datetime.now()

    for galaxy in slice:
        z = galaxy["zb_1"]
        RA = galaxy["RA"]
        Dec = galaxy["Dec"]

        Mpc_per_deg = cosmo.kpc_comoving_per_arcmin(z).to(u.Mpc/u.deg)
        dRA = radius / Mpc_per_deg
        dDec = dRA

        slicetmp = large_slice

        slicetmp = slice_RA(slicetmp, RA, dRA)
        slicetmp = slice_Dec(slicetmp, Dec, dDec)
        #slicetmp = slice_z_cosmo(slicetmp, z, radius)
        velocity_limit = 500 #km/s

        #print zobs(-(velocity_limit), z) - zobs(velocity_limit, z)

        slicetmp = slice_z_minmax(slicetmp, zobs(-(velocity_limit), z), zobs(velocity_limit, z))

        meanmass_comp = np.mean(slicetmp["Stell_Mass_1"][np.where(slicetmp["ID"] != galaxy["ID"])])
        medianmass_comp = np.median(slicetmp["Stell_Mass_1"][np.where(slicetmp["ID"] != galaxy["ID"])])

        densities = np.append(densities, len(slicetmp)-1)
        meanmasses = np.append(meanmasses, meanmass_comp)
        medianmasses = np.append(medianmasses, medianmass_comp)

        #print len(slicetmp)-1

    slice.add_column(Column(data=densities, name="densities"))
    slice.add_column(Column(data=meanmasses, name="meanmasses"))
    slice.add_column(Column(data=medianmasses, name="medianmasses"))

    print datetime.datetime.now()

    return slice

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import kpc_comoving_per_arcmin [as 別名]
from astropy.constants import c as speed_of_light



# prefix, target = 'M2FS16', 'A04'
# outfile = os.path.abspath(os.path.join(os.curdir,'..','..', 'OneDrive - umich.edu','Research','M2FSReductions','catalogs','merged_target_lists', 'mtlz_{}_{}_full.fits'.format(prefix,target)))
prefix, target = 'M2FS17', 'B09'
outfile = os.path.abspath(os.path.join(os.curdir,'..','data','catalogs','merged_target_lists', 'mtlz_{}_{}_full.fits'.format(prefix,target)))


complete_table = Table.read(outfile,format='fits')

ra_clust = float(complete_table.meta['RA_TARG'])
dec_clust= float(complete_table.meta['DEC_TARG'])
z_clust  = float(complete_table.meta['Z_TARG'])
kpc_p_amin = Planck13.kpc_comoving_per_arcmin(z_clust)
cluster = SkyCoord(ra=ra_clust*u.deg,dec=dec_clust*u.deg)
scalar_c = consts.c.to(u.km/u.s).value

measured_table = complete_table[np.bitwise_not(complete_table['SDSS_only'])]

corcut = 0.3

spec_overlap = measured_table[np.logical_not(measured_table['sdss_zsp'].mask)]
spec_overlap = spec_overlap[((spec_overlap['cor']>corcut))]
if len(spec_overlap)==0:
    raise(TypeError,"There were no overlapping spectra and SDSS")
blue_cam = np.array([fib[0]=='b' for fib in spec_overlap['FIBNAME']])
red_cam = np.bitwise_not(blue_cam)

plt.figure()

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import kpc_comoving_per_arcmin [as 別名]
def grab_sdss_spectra(radec, radius=0.1*u.deg, outfil='tmp.fits', debug=False, maxsep=None):
    """ Grab SDSS spectra

    radius: float (0.1)
      Search radius in deg -- Astroquery actually makes a box, not a circle
    outfil: str ('tmp.fits')
      Name of output file for FITS table
    maxsep: float (None) :: Mpc
      Maximum separation to include 

    JXP on 01 Jan 2015
    """
    from astroquery.sdss import SDSS
    from astropy import coordinates as coords
    from astropy.cosmology import Planck13 as cosmo 
    from astropy.coordinates import SkyCoord



    cC = coords.SkyCoord(ra=radec[0]*u.degree, dec=radec[1]*u.degree)

    # Query
    photoobj_fs = ['ra', 'dec', 'objid', 'run', 'rerun', 'camcol', 'field']
    mags = ['petroMag_u', 'petroMag_g', 'petroMag_r', 'petroMag_i', 'petroMag_z'] 
    magsErr = ['petroMagErr_u', 'petroMagErr_g', 'petroMagErr_r', 'petroMagErr_i', 'petroMagErr_z']

    phot_catalog = SDSS.query_region(cC,spectro=True,radius=radius,
                                     photoobj_fields=photoobj_fs+mags+magsErr) # Unique
    spec_catalog = SDSS.query_region(cC,spectro=True, radius=radius) # Duplicates exist
    nobj = len(phot_catalog)

    #
    print('grab_sdss_spectra: Found {:d} sources in the search box.'.format(nobj))

    # Coordinates
    cgal = SkyCoord(ra=phot_catalog['ra']*u.degree, dec=phot_catalog['dec']*u.degree)
    sgal = SkyCoord(ra=spec_catalog['ra']*u.degree, dec=spec_catalog['dec']*u.degree)
    sepgal = cgal.separation(cC) #in degrees

    # Check for problems and parse z
    zobj = np.zeros(nobj)
    idx, d2d, d3d = coords.match_coordinates_sky(cgal, sgal, nthneighbor=1)
    if np.max(d2d)*3600. > 1.*u.Unit('arcsec'):
        print('No spectral match!')
        xdb.set_trace()
    else:
        zobj = spec_catalog['z'][idx]

    idx, d2d, d3d = coords.match_coordinates_sky(cgal, cgal, nthneighbor=2)
    if np.min(d2d.to('arcsec')) < 1.*u.Unit('arcsec'):
        print('Two photometric sources with same RA/DEC')
        xdb.set_trace()

    #xdb.set_trace()


    # Cut on Separation
    if not maxsep is None:
        print('grab_sdss_spectra: Restricting to {:g} Mpc separation.'.format(maxsep))
        sepgal_kpc = cosmo.kpc_comoving_per_arcmin(zobj) * sepgal.to('arcmin')
        sepgal_mpc = sepgal_kpc.to('Mpc')
        gdg = np.where( sepgal_mpc < (maxsep * u.Unit('Mpc')))[0]
        phot_catalog = phot_catalog[gdg]
        #xdb.set_trace()

    nobj = len(phot_catalog)
    print('grab_sdss_spectra: Grabbing data for {:d} sources.'.format(nobj))

    # Grab Spectra from SDSS

    # Generate output table
    attribs = galaxy_attrib()
    npix = 5000 #len( spec_hdus[0][1].data.flux )
    spec_attrib = [(str('FLUX'), np.float32, (npix,)),
                   (str('SIG'), np.float32, (npix,)),
                   (str('WAVE'), np.float64, (npix,))]
    tbl = np.recarray( (nobj,), dtype=attribs+spec_attrib)

    tbl['RA'] = phot_catalog['ra']
    tbl['DEC'] = phot_catalog['dec']
    tbl['TELESCOPE'] = str('SDSS 2.5-M')

    # Deal with spectra separately (for now)
    npix = 5000 #len( spec_hdus[0][1].data.flux )
    
    for idx,obj in enumerate(phot_catalog): 
        #print('idx = {:d}'.format(idx))

        # Grab spectra (there may be duplicates)
        mt = np.where( sgal.separation(cgal[idx]).to('arcsec') < 1.*u.Unit('arcsec'))[0] 
        if len(mt) > 1:
            # Use BOSS if you have it
            mmt = np.where( spec_catalog[mt]['instrument'] == 'BOSS')[0]
            if len(mmt) > 0:
                mt = mt[mmt[0]]
            else:
                mt = mt[0]
        elif len(mt) == 0:
            xdb.set_trace()
        else:
#.........這裏部分代碼省略.........