Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
exit(1)
fname=sys.argv[1]
fmt=sys.argv[2]
if fmt=='ASCII' :
ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_ascii(fname)
elif fmt=='FITS' :
ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_fits(fname)
elif fmt=='HDF5' :
ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_hdf5(fname,1)
nside=64
npix=hp.nside2npix(nside)
mp=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix])[0]
mpr=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=rsd_arr)[0]
mpe1=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=e1_arr)[0]
mpe2=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=e2_arr)[0]
plt.figure();
plt.hist(z_arr,bins=100,histtype='step');
plt.xlabel('$z$',fontsize=16); plt.ylabel('$N(z)$',fontsize=16);
hp.mollview(mp,title='$N_g(\\hat{\\bf n})$')
hp.mollview(mpr/mp,title='$v_r$')
hp.mollview(mpe1/mp,title='$e_1$')
hp.mollview(mpe2/mp,title='$e_2$')
if dskw_arr is not None :
#Plot skewers
nside_skw=64
ipix_skw=1000
id_in_pix=np.where(hp.ang2pix(nside_skw,(90-dec_arr)*np.pi/180,ra_arr*np.pi/180)==ipix_skw)[0]
cols=['r','g','b','y','k']
plt.figure(); plt.title('Skewers'); plt.xlabel('$z$'); plt.ylabel('$\\delta$');
self.richness_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
self.richness_lower_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
self.richness_upper_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
self.richness_upper_limit_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
self.stellar_mass_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
self.fraction_observable_sparse_array = numpy.zeros([len_distance_modulus, len_pixels_target])
## Calculate the average stellar mass per star in the ischrone once
#stellar_mass_conversion = self.isochrone.stellarMass()
# Specific pixel
if coords is not None and distance_modulus_index is not None:
lon, lat = coords
theta = numpy.radians(90. - lat)
phi = numpy.radians(lon)
pix_coords = healpy.ang2pix(self.config.params['coords']['nside_pixel'], theta, phi)
lon, lat = self.roi.centers_lon_target, self.roi.centers_lat_target
logger.info('Begin loop over distance moduli in likelihood fitting ...')
for ii, distance_modulus in enumerate(self.distance_modulus_array):
# Specific pixel
if distance_modulus_index is not None and ii!=distance_modulus_index:
continue
logger.info(' (%i/%i) distance modulus = %.2f ...'%(ii+1, len_distance_modulus, distance_modulus))
self.u_color = self.u_color_array[ii]
self.observable_fraction_sparse = self.observable_fraction_sparse_array[ii]
for jj in range(0, len_pixels_target):
# Specific pixel
def view_observed_gsm(self, logged=False, show=False, **kwargs):
""" View the GSM (Mollweide), with below-horizon area masked. """
sky = self.observed_sky
if logged:
sky = np.log2(sky)
# Get RA and DEC of zenith
ra_rad, dec_rad = self.radec_of(0, np.pi / 2)
ra_deg = ra_rad / np.pi * 180
dec_deg = dec_rad / np.pi * 180
# Apply rotation
derotate = hp.Rotator(rot=[ra_deg, dec_deg])
g0, g1 = derotate(self._theta, self._phi)
pix0 = hp.ang2pix(self._n_side, g0, g1)
sky = sky[pix0]
coordrotate = hp.Rotator(coord=['C', 'G'], inv=True)
g0, g1 = coordrotate(self._theta, self._phi)
pix0 = hp.ang2pix(self._n_side, g0, g1)
sky = sky[pix0]
hp.mollview(sky, coord='G', **kwargs)
if show:
if not plt:
import pylab as plt
plt.show()
return sky
def __call__(self):
theta, phi = eq2rad(self.ctx.coords.ra, self.ctx.coords.dec)
pix_numbers = hp.ang2pix(self.ctx.params.nside, theta, phi)
cntr = Counter(pix_numbers)
tod = self.ctx.tod_vx
path = tempfile.mktemp()
with h5py.File(path, "w") as fp:
restructure(fp, tod, cntr, pix_numbers)
self.ctx.restructured_tod_path = path
self.ctx.restructured_tod_pixels = cntr.keys()
del self.ctx.tod_vx
del self.ctx.tod_vy
del self.ctx.ref_channel
del self.ctx.coords
def distance_from_healpix(nside, points, omap=None, odomains=None, domains=False, rmax=None, method="bubble"):
"""Find the distance from each pixel in healpix map with nside nside to the
nearest of the points[{dec,ra},npoint], returning a [ny,nx] map of distances.
If domains==True, then it will also return a [ny,nx] map of the index of the point
that was closest to each pixel. If rmax is specified, then distances will only be
computed up to rmax. Beyond that distance will be set to rmax and domains to -1.
This can be used to speed up the calculation when one only cares about nearby areas."""
import healpy
from pixell import distances
info = distances.healpix_info(nside)
if omap is None: omap = np.empty(info.npix)
if domains and odomains is None: odomains = np.empty(info.npix, np.int32)
pixs = utils.nint(healpy.ang2pix(nside, np.pi/2-points[0], points[1]))
return distances.distance_from_points_healpix(info, points, pixs, rmax=rmax, omap=omap, odomains=odomains, domains=domains, method=method)
This is used in case one wants to model the baseline only based on the
low Galactic emission regions. There is one mask template under the
data directory that can be used.
:param nside: Nside of healpix mask
:param mask_original: the original TOD mask
:param mask_gal: the healpix mask
:param ra: RA coordinates corresponding to the TOD
:param dec: DEC coordinates corresponding to the TOD
:return: final TOD mask
"""
theta = (np.pi / 2 - dec)
phi = ra
pix = hp.ang2pix(nside, theta, phi, nest=False)
new_mask = np.zeros_like(mask_original)
new_mask[:, mask_gal[pix]] = True
return new_mask
def superpixel(subpix, nside_subpix, nside_superpix):
"""
Return the indices of the super-pixels which contain each of the sub-pixels.
"""
if nside_subpix==nside_superpix: return subpix
theta, phi = hp.pix2ang(nside_subpix, subpix)
return hp.ang2pix(nside_superpix, theta, phi)
star = Star(hip=hip, hd=hd, vmag=vmag,
ra=ra, de=de, plx=plx, pra=0, pde=0, bv=bv)
stars[hd] = star
stars = sorted(stars.values(), key=lambda x: (x.vmag, x.hd))
# Create the tiles.
# For each star, try to add to the lowest tile that is not already full.
print 'create tiles'
max_vmag_order0 = 0
tiles = {}
for s in stars:
for order in range(8):
pix = healpy.ang2pix(1 << order, pi / 2.0 - s.de, s.ra, nest=True)
nuniq = pix + 4 * (1L << (2 * order))
assert int(log(nuniq / 4, 2) / 2) == order
tile = tiles.setdefault(nuniq, [])
if len(tile) >= MAX_SOURCES_PER_TILE: continue # Try higher order
tile.append(s)
if order == 0: max_vmag_order0 = max(max_vmag_order0, s.vmag)
break
print 'save tiles'
out_dir = './data/stars/'
COLUMNS = [
{'id': 'hip', 'type': 'i'},
{'id': 'hd', 'type': 'i'},
{'id': 'vmag', 'type': 'f', 'unit': eph.UNIT_VMAG, 'zerobits': 16},
{'id': 'ra', 'type': 'f', 'unit': eph.UNIT_RAD, 'zerobits': 8},
{'id': 'de', 'type': 'f', 'unit': eph.UNIT_RAD, 'zerobits': 8},
def get_map(ctx):
"""
Function for creating maps from TOD by a simple averaging in healpy
pixels.
:param ctx: context that contains TOD, frequencies, coordinate
information, and an RFI mask
"""
npix = hp.nside2npix(ctx.params.nside)
nfreq = ctx.frequencies.shape[0]
theta, phi = eq2rad(ctx.coords.ra, ctx.coords.dec)
inds = hp.ang2pix(ctx.params.nside, theta, phi)
maps = np.zeros((nfreq, 2, npix))
counts = np.zeros(maps.shape)
varflag = ctx.params.variance
varmaps = np.zeros((nfreq, 2, npix)) if varflag else None
for tod_ind in range(min(theta.shape[0], ctx.tod_vx.shape[1])):
ind = inds[tod_ind]
updateMaps(ctx.tod_vx, maps, counts, tod_ind, ind, 0, varmaps)
updateMaps(ctx.tod_vy, maps, counts, tod_ind, ind, 1, varmaps)
redshifts = 1420.40575177 / ctx.frequencies - 1
return_maps = varmaps if varflag else maps * counts
return return_maps, redshifts, counts
def AzimuthalRotation(hmap):
"""
Azimuthal rotation of a healpix map by pi/2 about the z-axis
"""
npix = len(hmap)
nside= hp.npix2nside(npix)
hpxidx = np.arange(npix)
t2,p2 = hp.pix2ang(nside, hpxidx)
p = p2 - np.pi/2
p[p < 0] += 2. * np.pi
t = t2
idx = hp.ang2pix(nside, t, p)
hout = hmap[idx]
return hout