Monday, October 25, 2010

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


  1. Computationally, there is a faster way. It turns out that the points are equally distributed on the z-axis, such that
    z = [-R,R],
    t = [0, 2*pi],
    r = sqrt(R*R - z*z),
    x = r*cos(t),
    y = r*sin(t)
    (from Knuth vol.2)

    It's based on the following proof that the surface area of a segment is equal on the height division:

    Radius of circle plane at z, r = sqrt(R*R-z*z)
    Differential radius, f'(z) = -z/sqrt(R*R-z*z)
    f'(z)^2 = z*z / (R*R-z*z)
    1+f'(z)^2 = R*R / (R*R-z*z)

    Surface of revolution between a + b:
    S(a,b) = 2*pi*int[b,a](f(z)*sqrt(1+f'(z)*f'(z)))dz
    S(a,b) = 2*pi*int[b,a]()dz
    S(a,b) = 2*pi*(b-a)

    As the result is b-a, the surface area only depends on where the the difference in height, rather than where in the sphere it is cut.

    As a result, you can reduce this computationally to 2 RND, 2 trig + 1 sqrt, rather than 2 RND, 6 trig.