Home
Manual
Packages
Global Index
Keywords
Quick Reference
|
/*
* healpix.i - Read/write A(l,m) files.
*
*-----------------------------------------------------------------------------
* $Id$
* $Log$
*/
local ALM_MAGIC ;
ALM_MAGIC = 0xDEADC0DE;
func alm_lmax (filename)
/* DOCUMENT alm_lmax(filename)
Get value of LMAX in A(l,m) file FILENAME (can be used to assert
validity of file).
SEE ALSO: alm_read, alm_write. */
{
hdr = array(long, 2);
_read, open(filename, "rb"), 0, hdr;
if (hdr(1) != ALM_MAGIC) error, "bad magic number in \""+filename+"\"";
lmax = hdr(2);
if (lmax < 0) error, "bad value of LMAX in \""+filename+"\"";
return lmax;
}
func alm_read (filename, unpack=)
/* DOCUMENT alm_read(filename)
Read data in A(l,m) file FILENAME. If keyword UNPACK is true, the
result is a 2×(LMAX+1)×(LMAX+1) array ortherwise the result is the
lower triangular part of A(l,m) packed as a 2×NUMBER array where:
NUMBER = 1 + LMAX*(LMAX + 3)/2
In any case, the result is such that:
RESULT(1,..) = real part of the A(l,m)
RESULT(2,..) = imaginary part of the A(l,m)
SEE ALSO: alm_lmax, alm_write. */
{
f = open(filename, "rb");
hdr = array(long, 2);
_read, f, 0, hdr;
if (hdr(1) != ALM_MAGIC) error, "bad magic number in \""+filename+"\"";
lmax = hdr(2);
if (lmax < 0) error, "bad value of LMAX in \""+filename+"\"";
number = 1 + lmax*(lmax + 3)/2;
data = array(float, 2, number);
_read, f, sizeof(hdr), data;
close, f;
if (unpack) return alm_unpack(data);
return double(data);
}
func alm_write (filename, a, overwrite=)
/* DOCUMENT alm_write, filename, a;
Write A(l,m) data stored in array A into file FILENAME. If keyword
OVERWRITE is true, FILENAME is overwritten if it already exists.
Input array A can be a real 2×(LMAX+1)×(LMAX+1) array or a complex
(LMAX+1)×(LMAX+1) array or only the lower triangular part of A(l,m)
packed as a complex vector with NUMBER elements or a 2×NUMBER real
array. If A is real, then A(1,..) and A(2,..) are the real and
imaginary parts of the A(l,m).
SEE ALSO: alm_lmax, alm_read, alm_pack. */
{
/* Make A into packed form. */
local lmax;
a = alm_pack(a, lmax);
/* open binary file and write data
(take care of log file and existing target) */
if (! overwrite && open(filename, "r", 1))
error, "file \""+filename+"\" already exists";
logfile = filename + "L";
if (open(logfile, "r", 1)) logfile = string(0);
f = open(filename, "wb");
if (logfile) remove, logfile;
hdr = array(long, 2);
hdr(1) = ALM_MAGIC;
hdr(2) = lmax; /* LMAX */
_write, f, 0, hdr;
_write, f, sizeof(hdr), a;
}
func alm_pack (a, &lmax)
/* DOCUMENT alm_pack(a)
-or- alm_pack(a, lmax)
Returns packed A(l,m) data from array A. Input array A can be a real
2×(LMAX+1)×(LMAX+1) array or a complex (LMAX+1)×(LMAX+1) array or only
the lower triangular part of A(l,m) packed as a complex vector with
NUMBER elements or a 2×NUMBER real array. If A is real, then A(1,..)
and A(2,..) are the real and imaginary parts of the A(l,m). The
result is an array of float's with dimension 2×NUMBER such that:
RESULT(1, ) = real part of the A(l,m)
RESULT(2, ) = imaginary part of the A(l,m)
and with:
NUMBER = 1 + LMAX*(LMAX + 3)/2.
Optional argument LMAX is an output variable to store the value of
the parameter.
SEE ALSO: alm_unpack, alm_write. */
{
/* Get significant triangular part of array (notation: N = LMAX + 1). */
dims = dimsof(a);
s = structof(a);
local data;
if (s == complex) {
/* Input is COMPLEX. */
if (dims(1) == 1) {
number = dims(2);
n = (long(sqrt(8*number + 1) + 0.5) - 1)/2;
if (number == n*(n + 1)/2) {
data = array(float, 2, number);
data(1, ) = float(a);
data(2, ) = a.im;
}
} else if (dims(1) == 2) {
if ((n = dims(2)) == dims(3)) {
a_re = float(a);
a_im = float(a.im);
data = array(float, 2, n*(n + 1)/2);
i = 1;
for (l=1 ; l<=n ; ++l) {
data(1, i:i+l-1) = a_re(l, 1:l);
data(2, i:i+l-1) = a_im(l, 1:l);
i += l;
}
}
}
} else {
/* Input is REAL. */
if (dims(1) == 2) {
if (dims(2) == 2) {
number = dims(3);
n = (long(sqrt(8*number + 1) + 0.5) - 1)/2;
if (number == n*(n + 1)/2) data = float(a);
}
} else if (dims(1) == 3) {
if (dims(2) == 2 && (n = dims(3)) == dims(4)) {
data = array(double, 2, n*(n + 1)/2);
i = 1;
for (l=1 ; l<=n ; ++l) {
data( , i:i+l-1) = a( , l, 1:l);
i += l;
}
data = float(data);
}
}
}
if (is_void(data)) error, "bad type/dimension list for input array";
lmax = n - 1;
return data;
}
func alm_unpack (ap)
/* DOCUMENT alm_unpack(ap)
Unpack a 2×NUMBER array AP into a 2×(LMAX+1)×(LMAX+1) array
where:
NUMBER = 1 + LMAX*(LMAX + 3)/2
In any case, the result is such that:
RESULT(1,..) = real part of the A(l,m)
RESULT(2,..) = imaginary part of the A(l,m)
SEE ALSO: alm_pack, alm_read. */
{
dims = dimsof(ap);
if (structof(ap) == complex) error, "complex type not yet implemented";
if (dims(1) == 2 && dims(2) == 2) {
number = dims(3);
n = (long(sqrt(8*number + 1) + 0.5) - 1)/2;
if (number == n*(n + 1)/2) {
ap = double(ap);
a = array(double, 2, n, n);
i = 1;
for (l=1 ; l<=n ; ++l) {
a( , l, 1:l) = ap( , i:i+l-1);
i += l;
}
return a;
}
}
error, "bad dimension list for packed A(l,m) data";
}
/*---------------------------------------------------------------------------*/
func hmap_read (filename)
/* DOCUMENT hmap_read(filename)
Read Healpix map from FITS file FILENAME.
SEE ALSO: hmap_write. */
{
/* I mus do something more sophisticated... */
h = fits_open(filename, 'r');
fits_next_hdu, h;
if (fits_get(h, "XTENSION") != "BINTABLE")
error, "expecting a FITS binary table";
if (fits_toupper(fits_get(h, "PIXTYPE")) != "HEALPIX")
error, "expecting PIXTYPE=HEALPIX";
if (fits_toupper(fits_get(h, "ORDERING")) != "RING")
error, "expecting ORDERING=RING";
if (is_void((nside = fits_get(h, "NSIDE"))) || nside <= 0)
error, "missing/invalid NSIDE";
if (is_void((firstpix = fits_get(h, "FIRSTPIX"))) || firstpix != 0)
error, "missing/invalid FIRSTPIX";
if (is_void((lastpix = fits_get(h, "LASTPIX"))) || lastpix < firstpix)
error, "missing/invalid LASTPIX";
if (is_void((lmax = fits_get(h, "MAX-LPOL"))) || lmax < 0)
error, "missing/invalid LMAX";
if (nside != lmax/2) error, "expecting NSIDE = LMAX/2";
ptr = fits_read_bintable(h);
if (numberof(ptr) != 1) error, "bad number of columns in FITS binary table";
return (*ptr(1))(firstpix+1:lastpix+1);
}
func hmap_write (filename, map, overwrite=)
/* DOCUMENT alm_save_map, filename, map;
Write Healpix map into FITS file FILENAME.
SEE ALSO: hmap_read. */
{
require, "fits2.i";
/* Compute NSIDE and LMAX based on:
* numberof(MAP) = NPIXTOT = 12*NSIDE*NSIDE
* NSIDE = LMAX/2
*/
if (structof(map) == complex) error, "complex type not yet implemented";
npix = numberof(map);
nside = long(sqrt(npix/12.0) + 0.5);
if (npix != 12*nside*nside) error, "bad number of elements in map";
lmax = 2*nside;
/* Put the map into a 1024×N array. */
m = 1024;
(data = array(float, m, (npix + m - 1)/m))(1:npix) = map(*);
/* create new FITS file with empty primary HDU */
h = fits_create(filename, extend=1, overwrite=overwrite);
fits_write_header, h;
/* append new HDU with BINTABLE data */
fits_new_bintable, h, "A(l,m) data for Healpix";
fits_set,h,"COMMENT" ,"-----------------------------------------------";
fits_set,h,"COMMENT" ," Sky Map Pixelisation Specific Keywords ";
fits_set,h,"COMMENT" ,"-----------------------------------------------";
fits_set,h,"PIXTYPE" ,"HEALPIX","HEALPIX Pixelisation";
fits_set,h,"ORDERING","RING","Pixel ordering scheme, either RING or NESTED";
fits_set,h,"NSIDE" ,nside,"Resolution parameter for HEALPIX";
fits_set,h,"FIRSTPIX",0,"First pixel # (0 based)";
fits_set,h,"LASTPIX", npix-1,"Last pixel # (0 based)";
fits_set,h,"COMMENT", "";
fits_set,h,"COMMENT", "-----------------------------------------------";
fits_set,h,"COMMENT", " Planck Simulation Specific Keywords ";
fits_set,h,"COMMENT", "-----------------------------------------------";
fits_set,h,"EXTNAME", "SIMULATION";
fits_set,h,"CREATOR", "Yorick", "Software creating the FITS file";
fits_set,h,"MAX-LPOL",lmax,"Maximum multipole L used in map synthesis";
fits_set,h,"COMMENT", "";
fits_set,h,"TTYPE1", "TEMPERATURE","temperature map";
fits_set,h,"TUNIT1", "K", "physical unit of map";
/* Write binary table header and data. */
fits_write_bintable, h, &data;
return h;
}
/*---------------------------------------------------------------------------*/
|