Healpix C++  3.50
compute_weights_module.cc
1 /*
2  * This file is part of Healpix_cxx.
3  *
4  * Healpix_cxx is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * Healpix_cxx is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with Healpix_cxx; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  *
18  * For more information about HEALPix, see http://healpix.sourceforge.net
19  */
20 
21 /*
22  * Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
23  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
24  * (DLR).
25  */
26 
27 /*
28  * Copyright (C) 2016-2018 Max-Planck-Society
29  * Author: Martin Reinecke
30  */
31 
32 #include "paramfile.h"
33 #include "fitshandle.h"
34 #include "announce.h"
35 #include "weight_utils.h"
36 
37 using namespace std;
38 
39 int compute_weights_module (int argc, const char **argv)
40  {
41  module_startup ("compute_weights", argc, argv);
42  paramfile params (getParamsFromCmdline(argc,argv));
43 
44  planck_assert (params.param_present("outfile_ring")
45  || params.param_present("outfile_full"), "no job requested");
46  int nside = params.find<int>("nside");
47  int nlmax = params.find<int>("nlmax");
48  if (nlmax&1)
49  {
50  cout << "Warning: specified nlmax is odd. Reducing by 1" << endl;
51  --nlmax;
52  }
53  double epsilon = params.find<double>("epsilon",1e-7);
54  int itmax = params.find<int>("max_iter",10000);
55 
56  if (params.param_present("outfile_ring"))
57  {
58  double epsilon_out;
59  vector<double> wgt=get_ringweights(nside,nlmax,epsilon,itmax,epsilon_out);
60  fitshandle out;
61  out.create (params.find<string>("outfile_ring"));
62  vector<fitscolumn> cols;
63  int repcount = 1;
64  cols.push_back (fitscolumn ("TEMPERATURE WEIGHTS","1",repcount,
65  PLANCK_FLOAT64));
66  cols.push_back (fitscolumn ("Q-POLARISATION WEIGHTS","1",repcount,
67  PLANCK_FLOAT64));
68  cols.push_back (fitscolumn ("U-POLARISATION WEIGHTS","1",repcount,
69  PLANCK_FLOAT64));
70  out.insert_bintab(cols);
71  out.set_key ("CREATOR",string("compute_weights"),
72  "Software creating the FITS file");
73  out.set_key ("NSIDE",nside,"Resolution parameter for HEALPIX");
74  out.set_key ("MAX_LPOL",nlmax,"Maximum multipole l used in map synthesis");
75  double val=*max_element(wgt.begin(),wgt.end());
76  out.set_key("MAXVAL1",val,"maximum value of T weights");
77  out.set_key("MAXVAL2",val,"maximum value of Q weights");
78  out.set_key("MAXVAL3",val,"maximum value of U weights");
79  val=*min_element(wgt.begin(),wgt.end());
80  out.set_key("MINVAL1",val,"minimum value of T weights");
81  out.set_key("MINVAL2",val,"minimum value of Q weights");
82  out.set_key("MINVAL3",val,"minimum value of U weights");
83  out.set_key("EPSILON",epsilon_out,"epsilon reached after minimization");
84  out.write_column(1,wgt);
85  out.write_column(2,wgt);
86  out.write_column(3,wgt);
87  }
88  if (params.param_present("outfile_full"))
89  {
90  double epsilon_out;
91  vector<double> wgt=get_fullweights(nside,nlmax,epsilon,itmax,epsilon_out);
92  fitshandle out;
93  out.create (params.find<string>("outfile_full"));
94  vector<fitscolumn> cols;
95  int repcount = 1;
96  cols.push_back (fitscolumn ("COMPRESSED PIXEL WEIGHTS","1",repcount,
97  PLANCK_FLOAT64));
98  out.insert_bintab(cols);
99  out.set_key ("CREATOR",string("compute_weights"),
100  "Software creating the FITS file");
101  out.set_key ("NSIDE",nside,"Resolution parameter for HEALPix");
102  out.set_key ("MAX_LPOL",nlmax,"Maximum l multipole");
103  double val=*max_element(wgt.begin(),wgt.end());
104  out.set_key("MAXVAL1",val,"maximum value of pixel weights");
105  val=*min_element(wgt.begin(),wgt.end());
106  out.set_key("MINVAL1",val,"minimum value of pixel weights");
107  out.set_key("EPSILON",epsilon_out,"epsilon reached after minimization");
108  out.write_column(1,wgt);
109  }
110  return 0;
111  }
void create(const std::string &fname)
vector< double > get_fullweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)
void insert_bintab(const std::vector< fitscolumn > &cols, const std::string &extname="xtension")
void write_column(int colnum, const arr< T > &data, int64 offset=0)
void set_key(const std::string &name, const T &value, const std::string &comment="")
paramfile getParamsFromCmdline(int argc, const char **argv, bool verbose=true)
void module_startup(const std::string &name, bool argc_valid, const std::string &usage, bool verbose=true)
#define planck_assert(testval, msg)
vector< double > get_ringweights(int nside, int lmax, double epsilon, int itmax, double &epsilon_out)

Generated on Mon Dec 10 2018 10:24:22 for Healpix C++