powspec_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 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include "powspec_fitsio.h"
00033 #include "powspec.h"
00034 #include "fitshandle.h"
00035 
00036 using namespace std;
00037 
00038 void read_powspec_from_fits (const string &infile, PowSpec &powspec,
00039   int nspecs, int lmax)
00040   {
00041   planck_assert ((nspecs==1)||(nspecs==4), "wrong number of spectra");
00042   fitshandle inp;
00043   inp.open(infile);
00044   inp.goto_hdu(2);
00045   if (inp.key_present("PDMTYPE"))    // official Planck format
00046     {
00047     inp.assert_pdmtype("POWERSPEC");
00048     arr<double> tmp;
00049     arr<int> itmp;
00050     if ((inp.coltype(1)==TINT)||(inp.coltype(1)==TLONG))
00051       inp.read_entire_column(1,itmp);
00052     else
00053       {
00054       cerr << "Warning: column containing l values is not of integer type!"
00055            << endl;
00056       inp.read_entire_column(1,tmp);
00057       itmp.alloc(tmp.size());
00058       for (int m=0; m<tmp.size(); ++m) itmp[m] = nearest<int>(tmp[m]);
00059       }
00060 
00061     if (nspecs==1)
00062       {
00063       arr<double> tt(lmax+1);
00064       tt.fill(0);
00065       inp.read_entire_column(2,tmp);
00066       for (int m=0; m<itmp.size(); ++m)
00067         if (itmp[m]<=lmax) tt[itmp[m]] = tmp[m];
00068       powspec.Set(tt);
00069       }
00070     else
00071       {
00072       arr<double> tt(lmax+1), gg(lmax+1), cc(lmax+1), tg(lmax+1);
00073       tt.fill(0); gg.fill(0); cc.fill(0); tg.fill(0);
00074       inp.read_entire_column(2,tmp);
00075       for (int m=0; m<itmp.size(); ++m)
00076         if (itmp[m]<=lmax) tt[itmp[m]] = tmp[m];
00077       inp.read_entire_column(3,tmp);
00078       for (int m=0; m<itmp.size(); ++m)
00079         if (itmp[m]<=lmax) gg[itmp[m]] = tmp[m];
00080       inp.read_entire_column(4,tmp);
00081       for (int m=0; m<itmp.size(); ++m)
00082         if (itmp[m]<=lmax) cc[itmp[m]] = tmp[m];
00083       inp.read_entire_column(5,tmp);
00084       for (int m=0; m<itmp.size(); ++m)
00085         if (itmp[m]<=lmax) tg[itmp[m]] = tmp[m];
00086       powspec.Set(tt,gg,cc,tg);
00087       }
00088     }
00089   else
00090     {
00091     int lmax_file = inp.nelems(1)-1;
00092     if (lmax_file<lmax)
00093       cerr << "warning: lmax in file smaller than expected; padding with 0."
00094            << endl;
00095     int lmax_read = min (lmax,lmax_file);
00096     if (nspecs==1)
00097       {
00098       arr<double> tt(lmax+1);
00099       inp.read_column_raw (1,&tt[0],lmax_read+1);
00100       for (int l=lmax_read+1; l<=lmax; ++l)
00101         tt[l]=0;
00102       powspec.Set(tt);
00103       }
00104     else
00105       {
00106       arr<double> tt(lmax+1), gg(lmax+1), cc(lmax+1), tg(lmax+1);
00107       inp.read_column_raw (1,&tt[0],lmax_read+1);
00108       inp.read_column_raw (2,&gg[0],lmax_read+1);
00109       inp.read_column_raw (3,&cc[0],lmax_read+1);
00110       inp.read_column_raw (4,&tg[0],lmax_read+1);
00111       for (int l=lmax_read+1; l<=lmax; ++l)
00112         tt[l]=gg[l]=cc[l]=tg[l]=0;
00113       powspec.Set(tt,gg,cc,tg);
00114       }
00115     }
00116   }
00117 
00118 void write_powspec_to_fits (fitshandle &out,
00119   const PowSpec &powspec, int nspecs)
00120   {
00121   planck_assert ((nspecs==1)||(nspecs==4), "wrong number of spectra");
00122   vector<fitscolumn> cols;
00123   cols.push_back(fitscolumn("l","multipole order",1,TINT32BIT));
00124   cols.push_back(fitscolumn("Temperature C_l","Kelvin-squared",1,TDOUBLE));
00125   if (nspecs>1)
00126     {
00127     cols.push_back(fitscolumn("E-mode C_l","Kelvin-squared",1,TDOUBLE));
00128     cols.push_back(fitscolumn("B-mode C_l","Kelvin-squared",1,TDOUBLE));
00129     cols.push_back(fitscolumn("T-E cross-corr.","Kelvin-squared",1,TDOUBLE));
00130     }
00131   out.insert_bintab(cols);
00132   out.add_key("PDMTYPE",string("POWERSPEC"),"Planck data model type");
00133   arr<int> tmparr(powspec.Lmax()+1);
00134   for (int l=0; l<=powspec.Lmax(); ++l) tmparr[l]=l;
00135   out.write_column(1,tmparr);
00136   out.write_column(2,powspec.tt());
00137   if (nspecs>1)
00138     {
00139     out.write_column(3,powspec.gg());
00140     out.write_column(4,powspec.cc());
00141     out.write_column(5,powspec.tg());
00142     }
00143   }