Code for Random Catalog with Mask

From Wolfram Math World:

To pick a random point on the surface of a sphere, it is incorrect to select spherical coordinates θ and φ from uniform distributions θ ∈ [0,2π) and φ ∈ [0,π], since the area element dΩ=sinφdθdφ is a function of φ, and hence points picked in this way will be "bunched" near the poles (left figure above).

To obtain points such that any small area on the sphere is expected to contain the same number of points (right figure above), choose u and v to be random variates on (0,1). Then

θ = 2πu

φ = cos^-1(2v - 1)

gives the spherical coordinates for a set of points which are uniformly distributed over the surface of the sphere. This works since the differential element of solid angle is given by

dΩ = sinφdθdφ = -dθd(cosφ)

So for us to pick a random ra and dec on the sphere we want to use the following in idl:

theta = 2*!pi*RANDOMU(S, 1)
phi = acos(2*RANDOMU(S, 1) - 1)

This should give us a random theta and phi where θ ∈ [0,2π) and φ ∈ [0,π].

Unfortunately, just to be difficult HEALPix expects the following as inputs for ang2pix_nest:

theta - angle (along meridian), in [0,Pi], theta=0 : north pole
phi - angle (along parallel), in [0,2*Pi]

So their convention is the opposite of mine above so we should use the following:
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

Here is the code:

file = "maskfile.fits"
; nside_in is the number of pixels n the mask
read_fits_map, file, weight, nside=nside_in, ordering=order

; make random phi and theta
phi = 2*!pi*RANDOMU(S, 1)
theta = acos(2*RANDOMU(S, 1) - 1)

; Find pixel ipnest
ang2pix_nest, nside_in, theta, phi, ipnest

; Find weight of pixel
thisweight = weight[ipnest]

; check to make sure pixel center (thistheta, thisphi) is close to phi and theta
pix2ang_nest, nside_in, ipnest, thistheta, thisphi

(Checked this for several theta/phi and it looked good)


So here is the code to make the random catalog of size randomsize:

randomsize = 100
randomtheta = findgen(randomsize)*0.0
randomphi = findgen(randomsize)*0
index = 0

while index lt randomsize do begin
    thisindex = index 
    while index eq thisindex do begin
       phi = 2*!pi*RANDOMU(S, 1)
       theta = acos(2*RANDOMU(S, 1) - 1)
       ang2pix_nest, nside_in, theta, phi, ipnest
       thisweight = weight[ipnest]
       if thisweight gt 2.0 then begin
          randomtheta[index] = theta
          randomphi[index] = phi
          index = index+1
       endif  
    endwhile
endwhile