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

Davis to the Rescue

I was at the cosmology seminar today and I was complaining to Nic Ross about my correlation function woos and Marc Davis happened to be in earshot and told me that I was thinking about this all wrong. Basically what I should do is take my galaxies with coordinates of Ra and Dec and calculate their x,y,z coordinates (projecting them onto a unit sphere). Then I can calculate the 3D correlation function on these objects (which I already have code to do), where I can translate between their physical separation in x, y, z coordinates to their angular separation by the following:

Ra₁, Dec₁ → x₁, y₁, z₁
Ra₂, Dec₂ → x₂, y₂, z₂

Separation of objects (s) in Cartesian coordinates:

s = [(x₁ - x₂)² + (y₁ - y₂)² + (z₁ - z₂)²]^½

Geometry relates the separation angle (γ) to Cartesian separation (s):

s = 2sin(γ/2)

or using the dot product of the two vectors:

x₁x₂ + y₁y₂ + z₁z₂ = cos( γ)

Therefore, the only change I need to make is after I grid objects, and select which objects are 'close by' using above conversion from s to γ, is to convert the 3D separation of the objects to the angular separation. This shouldn't be very difficult to do using the 3D code I already have.

Thank you Marc Davis for making my life a lot better today!

Useful links for writing equations in HTML (used above):
http://www.w3schools.com/tags/ref_entities.asp
http://www.tizag.com/htmlT/htmlsuperscript.php
http://htmlhelp.com/reference/html40/entities/symbols.html

Cartesian to Spherical

I've spent most of the day pouring over the 3D correlation function code trying to figure out the best way to convert it to an angular correlation function. Here is how the 3D correlation function works and my plan to change it:

There are two galaxy files which we are going to correlate, file S and file P. The angular region of the galaxies in each file is the same, and for simplicity sake let's assume it is the whole sky (in actuality it is some mask of Ra and Dec region(s)). The sky is divided up into a grid using lines that are evenly spaced in RA and Dec. Number of grid lines is set by the user.

All the galaxies in the S file are split up into the grid, such that for every grid section there is a list (or chain) of the galaxies contained in that grid.

A galaxy (X) is chosen from the P file and then placed into the S file grid. The distances between X and the S galaxies in that same grid section, as well as the S galaxies in close by grid sections are calculated. 'Close by' sections is defined as grids which are within the maximum correlation distance to the grid. So for instance if we are correlating out to 20 Mpc/h and the grids are spaced at 5 Mpc/h, then you could go out 4 grid sections in each direction, and thus correlate in a box which is 8X8 grids (in two dimensions). This prevents you from calculating correlations between galaxies that are further away than your maximum correlation distance.

Finding the 'close by' grids is very easy in Cartesian coordinates (as I described above), but very difficult in spherical coordinates. For example, if you have two objects on the north pole, separated by a small Dec, then it essentially doesn't matter what their RA separation is, they will still have a small angular separation. Whereas on the equator if objects have a small Dec deparation, but a large RA separation, then this means they are on different sides of the sky. So you cant simply look at grids with angle separation of +/- the maximum angular correlation because this breaks down as you go to the poles.

What I suggest to do is to create a distance matrix which has the angular distance (γ) from the center of each grid point to the center of every other grid point using this formula:

cos γ = cos(90 - Dec1)cos(90 - Dec2) + sin(90 - Dec1)sin(90 - Dec2)cos(RA1 - RA2)

This matrix only needs to be calculated once per grid, and can be saved to a file, and then read in if the same grid is used again.

Then I follow the same correlation strategy as above expect when finding the 'close by' grid sections, I go to this matrix and find all the point which are within the maximum correlation angle to the current grid and then calculate distances to all galaxies within these grids. I can speed things up further, but inputting objects from the P file grid by grid, so that once I find all the 'close by' S galaxies to galaxy X, I can then calculate distance to all other P objects in the same grid as X with the same S galaxies.

The main increase in time here will be due to making the matrix and then scanning the matrix to figure out which objects are closest, but I did a back of the envelope calculation, and I think the processing time used to do this is many orders of magnitude less than the processing time needed to correlate the galaxies with each other (10^8 vs 10^12) so hopefully this extra processing time will be in the noise.

Everything Working (finally)

Princeton (and Alexia) makes everything better! It seems I have finally got a working 3D correlation function. I don't know why it took me so long to get this thing to work, it seems like it should be simple enough to do. Anyway, here is a summary what I've done...

The inputs to the function are a set of mock "data" point in Cartesian coordinates. For Sloan data, these will be converted from ra/dec (spherical coordinates) to Cartesian in python. There are also mask inputs, both in spherical and Cartesian coordinates.

For the mock data I simply applied a mask cut on ra/dec/redshift and then converted those points back to x/y/z:

The mask in ra/dec/redshift is a contiguous box.

Converted to x/y/z

The mask has the minimum and maximum values that the data can fall in for each coordinate. For example in Cartesian coordinates the mask contains:

In [541]: minX
Out[541]: 173.568011004

In [542]: maxX
Out[542]: 449.984440618

In [543]: minY
Out[543]: -289.251614958

In [544]: maxY
Out[544]: 289.251614958

In [545]: minZ
Out[545]: -190.178217783

In [546]: maxZ
Out[546]: 190.178217783

The data is then scaled down to a 1x1x1 box (this is what the Alexia/Martin correlation function calculation code is expecting). This is done by doing the following to each dimension of each data point (i):

posX[i] = (dataX[i] - minX + padding*rmax/2)/maxBoxside

