48 class PolarizationHolder
56 void load(
const std::string &filename)
62 for (
int i=0; i<Q.
Npix(); ++i)
63 mag[i]=sqrt(Q[i]*Q[i] + U[i]*U[i]);
65 mag.minmax(mmin,mmax);
66 for (
int i=0; i<Q.
Npix(); ++i)
74 void getQU(
const pointing &p,
float &q,
float &u)
const 80 for (
tsize i=0;i<4;++i) {q+=Q[pix[i]]*wgt[i];u+=U[pix[i]]*wgt[i]; }
88 if (abs(loc.
x)+abs(loc.
y) > 0.0)
89 east =
vec3(-loc.
y,loc.
x,0).Norm();
92 return north*(-cos(angle)) + east*sin(angle);
96 float getQUMagnitude(
const pointing& p)
const 100 return sqrt(q*q + u*u);
107 loc=(loc+dir*theta).Norm();
108 vec3 tdir=ph.getQUDir(loc);
109 dir = (
dotprod(dir,tdir)<0) ? -tdir : tdir;
131 vec3 first_dir=ph.getQUDir(location);
132 vec3 dir = first_dir;
135 locs[locs.
size()/2] = loc;
137 for(
int i = 1 + locs.
size()/2; i<int(locs.
size()); i++)
145 for(
int i = -1 + locs.
size()/2; i>=0; i--)
157 double sinx = sin(pi*(i+1) / (kernel.
size()+1));
158 kernel[i] = sinx*sinx;
166 for(
tsize i=0; i<convolution.
size(); i++)
170 total += kernel[j] * raw[i+j];
171 convolution[i] = total;
177 int kernel_steps,
double step_radian)
179 arr<double> kernel(kernel_steps), convolution, rawtexture;
186 for(
int i=0; i<texture.
Npix(); i++)
192 rawtexture.
alloc(curve.size());
193 for (
tsize i2=0; i2<curve.size(); i2++)
195 convolve(kernel, rawtexture, convolution);
196 for (
tsize j=0; j<convolution.
size(); j++)
198 int k = texture.
vec2pix(curve[j+kernel.size()/2]);
199 texture[k] += convolution[j];
207 int main(
int argc,
const char** argv)
212 PolarizationHolder ph;
213 ph.load(params.find<
string>(
"in"));
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");
224 if (params.param_present(
"texture"))
228 planck_rng rng(params.find<
int>(
"rand_seed",42));
229 if (params.param_present(
"ell"))
231 int ell = params.find<
int>(
"ell");
232 if (ell<0) ell=2*nside;
235 cout <<
"Background texture using ell = " << ell << endl;
238 for (
int m=0; m<=ell; m++)
244 for (
int i=0; i<th.
Npix(); i++)
250 tex(nside,
RING,SET_NSIDE),
251 mag(nside,
RING,SET_NSIDE);
255 for (
int i=0; i<mag.Npix(); i++)
259 mag[i] = min(polmax,max(polmin,ph.getQUMagnitude(p)));
265 lic_function(hit, tex, ph, th, steps, kernel_steps, step_radian);
267 for (
int i=0; i<tex.Npix(); ++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)
274 mag[i]*=(tex[i]-tmin);
275 tex[i]=1.0-(tex[i]-tmin)/(tmax-tmin);
277 mag.minmax(mmin,mmax);
278 for (
int i=0; i<mag.Npix(); ++i)
279 mag[i]=1.0-(mag[i]-mmin)/(mmax-mmin);
vec3 pix2vec(I pix) const
I vec2pix(const vec3 &vec) const
void make_kernel(arr< double > &kernel)
void runge_kutta_2(const vec3 &location, const PolarizationHolder &ph, double theta, arr< vec3 > &locs)
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)
T interpolated_value(const pointing &ptg) const
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
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)
void get_step(const PolarizationHolder &ph, vec3 &loc, vec3 &dir, double theta)
void convolve(const arr< double > &kernel, const arr< double > &raw, arr< double > &convolution)
void runge_kutta_step(vec3 &loc, vec3 &dir, const PolarizationHolder &ph, double theta)