Millimeter-wave map filtering and Sunyaev-Zel'dovich galaxy cluster/source detection package. Originally developed for the Atacama Cosmology Telescope project.
I think the file nemo/examples/SOSims/readWebSkyInputCatalog.py is not synchronized with latest Nemo's commits (makeACTName not found), and also I had to copy paste the web sky/cosmology.py file to get zofchi.
Attached is the modified file that runs for me:
"""
Reads in the halos.pksc file and outputs .fits and .reg format catalogs for
the 1 million most massive halos (M200m > 1e14 MSun). Cuts down to SO dec
range.
"""
import os
import sys
import healpy as hp
import astropy.table as atpy
from nemo import catalogs
# from cosmology import *
import numpy as np
import IPython
from scipy.interpolate import *
import numpy as np
omegab = 0.049
omegac = 0.261
omegam = omegab + omegac
h = 0.68
ns = 0.965
sigma8 = 0.81
c = 3e5
H0 = 100*h
nz = 100000
z1 = 0.0
z2 = 6.0
za = np.linspace(z1,z2,nz)
dz = za[1]-za[0]
H = lambda z: H0*np.sqrt(omegam*(1+z)**3+1-omegam)
dchidz = lambda z: c/H(z)
chia = np.cumsum(dchidz(za))*dz
zofchi = interp1d(chia,za)
def makeACTName(RADeg, decDeg, prefix = 'ACT-CL'):
"""Makes ACT cluster name from RADeg, decDeg
"""
actName=prefix+" J"+makeRA(RADeg)+makeDec(decDeg)
return actName
def makeRA(myRADeg):
"""Makes RA part of ACT names.
"""
hours=(myRADeg/360)*24
if hours<10:
sHours="0"+str(hours)[0]
else:
sHours=str(hours)[:2]
mins=float(str(hours)[str(hours).index("."):])*60
if mins<10:
sMins="0"+str(mins)[:3]
else:
sMins=str(mins)[:4]
return (sHours+sMins)#[:-2] # Trims off .x as not used in ACT names
#------------------------------------------------------------------------------------------------------------
def makeDec(myDecDeg):
"""Makes dec part of ACT names
"""
# Positive
if myDecDeg>0:
if myDecDeg<10:
sDeg="0"+str(myDecDeg)[0]
else:
sDeg=str(myDecDeg)[:2]
mins=float(str(myDecDeg)[str(myDecDeg).index("."):])*60
if mins<10:
sMins="0"+str(mins)[:1]
else:
sMins=str(mins)[:2]
return "+"+sDeg+sMins
else:
if myDecDeg>-10:
sDeg="-0"+str(myDecDeg)[1]
else:
sDeg=str(myDecDeg)[:3]
mins=float(str(myDecDeg)[str(myDecDeg).index("."):])*60
if mins<10:
sMins="0"+str(mins)[:1]
else:
sMins=str(mins)[:2]
return str(sDeg+sMins)
rho = 2.775e11*omegam*h**2 # Msun/Mpc^3
# f=open('../WebSky/halos.pksc', 'rb')
f=open('/path/to/halos.pksc','rb')
N=np.fromfile(f,count=3,dtype=np.int32)[0]
# only take first five entries for testing (there are ~8e8 halos total...)
# comment the following line to read in all halos
N = 1000000
catalog=np.fromfile(f,count=N*10,dtype=np.float32)
catalog=np.reshape(catalog,(N,10))
x = catalog[:,0]; y = catalog[:,1]; z = catalog[:,2] # Mpc (comoving)
vx = catalog[:,3]; vy = catalog[:,4]; vz = catalog[:,5] # km/sec
R = catalog[:,6] # Mpc
# convert to mass, comoving distance, radial velocity, redshfit, RA and DEc
M = 4*np.pi/3.*rho*R**3 # Msun (M200m)
chi = np.sqrt(x**2+y**2+z**2) # Mpc
vrad = (x*vx + y*vy + z*vz) / chi # km/sec
redshift = zofchi(chi)
theta, phi = hp.vec2ang(np.column_stack((x,y,z))) # in radians
decDeg=-1*(np.degrees(theta)-90) # Because HEALPix
RADeg=np.degrees(phi)
names=[]
for ra, dec in zip(RADeg, decDeg):
names.append(makeACTName(ra, dec, prefix= 'MOCK-CL'))
outFileName="halos.fits"
tab=atpy.Table()
tab.add_column(atpy.Column(names, 'name'))
tab.add_column(atpy.Column(RADeg, 'RADeg'))
tab.add_column(atpy.Column(decDeg, 'decDeg'))
tab.add_column(atpy.Column(M, 'M200m'))
tab.add_column(atpy.Column(redshift, "z"))
tab=tab[np.where(tab['decDeg'] < 30)]
tab=tab[np.where(tab['decDeg'] > -70)]
tab.write(outFileName, overwrite = True)
catalogs.catalog2DS9(tab, "halos.reg")
#IPython.embed()
#sys.exit()
### e.g. project to a map, matching the websky orientations
#nside = 1024
#map = np.zeros((hp.nside2npix(nside)))
#pix = hp.vec2pix(nside, x, y, z)
#pix = hp.ang2pix(nside, theta, phi) does the same
#weight = 1. #1 for number density, array of size(x) for arbitrary
#np.add.at(map, pix, weight)
Hi,
I think the file nemo/examples/SOSims/readWebSkyInputCatalog.py is not synchronized with latest Nemo's commits (makeACTName not found), and also I had to copy paste the web sky/cosmology.py file to get zofchi. Attached is the modified file that runs for me: