alm_fitsio.cc

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Copyright (C) 2003, 2004, 2005 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include <string>
00033 #include "alm_fitsio.h"
00034 #include "alm.h"
00035 #include "fitshandle.h"
00036 #include "xcomplex.h"
00037 
00038 using namespace std;
00039 
00040 void get_almsize(fitshandle &inp, int &lmax, int &mmax)
00041   {
00042   if (inp.key_present("MAX-LPOL") && inp.key_present("MAX-MPOL"))
00043     {
00044     inp.get_key ("MAX-LPOL",lmax);
00045     inp.get_key ("MAX-MPOL",mmax);
00046     return;
00047     }
00048 
00049   int n_alms = inp.nelems(1);
00050   arr<int> index;
00051   const int chunksize=1024*256;
00052   int offset=0;
00053   lmax=-1;
00054   mmax=-1;
00055   while (offset<n_alms)
00056     {
00057     int ppix=min(chunksize,n_alms-offset);
00058     index.alloc(ppix);
00059     inp.read_column(1,index,offset);
00060 
00061     for (int i=0; i<ppix; ++i)
00062       {
00063       int l = isqrt(index[i]-1);
00064       int m = index[i] - l*l - l - 1;
00065       if (l>lmax) lmax=l;
00066       if (m>mmax) mmax=m;
00067       }
00068     offset+=chunksize;
00069     }
00070   }
00071 
00072 void get_almsize(const string &filename, int &lmax, int &mmax, int hdunum)
00073   {
00074   fitshandle inp;
00075   inp.open (filename);
00076   inp.goto_hdu(hdunum);
00077   get_almsize (inp, lmax, mmax);
00078   }
00079 
00080 void get_almsize_pol(const string &filename, int &lmax, int &mmax)
00081   {
00082   int tlmax, tmmax;
00083   fitshandle inp;
00084   inp.open (filename);
00085   lmax=mmax=0;
00086   for (int hdu=2; hdu<=4; ++hdu)
00087     {
00088     inp.goto_hdu(hdu);
00089     get_almsize (inp,tlmax,tmmax);
00090     if (tlmax>lmax) lmax=tlmax;
00091     if (tmmax>mmax) mmax=tmmax;
00092     }
00093   }
00094 
00095 template<typename T> void read_Alm_from_fits
00096   (fitshandle &inp, Alm<xcomplex<T> >&alms, int lmax, int mmax)
00097   {
00098   int n_alms = inp.nelems(1);
00099   arr<int> index;
00100   arr<T> re, im;
00101 
00102   alms.Set(lmax, mmax);
00103   alms.SetToZero();
00104   int max_index = lmax*lmax + lmax + mmax + 1;
00105   const int chunksize=1024*256;
00106   int offset=0;
00107   while (offset<n_alms)
00108     {
00109     int ppix=min(chunksize,n_alms-offset);
00110     index.alloc(ppix);
00111     re.alloc(ppix); im.alloc(ppix);
00112     inp.read_column(1,index,offset);
00113     inp.read_column(2,re,offset);
00114     inp.read_column(3,im,offset);
00115 
00116     for (int i=0; i<ppix; ++i)
00117       {
00118       if (index[i]>max_index) return;
00119 
00120       int l = isqrt(index[i]-1);
00121       int m = index[i] - l*l - l - 1;
00122       planck_assert(m>=0,"negative m encountered");
00123       planck_assert(l>=m, "wrong l,m combination");
00124       if ((l<=lmax) && (m<=mmax))
00125         alms(l,m).Set (re[i], im[i]);
00126       }
00127     offset+=chunksize;
00128     }
00129   }
00130 
00131 template void read_Alm_from_fits (fitshandle &inp,
00132   Alm<xcomplex<double> > &alms, int lmax, int mmax);
00133 template void read_Alm_from_fits (fitshandle &inp,
00134   Alm<xcomplex<float> > &alms, int lmax, int mmax);
00135 
00136 
00137 template<typename T> void read_Alm_from_fits
00138   (const string &filename, Alm<xcomplex<T> >&alms, int lmax, int mmax,
00139   int hdunum)
00140   {
00141   fitshandle inp;
00142   inp.open (filename);
00143   inp.goto_hdu(hdunum);
00144   read_Alm_from_fits(inp,alms,lmax,mmax);
00145   }
00146 
00147 template void read_Alm_from_fits (const string &filename,
00148   Alm<xcomplex<double> > &alms, int lmax, int mmax, int hdunum);
00149 template void read_Alm_from_fits (const string &filename,
00150   Alm<xcomplex<float> > &alms, int lmax, int mmax, int hdunum);
00151 
00152 
00153 template<typename T> void write_Alm_to_fits
00154   (fitshandle &out, const Alm<xcomplex<T> > &alms, int lmax, int mmax,
00155   int datatype)
00156   {
00157   vector<fitscolumn> cols;
00158   cols.push_back (fitscolumn("index","l*l+l+m+1",1,TINT32BIT));
00159   cols.push_back (fitscolumn("real","unknown",1,datatype));
00160   cols.push_back (fitscolumn("imag","unknown",1,datatype));
00161   out.insert_bintab(cols);
00162   arr<int> index;
00163   arr<double> re, im;
00164 
00165   int lm=alms.Lmax(), mm=alms.Mmax();
00166   int n_alms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
00167 
00168   int l=0, m=0;
00169   const int chunksize=1024*256;
00170   int offset=0;
00171   while (offset<n_alms)
00172     {
00173     int ppix=min(chunksize,n_alms-offset);
00174     index.alloc(ppix);
00175     re.alloc(ppix); im.alloc(ppix);
00176     for (int i=0; i<ppix; ++i)
00177       {
00178       index[i] = l*l + l + m + 1;
00179       if ((l<=lm) && (m<=mm))
00180         { re[i] = alms(l,m).re; im[i] = alms(l,m).im; }
00181       else
00182         { re[i] = 0; im[i] = 0; }
00183       ++m;
00184       if ((m>l) || (m>mmax)) { ++l; m=0; }
00185       }
00186     out.write_column(1,index,offset);
00187     out.write_column(2,re,offset);
00188     out.write_column(3,im,offset);
00189 
00190     offset+=chunksize;
00191     }
00192   out.add_key("MAX-LPOL",lmax,"highest l in the table");
00193   out.add_key("MAX-MPOL",mmax,"highest m in the table");
00194   }
00195 
00196 template void write_Alm_to_fits
00197   (fitshandle &out, const Alm<xcomplex<double> > &alms, int lmax,
00198    int mmax, int datatype);
00199 template void write_Alm_to_fits
00200   (fitshandle &out, const Alm<xcomplex<float> > &alms, int lmax,
00201    int mmax, int datatype);
00202 
00203 
00204 template<typename T> void write_compressed_Alm_to_fits
00205   (fitshandle &out, const Alm<xcomplex<T> > &alms, int lmax, int mmax,
00206   int datatype)
00207   {
00208   vector<fitscolumn> cols;
00209   cols.push_back (fitscolumn("index","l*l+l+m+1",1,TINT32BIT));
00210   cols.push_back (fitscolumn("real","unknown",1,datatype));
00211   cols.push_back (fitscolumn("imag","unknown",1,datatype));
00212   out.insert_bintab(cols);
00213   arr<int> index;
00214   arr<double> re, im;
00215 
00216   int n_alms = 0;
00217   for (int m=0; m<=mmax; ++m)
00218     for (int l=m; l<=lmax; ++l)
00219       if (alms(l,m).norm()>0) ++n_alms;
00220 
00221   int l=0, m=0;
00222   const int chunksize=1024*256;
00223   int real_lmax=0, real_mmax=0;
00224   int offset=0;
00225   while (offset<n_alms)
00226     {
00227     int ppix=min(chunksize,n_alms-offset);
00228     index.alloc(ppix);
00229     re.alloc(ppix); im.alloc(ppix);
00230     for (int i=0; i<ppix; ++i)
00231       {
00232       while (alms(l,m).norm()==0)
00233         {
00234         ++m;
00235         if ((m>l) || (m>mmax)) { ++l; m=0; }
00236         }
00237       index[i] = l*l + l + m + 1;
00238       re[i] = alms(l,m).re;
00239       im[i] = alms(l,m).im;
00240       if (l>real_lmax) real_lmax=l;
00241       if (m>real_mmax) real_mmax=m;
00242       ++m;
00243       if ((m>l) || (m>mmax)) { ++l; m=0; }
00244       }
00245     out.write_column(1,index,offset);
00246     out.write_column(2,re,offset);
00247     out.write_column(3,im,offset);
00248 
00249     offset+=chunksize;
00250     }
00251   out.add_key("MAX-LPOL",real_lmax,"highest l in the table");
00252   out.add_key("MAX-MPOL",real_mmax,"highest m in the table");
00253   }
00254 
00255 template void write_compressed_Alm_to_fits
00256   (fitshandle &out, const Alm<xcomplex<double> > &alms, int lmax,
00257    int mmax, int datatype);
00258 template void write_compressed_Alm_to_fits
00259   (fitshandle &out, const Alm<xcomplex<float> > &alms, int lmax,
00260    int mmax, int datatype);