Reading Healpix and FITS files

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Jamie Portsmouth

Reading Healpix and FITS files

Post by Jamie Portsmouth » October 10 2004

I'm unfamiliar with dealing with CMB maps.
I would like to be able to:

a) Take say one of the publicly available WMAP map files, and
extract the sky temperature as a function of pixel (and I'd really like
the pixel coordinates too..)

b) Take the FITS format output of say the Healpix routine anafast
and read the C_l's and alm's. I realize there are f90 routines
provided with Healpix which do this, but i'm using OS X and was
never able to get the relevant f90 routines to compile with all the
dependencies in the makefile.

What is the best way to go about doing this (short of being a FITS guru master)?

Anze Slosar
Posts: 183
Joined: September 24 2004
Affiliation: Brookhaven National Laboratory
Contact:

Post by Anze Slosar » October 11 2004

Try something like this:

Code: Select all

SUBROUTINE LOAD_HEAL_FILE(filename, DATAR, nside, npixtot)
use healpix_types
use fitstools
use pix_tools
use alm_tools
use utilities 

implicit none
CHARACTER(LEN=*), PARAMETER :: code = 'tlike'
CHARACTER (LEN=*) filename
INTEGER, INTENT (OUT), OPTIONAL :: nside, npixtot
INTEGER(I4B) :: npixtoto, nmapso, info, status
INTEGER(I4B) :: ordering, obs_npix, nsideo, mlpol, type, polarisation

REAL(SP), DIMENSION(:,:), POINTER :: DATAR


     IF (debug.gt.0) PRINT *,'Reading file ', filename
     npixtoto = getsize_fits(filename, nmapso, ordering, obs_npix, nsideo, mlpol, type, polarisation )
     IF (debug.gt.1)   PRINT *,'nsideo=', nsideo
     ALLOCATE(DATAR(0:12*nsideo**2-1,1:2), stat = status)
     IF (status /= 0) CALL die_alloc(code,'alloc '//filename)
     CALL input_map(filename, DATAR, npixtoto, nmapso, fmissval=-1.6375e30_sp)

     IF (ordering == 2 ) then 
        IF (debug.gt.0) PRINT *,'Converting NEST to RING...'
        CALL convert_nest2ring(nsideo,DATAR(:,1))
        CALL convert_nest2ring(nsideo,DATAR(:,2))
        ordering =1 
     ENDIF
     if (present(npixtot)) npixtot=npixtoto
     if (present(nside)) nside = nsideo
end
with

Code: Select all

   REAL(SP), DIMENSION(:,:), POINTER :: map_t
and to get positions you should do something like:

Code: Select all

   DO II=0,NPIXTOT-1
         CALL PIX2VEC_RING(nside, ii, VEC1)
           !! Now VEC1 contains point on a sphere oorresponding to the
           !!   current pixel

Ben Gold
Posts: 81
Joined: September 25 2004
Affiliation: University of Minnesota
Contact:

Post by Ben Gold » October 12 2004

If you're desperate for a FORTRAN 95 compiler for MacOS X, you can get a binary for one here:

http://www.g95.org/

I haven't used it much but CAMB at least compiles and seems to run OK. If you want to wait for a less experimental compiler, gcc 4.0.0 is scheduled to include FORTRAN 95; it's tentatively scheduled for release early next year.

Or am I misinterpreting your compiler difficulties?

Post Reply