How to use the healpy.ang2pix function in healpy

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