where dataX[i] is the x value of the ith data point, minX is the minimum value that x can be (from the mask), padding is how much padding you want around the edge of your data (this prevents power from being "wrapped" around as the correlation calculation uses periodic boundary conditions), rmax is the maximum distance you are calculating the correlation function out to, and maxBoxside is the length of the longest side of the databox [i.e. max(maxX - minX, maxY - minY, maxZ - minZ)].

Because the data falls within a contiguous ra/dec/redshift, I populate the randoms in the ra/dec/redshift mask and then convert them x/y/z using the same conversion as I do on the data. I then apply the same scaling (as described in previous paragraph) to the x/y/z randoms. The result is data which falls on top of randoms and is contained in a padded 1x1x1 box:

As you can see all the data falls between 0 and 1 and the data falls on top of the randoms.

It is hard to see from the above plots but I would also like the redshift distribution of the randoms to follow that of the data. This is done by binning the data into redshift bins (20 in the example I am plotting here) and then for every data point in a particular bin, I generate 10 random points in the same bin.

Histogram of number of point in each redshift bin.
I multiplied the data by 10 so that the scale is the same as the randoms.

Once I have both the data and the randoms in a padded 1x1x1 box (making sure the randoms follow the same redshift distribution as the data), then the 3D correlation function can be calculated. This calculation is done in Cartesian coordinates, but this shouldn't matter because we are just looking at distances of points from each other and so as long as each side of the box is in the same units we are good.

Here is a comparison of my 3D correlation function with Alexia's working 3D function. The reason they don't fall exactly on top of each other is because Alexia's is calculated on different data points (in the same mock catalog) due to her's requiring a Cartesian mask for the data:

My working 3D correlation function!

Now I get to run it on the Sloan data and see if the reconstruction still fails or if that fixes my problem. If it works, I am done with my PhD thesis (well not really, but it would be huge progress). Let's keep our fingers crossed!

Bad Blogger

I've been a bad blogger this past week. Apologies for not posting. I could give excuses, but that would break two of my blog rules, so I'll quit while I am ahead (or behind as it were).

Here is a summary of where I am with this darn 3D correlation function. The Sloan data is in ra/dec/redshift coordinates. However, it is easier to calculate the 3d correlation function in x-y-z comoving coordinates. So I convert the data into comoving coordinates to calculate the correlation function. However the data lives in a ra/dec/redshift space (mask) and this corresponds to a non-boxlike x-y-z space. Therefore I need to apply the mask in ra/dec/redshift for both the randoms and the data, and then convert both to x-y-z to calculate the correlation function.

But the mock data I am currently testing my 3D correlation function with is, actually in x-y-z coordinates to begin with. This is causing an issue, because I am converting it to ra/dec/redshift, which results in non-uniform data distribution in ra/dec/redshift space (because the mock data, unlike the Sloan data, is a contiguous box in x-y-z space). However the randoms are generated in uniform ra/dec/redshift space and converted to x-y-z space (to match the Sloan mask). When I compare the mock data with the randoms they don't have the same masks because of this issue:

I think what I need to do is take the mock data in x-y-z, convert it to ra/dec/redshift and then apply a mask in that coordinate system. Then convert it back to x-y-z and use those points as my data, and then apply the same mask to the randoms. This will be more similar to what I will be doing with the Sloan data and should get my data and randoms to fall in the same location on my vector space.

Coordinate Confusion

Time to run my 3D autocorrelation function on a mock catalog where we know the answer to try to see if there is a problem with my code or that the reconstruction is failing for another reason. I downloaded the following raw mock LRG catalog from Martin White's web page: halo_000_0.8000.dat.gz. I discovered that I can transfer files to riemann via the insure replacement by doing a two step scp. That is useful to know.

The code I am using to convert from ra, dec (degrees), comoving distance (r) to x, y, z and back:

%  x=r*sin(pi/2-pi/180*dec)*cos(pi/180*ra);
%  y=r*sin(pi/2-pi/180*dec)*sin(pi/180*ra);
%  z=r*cos(pi/2-pi/180*dec) ;

%  r = sqrt(x^2 + y^2 + z^2)
%  dec =  90 - 180/pi*arccos(z/r)
%  ra = 180/pi*arcsin(y/sqrt(x^2 + y^2))

I am a little confused because these mock files are in Cartesian coordinates and if I convert them using the above code into ra and dec, I don't get the objects populating all of ra/dec space:

I constrained the data to be in a sphere of radius 0.5 boxsize (936.0 Mpc/h), and I would think that the ra and dec would then go from 0 to 360 and 0 to 180 respectively.

After some digging, I discovered that by using the arctan2 function I get the proper range:

%  theta = arctan2(sqrt(x^2+y^2), z)
%  phi = arctan2(x, y)

However when I try to then convert back to Cartesian I get a strange answer when I plot the original z versus the converted back z:

This is distressing!

Sources:
http://astro.uchicago.edu/cosmus/tech/code/radecz2xxyyzz.m
http://www.atlasoftheuniverse.com/cosmodis.c
http://www.math.montana.edu/frankw/ccp/multiworld/multipleIVP/spherical/body.htm
http://www.daniweb.com/forums/post860497-2.html

Python tips for the day
1) If you are getting the following error when trying to plot:

RuntimeError: Agg rendering complexity exceeded.
Consider downsampling or decimating your data

Then exit and restart your session and try again and it for some reason works. Stupid python.

2) When converting from spherical use atan2 not asin.