powspec.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.h"
00033 
00034 using namespace std;
00035 
00036 void PowSpec::dealloc()
00037   {
00038   tt_.dealloc();
00039   gg_.dealloc();
00040   cc_.dealloc();
00041   tg_.dealloc();
00042   tc_.dealloc();
00043   gc_.dealloc();
00044   }
00045 
00046 void PowSpec::Set(arr<double> &tt_new)
00047   {
00048   dealloc();
00049   num_specs = 1;
00050   tt_.transfer(tt_new);
00051 //FIXME: temporarily relaxed to allow cross-spectra
00052   for (int l=0; l<tt_.size(); ++l)
00053     if (tt_[l]<0)
00054       {
00055       cerr << "Warning: negative values in TT spectrum" << endl;
00056       break;
00057       }
00058   }
00059 
00060 void PowSpec::Set(arr<double> &tt_new, arr<double> &gg_new,
00061   arr<double> &cc_new, arr<double> &tg_new)
00062   {
00063   dealloc();
00064   num_specs = 4;
00065   tt_.transfer(tt_new);
00066   gg_.transfer(gg_new);
00067   cc_.transfer(cc_new);
00068   tg_.transfer(tg_new);
00069   planck_assert((tt_.size()==gg_.size()) && (tt_.size()==cc_.size())
00070     && (tt_.size()==tg_.size()), "PowSpec::Set: size mismatch");
00071   for (int l=0; l<tt_.size(); ++l)
00072     {
00073     planck_assert (tt_[l]>=0, "negative TT spectrum at l="+dataToString(l));
00074     planck_assert (gg_[l]>=0, "negative GG spectrum at l="+dataToString(l));
00075     planck_assert (cc_[l]>=0, "negative CC spectrum at l="+dataToString(l));
00076     planck_assert (abs(tg_[l]<=sqrt(tt_[l]*gg_[l])),
00077       "Inconsistent T, E and TxE terms at l="+dataToString(l));
00078     }
00079   }
00080 
00081 void PowSpec::Smooth_with_Gauss (double fwhm)
00082   {
00083   planck_assert (num_specs<=4, "not yet implemented for num_specs>4");
00084   double sigma = fwhm*fwhm2sigma;
00085   double fact_pol = exp(2*sigma*sigma);
00086   for (int l=0; l<tt_.size(); ++l)
00087     {
00088     double f1 = exp(-.5*l*(l+1)*sigma*sigma);
00089     double f2 = f1*fact_pol;
00090     tt_[l] *= f1*f1;
00091     if (num_specs>1)
00092       {
00093       gg_[l] *= f2*f2;
00094       cc_[l] *= f2*f2;
00095       tg_[l] *= f1*f2;
00096       }
00097     }
00098   }