Healpix C++  3.50
alice3.cc
Go to the documentation of this file.
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 /*! \file alice3.cc
28  * Copyright (C) 2005-2015 David Larson, Max-Planck-Society
29  * \author David Larson \author Martin Reinecke
30  */
31 
32 #include <iostream>
33 #include "announce.h"
34 #include "paramfile.h"
35 #include "healpix_map_fitsio.h"
36 #include "lsconstants.h"
37 #include "arr.h"
38 #include "fitshandle.h"
39 #include "vec3.h"
40 #include "string_utils.h"
41 #include "alm.h"
42 #include "alm_healpix_tools.h"
43 #include "planck_rng.h"
44 #include "healpix_map.h"
45 
46 using namespace std;
47 
48 class PolarizationHolder
49  {
50  private:
51  Healpix_Map<float> Q, U;
52 
53  public:
54  // Load a polarized fits file, with Q and U as the second
55  // and third columns (the standard form).
56  void load(const std::string &filename)
57  {
58  read_Healpix_map_from_fits(filename, Q, 2, 2);
59  read_Healpix_map_from_fits(filename, U, 3, 2);
60 #if 0
61  Healpix_Map<float> mag(Q);
62  for (int i=0; i<Q.Npix(); ++i)
63  mag[i]=sqrt(Q[i]*Q[i] + U[i]*U[i]);
64  float mmin,mmax;
65  mag.minmax(mmin,mmax);
66  for (int i=0; i<Q.Npix(); ++i)
67  {
68  Q[i]/=mmax;
69  U[i]/=mmax;
70  }
71 #endif
72  }
73 
74  void getQU(const pointing &p, float &q, float &u) const
75  {
76  fix_arr<int,4> pix;
78  Q.get_interpol(p,pix,wgt);
79  q=u=0.f;
80  for (tsize i=0;i<4;++i) {q+=Q[pix[i]]*wgt[i];u+=U[pix[i]]*wgt[i]; }
81  }
82 
83  vec3 getQUDir(const vec3 &loc) const
84  {
85  float q,u;
86  getQU(loc,q,u);
87  vec3 east(1,0,0);
88  if (abs(loc.x)+abs(loc.y) > 0.0)
89  east = vec3(-loc.y,loc.x,0).Norm();
90  vec3 north = crossprod(loc, east);
91  double angle = 0.5*safe_atan2(u,q);
92  return north*(-cos(angle)) + east*sin(angle);
93  }
94 
95  // Return the magnitude of the polarization at some pointing.
96  float getQUMagnitude(const pointing& p) const
97  {
98  float q,u;
99  getQU(p,q,u);
100  return sqrt(q*q + u*u);
101  }
102  };
103 
104 /*! Steps from loc in direction dir for an angle theta and updates loc and dir. */
105 void get_step(const PolarizationHolder &ph, vec3 &loc, vec3 &dir, double theta)
106  {
107  loc=(loc+dir*theta).Norm();
108  vec3 tdir=ph.getQUDir(loc);
109  dir = (dotprod(dir,tdir)<0) ? -tdir : tdir;
110  }
111 
112 /*! Performs one Runge-Kutta second order step. Updates loc and dir. */
113 void runge_kutta_step(vec3 &loc, vec3 &dir, const PolarizationHolder &ph,
114  double theta)
115  {
116  // Take a half-theta step
117  vec3 tloc=loc;
118  get_step(ph, tloc, dir, theta/2.0);
119 
120  // Then take a full step with the new direction
121  get_step(ph, loc, dir, theta);
122  }
123 
124 /*! Second order Runge-Kutta integration on the sphere. Given a
125  starting location, a qu map of the sky, and a step size theta, this
126  subroutine returns an array of vectors extending in both
127  directions from the starting location. */
128 void runge_kutta_2(const vec3 &location, const PolarizationHolder &ph,
129  double theta, arr<vec3> &locs)
130  {
131  vec3 first_dir=ph.getQUDir(location);
132  vec3 dir = first_dir;
133  vec3 loc = location;
134 
135  locs[locs.size()/2] = loc;
136 
137  for(int i = 1 + locs.size()/2; i<int(locs.size()); i++)
138  {
139  runge_kutta_step(loc, dir, ph, theta);
140  locs[i] = loc;
141  }
142 
143  dir = -first_dir;
144  loc = location;
145  for(int i = -1 + locs.size()/2; i>=0; i--)
146  {
147  runge_kutta_step(loc, dir, ph, theta);
148  locs[i] = loc;
149  }
150  }
151 
152 /*! Create a sinusoidal kernel. */
154  {
155  for(tsize i=0; i<kernel.size(); i++)
156  {
157  double sinx = sin(pi*(i+1) / (kernel.size()+1));
158  kernel[i] = sinx*sinx;
159  }
160  }
161 
162 /*! Convolve an array with a kernel. */
163 void convolve(const arr<double> &kernel, const arr<double> &raw, arr<double> &convolution)
164  {
165  convolution.alloc(raw.size()-kernel.size()+1);
166  for(tsize i=0; i<convolution.size(); i++)
167  {
168  double total=0;
169  for (tsize j=0; j<kernel.size(); j++)
170  total += kernel[j] * raw[i+j];
171  convolution[i] = total;
172  }
173  }
174 
175 int lic_function(Healpix_Map<float> &hitcount, Healpix_Map<float> &texture,
176  const PolarizationHolder &ph, const Healpix_Map<float> &th, int steps,
177  int kernel_steps, double step_radian)
178  {
179  arr<double> kernel(kernel_steps), convolution, rawtexture;
180  make_kernel(kernel);
181  arr<vec3> curve(steps);
182 
183  texture.fill(0.);
184  int num_curves=0;
185 
186  for(int i=0; i<texture.Npix(); i++)
187  {
188  if (hitcount[i]<1.0)
189  {
190  num_curves++;
191  runge_kutta_2(texture.pix2vec(i), ph, step_radian, curve);
192  rawtexture.alloc(curve.size());
193  for (tsize i2=0; i2<curve.size(); i2++)
194  rawtexture[i2] = th.interpolated_value(curve[i2]);
195  convolve(kernel, rawtexture, convolution);
196  for (tsize j=0; j<convolution.size(); j++)
197  {
198  int k = texture.vec2pix(curve[j+kernel.size()/2]);
199  texture[k] += convolution[j];
200  hitcount[k] += 1.;
201  }
202  }
203  }
204  return num_curves;
205  }
206 
207 int main(int argc, const char** argv)
208  {
209  module_startup ("alice3", argc, argv);
210  paramfile params (getParamsFromCmdline(argc,argv));
211 
212  PolarizationHolder ph;
213  ph.load(params.find<string>("in"));
214 
215  int nside = params.find<int>("nside");
216  int steps = params.find<int>("steps",100);
217  double step_radian=arcmin2rad*params.find<double>("step_arcmin",10.);
218  int kernel_steps = params.find<int>("kernel_steps",50);
219  float polmin = params.find<float>("polmin",-1e30),
220  polmax = params.find<float>("polmax",1e30);
221  string out = params.find<string>("out");
222 
223  Healpix_Map<float> th(nside,RING,SET_NSIDE);
224  if (params.param_present("texture"))
225  read_Healpix_map_from_fits(params.find<string>("texture"),th);
226  else
227  {
228  planck_rng rng(params.find<int>("rand_seed",42));
229  if (params.param_present("ell"))
230  {
231  int ell = params.find<int>("ell");
232  if (ell<0) ell=2*nside;
233  Alm<xcomplex<float> > a(ell,ell);
234  a.SetToZero();
235  cout << "Background texture using ell = " << ell << endl;
236 
237  a(ell,0)=fcomplex(rng.rand_gauss(),0.);
238  for (int m=0; m<=ell; m++)
239  { a(ell,m).real(rng.rand_gauss()); a(ell,m).imag(rng.rand_gauss()); }
240  alm2map(a, th);
241  }
242  else
243  {
244  for (int i=0; i<th.Npix(); i++)
245  th[i] = rng.rand_uni() - 0.5;
246  }
247  }
248 
249  Healpix_Map<float> hit(nside,RING,SET_NSIDE),
250  tex(nside,RING,SET_NSIDE),
251  mag(nside,RING,SET_NSIDE);
252 
253  hit.fill(0.);
254 
255  for (int i=0; i<mag.Npix(); i++)
256  {
257  pointing p = mag.pix2ang(i);
258 
259  mag[i] = min(polmax,max(polmin,ph.getQUMagnitude(p)));
260  tex[i] = th.interpolated_value(p);
261  }
262 
263  write_Healpix_map_to_fits(out+"_background.fits",tex,PLANCK_FLOAT32);
264 
265  lic_function(hit, tex, ph, th, steps, kernel_steps, step_radian);
266 
267  for (int i=0; i<tex.Npix(); ++i)
268  tex[i]/=hit[i];
269  float tmin,tmax,mmin,mmax;
270  tex.minmax(tmin,tmax);
271  mag.minmax(mmin,mmax);
272  for (int i=0; i<tex.Npix(); ++i)
273  {
274  mag[i]*=(tex[i]-tmin);
275  tex[i]=1.0-(tex[i]-tmin)/(tmax-tmin);
276  }
277  mag.minmax(mmin,mmax);
278  for (int i=0; i<mag.Npix(); ++i)
279  mag[i]=1.0-(mag[i]-mmin)/(mmax-mmin);
280  write_Healpix_map_to_fits(out+"_texture.fits",tex,PLANCK_FLOAT32);
281  write_Healpix_map_to_fits(out+"_mod_texture.fits",mag,PLANCK_FLOAT32);
282  }
vec3 pix2vec(I pix) const
Definition: healpix_base.h:193
double rand_gauss()
I vec2pix(const vec3 &vec) const
Definition: healpix_base.h:162
void make_kernel(arr< double > &kernel)
Definition: alice3.cc:153
Definition: alm.h:88
void runge_kutta_2(const vec3 &location, const PolarizationHolder &ph, double theta, arr< vec3 > &locs)
Definition: alice3.cc:128
void read_Healpix_map_from_fits(fitshandle &inp, Healpix_Map< T > &map, int colnum)
double safe_atan2(double y, double x)
paramfile getParamsFromCmdline(int argc, const char **argv, bool verbose=true)
std::size_t tsize
T interpolated_value(const pointing &ptg) const
Definition: healpix_map.h:218
void module_startup(const std::string &name, bool argc_valid, const std::string &usage, bool verbose=true)
vec3_t< T > crossprod(const vec3_t< T > &a, const vec3_t< T > &b)
void get_interpol(const pointing &ptg, fix_arr< I, 4 > &pix, fix_arr< double, 4 > &wgt) const
vec3_t< float64 > vec3
void write_Healpix_map_to_fits(fitshandle &out, const Healpix_Map< T > &map, PDT datatype)
T dotprod(const vec3_t< T > &v1, const vec3_t< T > &v2)
void alm2map(const Alm< xcomplex< T > > &alm, Healpix_Map< T > &map, bool add_map)
double rand_uni()
void get_step(const PolarizationHolder &ph, vec3 &loc, vec3 &dir, double theta)
Definition: alice3.cc:105
tsize size() const
void convolve(const arr< double > &kernel, const arr< double > &raw, arr< double > &convolution)
Definition: alice3.cc:163
void runge_kutta_step(vec3 &loc, vec3 &dir, const PolarizationHolder &ph, double theta)
Definition: alice3.cc:113
I Npix() const
Definition: healpix_base.h:429
void fill(const T &val)
Definition: healpix_map.h:90

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