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

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)

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 =

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:

So their convention is the opposite of mine above so we should use the following:

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

theta =

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:

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

Computationally, there is a faster way. It turns out that the points are equally distributed on the z-axis, such that

ReplyDeletez = [-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.

Thank you borandi!

ReplyDelete