GLAMERDOC++
Gravitational Lensing Code Library
Loading...
Searching...
No Matches
multimap.h
1/*
2 * MOKAlens.h
3 *
4 */
5
6#ifndef MultiLENS_H_
7#define MultiLENS_H_
8
9#include "standard.h"
10#include "profile.h"
11#include "InputParams.h"
12#include "lens_halos.h"
13#include "grid_maintenance.h"
14#include "cpfits.h"
15
16#include <stdexcept>
17
18// this is just to make sure no plan is executed with the wrong dimensions
19struct my_fftw_plan{
20 fftw_plan plan_r2c;
21 fftw_plan plan_c2r;
22 size_t nx;
23 size_t ny;
24
25 void initialize(size_t _nx_,size_t _ny_){
26 nx = _nx_;
27 ny = _ny_;
28
29 std::vector<double> rdummy(nx*ny);
30 size_t Nkx = (nx/2+1);
31 fftw_complex *cdummy = new fftw_complex[ny*Nkx];
32
33 plan_r2c = fftw_plan_dft_r2c_2d(ny,nx,rdummy.data(),cdummy
34 ,FFTW_ESTIMATE);
35 plan_c2r = fftw_plan_dft_c2r_2d(ny,nx,cdummy,rdummy.data()
36 ,FFTW_MEASURE);
37 delete[] cdummy;
38 }
39
40 ~my_fftw_plan(){
41 fftw_destroy_plan(plan_r2c);
42 fftw_destroy_plan(plan_c2r);
43 }
44};
45
56struct LensMap{
57
58 LensMap():nx(0),ny(0),boxlMpc(0),angular_pixel_size(0){};
59
60 // move operators
61 LensMap(LensMap &&m);
62 LensMap& operator=(LensMap &&m);
63
65 std::valarray<double> surface_density; // Msun / Mpc^2
66 std::valarray<float> alpha1_bar; // Msun / Mpc
67 std::valarray<float> alpha2_bar; // Msun / Mpc
68 std::valarray<float> gamma1_bar; // Msun / Mpc^2
69 std::valarray<float> gamma2_bar; // Msun / Mpc^2
70 std::valarray<float> phi_bar; // Msun
71 int nx,ny;
72
73 double boxlMpc; // range in Mpc
74 double angular_pixel_size; // in radians
75 Point_2d center;
76 Point_2d lowerleft;
77 Point_2d upperright;
78
80 double x_resolution(){return boxlMpc / nx ;}
81 double y_resolution(){return (upperright[1]-lowerleft[1])/ny;}
82 // # of pixels times resolution
83 double x_range(){return boxlMpc;}
84 // # of pixels times resolution
85 double y_range(){return (upperright[1]-lowerleft[1]);}
86
87 bool evaluate(const PosType *x,KappaType &sigma,KappaType *gamma,PosType *alpha);
88
89 LensMap(std::string fits_input_file,double angDist){
90 read(fits_input_file,angDist);
91 }
92
94 void read(std::string input_fits,double angDist);//,float h,float z);
96 void Myread(std::string fits_input_file);
97
98
100 //void read_header(std::string input_fits,float h,float z);
101 void read_header(std::string input_fits,double angDist);
102
104// void read_sub(std::string input_fits
105// ,const std::vector<long> &first
106// ,const std::vector<long> &last
107// ,double Dist
108// );
109
110 void read_sub(CPFITS_READ &cpfits
111 ,std::vector<long> &first
112 ,std::vector<long> &last
113 ,double Dist
114 );
115
116 void write(std::string filename);
118 void write(std::string filename,LensingVariable quant);
119
120 // this calculates the other lensing quantities from the density map
121
122 //std::mutex mutex_lensmap;
123
124 // these are thread safe versions
125 template <class T>
126 void ProcessFFTs(float zerosize
127 ,T Wphi_of_k
128 ,my_fftw_plan &plan_padded // the plan must be initialized with thr right size arrays
129 ,bool do_alpha = true);
130
131 template <class T>
132 void ProcessFFTs(T Wphi_of_k
133 ,my_fftw_plan &plan // the plan must be initialized with thr right size arrays
134 ,bool do_alpha = true);
135
136
137 template <class T>
138 void ProcessFFTs(float zerosize
139 ,T Wphi_of_k
140 ,bool do_alpha = true);
141
142 template <class T>
143 void ProcessFFTs(T Wphi_of_k
144 ,bool do_alpha = true);
145
146 void make_fftw_plans(
147 my_fftw_plan &plans
148 ,double zerosize
149 ){
150
151 size_t Nnx=int(zerosize*nx);
152 size_t Nny=int(zerosize*ny);
153
154 plans.nx = Nnx;
155 plans.ny = Nny;
156
157 std::vector<double> rdummy(Nnx*Nny);
158 size_t Nkx = (Nnx/2+1);
159 fftw_complex *cdummy = new fftw_complex[Nny*Nkx];
160
161 plans.plan_r2c = fftw_plan_dft_r2c_2d(Nny,Nnx,rdummy.data(),cdummy
162 ,FFTW_ESTIMATE);
163 plans.plan_c2r = fftw_plan_dft_c2r_2d(Nny,Nnx,cdummy,rdummy.data()
164 ,FFTW_MEASURE);
165
166 delete[] cdummy;
167 }
168};
169
176class LensHaloMultiMap : public LensHalo
177{
178public:
179
180 LensHaloMultiMap(
181 std::string fitsfile
182 ,std::string dir_data
183 ,double redshift
184 ,double mass_unit
185 ,double npane
186 ,int number_of_subfields
187 ,COSMOLOGY &c
188 ,bool write_subfields = false
189 ,std::string dir_scratch = ""
190 ,bool subtract_ave = true
191 ,double ffactor = 5 // coarse grid size in units of smoothing size
192 ,double gfactor = 5 // ratio of the border size to the smoothing scale
193 ,double rsmooth2 = 0 // smoothing scale on potential
194 );
195
196 ~LensHaloMultiMap(){
197 --count;
198 if(count == 0){
199 plans_set_up = false;
200 }
201 };
202
204
205 static my_fftw_plan plan_short_range;
206 //static fftw_plan plan_c2r_short_range;
207 static my_fftw_plan plan_long_range;
208 //static fftw_plan plan_c2r_long_range;
209 static std::mutex mutex_multimap;
210 static bool plans_set_up;
211 static size_t nx_sub; // size of used highres field
212 static size_t ny_sub;
213 static size_t nx_long; // size of used highres field
214 static size_t ny_long;
215 static size_t nx_sub_extended; // size of highres field with borders
216 static size_t ny_sub_extended;
217 static int count;
218
219 //const double ffactor = 5,gfactor = 5;
220 double ffactor;
221 double gfactor;
222 //const double ffactor = 10,gfactor = 10;
223
224
229 void resetsubmap(
230 int i
231 ,const std::vector<long> &lower_left_pix
232 );
233
235 void resetsubmapPhys(int i,Point_2d ll);
236 void resetsubmapAngular(int i,Point_2d ll){
237 if(i >= short_range_maps.size()) throw std::invalid_argument("");
238 double D = getDist();
239 resetsubmapPhys(i,ll*D);
240 }
241
242 void force_halo(double *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,double const *xcm,bool subtract_point=false,PosType screening = 1.0);
243
244 void writeImage(std::string fn);
245
247 Point_2d getLowerLeft_lr() const { return long_range_map.lowerleft; }
249 Point_2d getUpperRight_lr() const { return long_range_map.upperright; }
251 Point_2d getCenter_lr() const { return long_range_map.center; }
252
254 Point_2d getLowerLeft_sr(int i) const { return short_range_maps[i].lowerleft; }
256 Point_2d getUpperRight_sr(int i) const { return short_range_maps[i].upperright; }
258 Point_2d getCenter_sr(int i) const { return short_range_maps[i].center; }
259
261 int NumberOfShortRangeMaps(){return short_range_maps.size();}
262
264 double getRangeMpc_lr() const { return long_range_map.boxlMpc; }
266 double getRangeMpc_sr() const { return nx_sub*resolution_mpc; }
267
269 size_t getNx_lr() const { return long_range_map.nx; }
271 size_t getNy_lr() const { return long_range_map.ny; }
272
274 size_t getNx_sr() const { return nx_sub; }
276 size_t getNy_sr() const { return ny_sub; }
277
279 size_t getNx() const { return Noriginal[0]; }
281 size_t getNy() const { return Noriginal[1]; }
282
283 double getMax() const {return max_pix;}
284 double getMin() const {return min_pix;}
285 double getResolutionMpc() const {return resolution_mpc;}
286 double getResolutionAngular() const {return angular_resolution;}
287
288 void operator =(LensHaloMultiMap &&m){
289 LensHalo::operator=(std::move(m));
290 cosmo = m.cosmo;
291 long_range_map = std::move(m.long_range_map);
292 short_range_maps = std::move(m.short_range_maps);
293 //single_grid = m.single_grid;
294 cpfits = std::move(m.cpfits);
295 //m.cpfits = nullptr;
296 max_pix = m.max_pix;
297 min_pix = m.min_pix;
298 mass_unit = m.mass_unit;
299 Noriginal[0] = m.Noriginal[0];
300 Noriginal[1] = m.Noriginal[1];
301 resolution_mpc = m.resolution_mpc;
302 angular_resolution = m.angular_resolution;
303 border_width_pix = m.border_width_pix;
304 fitsfilename = m.fitsfilename;
305 rs2 = m.rs2;
306 rsmooth2 = m.rsmooth2;
307 zerosize = m.zerosize;
308 unit = m.unit;
309 wsr = m.wsr;
310 wlr = m.wlr;
311 wsr_smooth = m.wsr_smooth;
312 ave_ang_sd = m.ave_ang_sd;
313 subfield_filename = m.subfield_filename;
314 write_shorts = m.write_shorts;
315 ffactor = m.ffactor;
316 gfactor = m.gfactor;
317 }
318
319 LensHaloMultiMap(LensHaloMultiMap &&m):
320 LensHalo(std::move(m)),cosmo(m.cosmo),cpfits(std::move(m.cpfits))
321 {
322 long_range_map = std::move(m.long_range_map);
323 short_range_maps = std::move(m.short_range_maps);
324 //single_grid = m.single_grid;
325 //cpfits = std::move(m.cpfits);
326 //cpfits = m.cpfits;
327 //m.cpfits = nullptr;
328 max_pix = m.max_pix;
329 min_pix = m.min_pix;
330 mass_unit = m.mass_unit;
331 Noriginal[0] = m.Noriginal[0];
332 Noriginal[1] = m.Noriginal[1];
333 resolution_mpc = m.resolution_mpc;
334 angular_resolution = m.angular_resolution;
335 border_width_pix = m.border_width_pix;
336 fitsfilename = m.fitsfilename;
337 rs2 = m.rs2;
338 rsmooth2 = m.rsmooth2;
339 zerosize = m.zerosize;
340 unit = m.unit;
341 wsr = m.wsr;
342 wlr = m.wlr;
343 wsr_smooth = m.wsr_smooth;
344 ave_ang_sd = m.ave_ang_sd;
345 subfield_filename = m.subfield_filename;
346 write_shorts = m.write_shorts;
347 ffactor = m.ffactor;
348 gfactor = m.gfactor;
349 }
350
351public:
352 LensMap long_range_map;
353 std::vector<LensMap> short_range_maps;
354
355private:
356
357 bool write_shorts;
358 //bool single_grid;
359 COSMOLOGY &cosmo;
360 CPFITS_READ cpfits;
361
362 double ave_ang_sd;
363
364 double max_pix = std::numeric_limits<double>::lowest();
365 double min_pix = std::numeric_limits<double>::max();
366 double rsmooth2;
367
368 double mass_unit;
369
370 size_t Noriginal[2]; // number of pixels in each dimension in original image
371 double resolution_mpc; // resolution of original image and short range image in Mpc
372 double angular_resolution; // angular resolution of original image
373 long border_width_pix; // width of short range maps padding
374 std::string fitsfilename;
375 std::string subfield_filename;
376
377 double rs2;
378
379 //const COSMOLOGY& cosmo;
380 int zerosize;
381
382 // setsup everything for the given short range map.
383 void setsubmap(LensMap &short_range_map
384 ,const std::vector<long> &lower_left
385 //,const std::vector<long> &upper_right
386 );
387 struct UNIT{
388 int operator()(float k2){return 1;}
389 };
390 struct WLR{
391 float rs2;
392 float operator()(float k2){return exp(-k2*rs2);}
393 };
394 struct WSR{
395 float rs2;
396 float operator()(float k2){return 1 - exp(-k2*rs2);}
397 };
398
399 struct WSRSMOOTH {
400 float rs2;
401 float r_sm2; // smoothing scale
402 float operator()(float k2){return (1 - exp(-k2*rs2)) * exp(-k2*r_sm2);}
403 };
404
405 UNIT unit;
406 WSRSMOOTH wsr_smooth;
407 WSR wsr;
408 WLR wlr;
409};
410
411
416template <typename T>
417void LensMap::ProcessFFTs(
418 float zerosize
419 ,T Wphi_of_k
420 ,bool do_alpha){
421 my_fftw_plan plan_padded;
422
423 make_fftw_plans(plan_padded,zerosize);
424
425 ProcessFFTs<T>(zerosize,Wphi_of_k,plan_padded,do_alpha);
426}
427
429template <typename T>
430void LensMap::ProcessFFTs(
431 float zerosize
432 ,T Wphi_of_k
433 ,my_fftw_plan &plan_padded // the plan must be initialized with thr right size arrays
434 ,bool do_alpha){
435
436 assert(surface_density.size() == nx*ny);
437
438 // size of the new map in x and y directions, factor by which each size is increased
439 int Nnx=int(zerosize*nx);
440 int Nny=int(zerosize*ny);
441 size_t NN = Nnx*Nny;
442 double boxlx = boxlMpc*zerosize;
443 double boxly = ny*boxlMpc*zerosize/nx;
444
445 assert( (plan_padded.nx == Nnx)*(plan_padded.ny == Nny) );
446
447 int imin = (Nnx-nx)/2;
448 int imax = (Nnx+nx)/2;
449 int jmin = (Nny-ny)/2;
450 int jmax = (Nny+ny)/2;
451
452 size_t Nkx = (Nnx/2+1);
453
454 double tmp = 2.*M_PI/boxlx;
455 std::vector<double> kxs(Nkx);
456 for( int i=0; i<Nkx; i++ ){
457 kxs[i] = i*tmp;
458 }
459 tmp = 2.*M_PI/boxly;
460 std::vector<double> kys(Nny);
461 for( int j=0; j<Nny; j++ ){
462 kys[j]=(j<Nny/2)?double(j):double(j-Nny);
463 kys[j] *= tmp;
464 }
465
466 std::vector<double> extended_map( NN );
467
468 // assume locate in a rectangular map and build up the new one
469 for( int j=0; j<Nny; j++ ){
470 for( int i=0; i<Nnx; i++ ){
471 if(i>=imin && i<imax && j>=jmin && j<jmax){
472 int ii = i-imin;
473 int jj = j-jmin;
474
475 if(ii>=nx || jj>=ny){
476 std::cerr << " 1 error mapping " << ii << " " << jj << std::endl;
477 exit(1);
478 }
479 if(ii<0 || jj<0){
480 std::cerr << " 2 error mapping " << ii << " " << jj << std::endl;
481 exit(1);
482 }
483 assert(ii+nx*jj < surface_density.size());
484 assert(i+Nnx*j < extended_map.size());
485 extended_map[i+Nnx*j] = surface_density[ii+nx*jj];
486 //float tmp = extended_map[i+Nnx*j];
487 assert(!isnan(extended_map[i+Nnx*j]));
488 if(isinf(extended_map[i+Nnx*j])){
489 extended_map[i+Nnx*j] = 0;
490 }
491 }else{
492 extended_map[i+Nnx*j] = 0;
493 }
494
495 }
496 }
497
498 //std::vector<fftw_complex> fNmap(Nny*(Nnx/2+1));
499 //std::vector<fftw_complex> fphi( Nny*(Nnx/2+1) );
500 fftw_complex *fphi = new fftw_complex[Nny*Nkx];
501
502 fftw_execute_dft_r2c(plan_padded.plan_r2c
503 ,extended_map.data(),fphi);
504
505
506 // fourier space
507 // std:: cout << " allocating fourier space maps " << std:: endl;
508
509 // build modes for each pixel in the fourier space
510 for( int i=0; i<Nkx; i++ ){
511 for( int j=0; j<Nny; j++ ){
512
513 double k2 = kxs[i]*kxs[i] + kys[j]*kys[j];
514 size_t k = i+(Nkx)*j;
515
516 // null for k2 = 0 no divergence
517 if(k2 == 0){
518 fphi[k][0] = 0.;
519 fphi[k][1] = 0.;
520 }else{
521
522 //assert(k < Nny*Nkx);
523 //assert(!isnan(fphi[k][0]));
524
525 // fphi
526 //fphi[i+(Nkx)*j][0]= -2.*fNmap[i+(Nkx)*j][0]/k2;
527 //fphi[i+(Nkx)*j][1]= -2.*fNmap[i+(Nkx)*j][1]/k2;
528
529 // fphi
530 fphi[k][0] *= -2./k2;
531 fphi[k][1] *= -2./k2;
532
533 // apply window function
534 double w = Wphi_of_k(k2);
535 fphi[k][0] *= w;
536 fphi[k][1] *= w;
537 }
538 }
539 }
540
541 fftw_complex *fft= new fftw_complex[Nny*(Nkx)];
542 //double *realsp = new double[Nnx*Nny];
543
544 //std::vector<fftw_complex> fft( Nny*(Nkx) );
545 std::vector<double> realsp(Nnx*Nny);
546
547 if(do_alpha){
548 // alpha1
549 {
550 // build modes for each pixel in the fourier space
551 for( int i=0; i<Nkx; i++ ){
552 for( int j=0; j<Nny; j++ ){
553
554 size_t k = i+(Nkx)*j;
555 fft[k][0] = -kxs[i]*fphi[k][1];
556 fft[k][1] = kxs[i]*fphi[k][0];
557 //assert(!isnan(fft[k][0]));
558 }
559 }
560
561 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
562
563 alpha1_bar.resize(nx*ny,0);
564
565 for( int j=jmin; j<jmax; j++ ){
566 int jj = j-jmin;
567 for( int i=imin; i<imax; i++ ){
568 int ii = i-imin;
569
570 alpha1_bar[ii+nx*jj] = -1*float(realsp[i+Nnx*j]/NN);
571 }
572 }
573 }
574
575 // alpha2
576 {
577 // build modes for each pixel in the fourier space
578 for( int j=0; j<Nny; j++ ){
579 for( int i=0; i<Nkx; i++ ){
580 size_t k = i+(Nkx)*j;
581
582 // alpha
583 fft[k][0] = -kys[j]*fphi[k][1];
584 fft[k][1] = kys[j]*fphi[k][0];
585
586 }
587 }
588
589 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
590
591 alpha2_bar.resize(nx*ny,0);
592
593 for( int j=jmin; j<jmax; j++ ){
594 int jj = j-jmin;
595 for( int i=imin; i<imax; i++ ){
596 int ii = i-imin;
597
598 alpha2_bar[ii+nx*jj] = -1*float(realsp[i+Nnx*j]/NN);
599 }
600 }
601 }
602 }
603 // gamma1
604 {
605 // build modes for each pixel in the fourier space
606 for( int i=0; i<Nkx; i++ ){
607 for( int j=0; j<Nny; j++ ){
608
609 size_t k = i+(Nkx)*j;
610 // gamma
611 fft[k][0] = 0.5*(kxs[i]*kxs[i]-kys[j]*kys[j])*fphi[k][0];
612 fft[k][1] = 0.5*(kxs[i]*kxs[i]-kys[j]*kys[j])*fphi[k][1];
613 }
614 }
615
616 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
617
618 gamma1_bar.resize(nx*ny,0);
619
620 for( int j=jmin; j<jmax; j++ ){
621 int jj = j-jmin;
622 for( int i=imin; i<imax; i++ ){
623 int ii = i-imin;
624
625 gamma1_bar[ii+nx*jj] = float( realsp[i+Nnx*j]/NN);
626 }
627 }
628 }
629 // gamma2
630 {
631 // build modes for each pixel in the fourier space
632 for( int i=0; i<Nkx; i++ ){
633 for( int j=0; j<Nny; j++ ){
634
635 size_t k = i+(Nkx)*j;
636
637 // gamma
638 fft[k][0] = kxs[i]*kys[j]*fphi[k][0];
639 fft[k][1] = kxs[i]*kys[j]*fphi[k][1];
640
641 }
642 }
643
644 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
645
646 gamma2_bar.resize(nx*ny);
647
648 for( int j=0; j<ny; j++ ){
649 int jj = j+jmin;
650 for( int i=0; i<nx; i++ ){
651 int ii = i+imin;
652
653 gamma2_bar[i+nx*j] = float(-realsp[ii+Nnx*jj]/NN);
654 }
655 }
656 }
657
658 // kappa - this is done over because of the window in Fourier space
659 {
660
661 // build modes for each pixel in the fourier space
662 for( int i=0; i<Nkx; i++ ){
663 for( int j=0; j<Nny; j++ ){
664
665 double k2 = -(kxs[i]*kxs[i] + kys[j]*kys[j])/2;
666
667 size_t k = i+(Nkx)*j;
668
669 // surface density
670 fft[k][0] = k2*fphi[k][0];
671 fft[k][1] = k2*fphi[k][1];
672
673 }
674 }
675
676 assert( (plan_padded.nx == Nnx)*(plan_padded.ny == Nny) );
677 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
678
679 for( int j=jmin; j<jmax; j++ ){
680 int jj = j-jmin;
681 for( int i=imin; i<imax; i++ ){
682 int ii = i-imin;
683
684 surface_density[ii+nx*jj] = float(realsp[i+Nnx*j]/NN);
685
686 }
687 }
688 }
689
690 // phi
691 {
692
693 // build modes for each pixel in the fourier space
694 for( int i=0; i<Nkx; i++ ){
695 for( int j=0; j<Nny; j++ ){
696 size_t k = i+(Nkx)*j;
697
698 // surface density
699 fft[k][0] = fphi[k][0];
700 fft[k][1] = fphi[k][1];
701
702 }
703 }
704
705 assert( (plan_padded.nx == Nnx)*(plan_padded.ny == Nny) );
706 fftw_execute_dft_c2r(plan_padded.plan_c2r,fft,realsp.data());
707
708 phi_bar.resize(nx*ny);
709 for( int j=jmin; j<jmax; j++ ){
710 int jj = j-jmin;
711 for( int i=imin; i<imax; i++ ){
712 int ii = i-imin;
713
714 phi_bar[ii+nx*jj] = float(realsp[i+Nnx*j]/NN);
715
716 }
717 }
718 }
719
720 // std:: cout << " remapping the map in the original size " << std:: endl;
721 delete[] fft;
722 delete[] fphi;
723}
724
726template <class T>
727 void LensMap::ProcessFFTs(T Wphi_of_k
728 ,bool do_alpha){
729 my_fftw_plan plan;
730
731 make_fftw_plans(plan,1.0);
732 ProcessFFTs<T>(Wphi_of_k,plan,do_alpha);
733 }
734
735template <class T>
736void LensMap::ProcessFFTs(T Wphi_of_k
737 ,my_fftw_plan &plan // the plan must be initialized with thr right size arrays
738 ,bool do_alpha){
739
740
741 assert(surface_density.size() == nx*ny);
742 assert( (plan.nx == nx)*(plan.ny == ny) );
743
744 // size of the new map in x and y directions, factor by which each size is increased
745
746 size_t Nkx = (nx/2+1);
747 size_t NN = nx*ny;
748
749 double tmp = 2.*M_PI/boxlMpc;
750
751 std::vector<double> kxs(Nkx);
752 for( int i=0; i<Nkx; i++ ){
753 kxs[i] = i*tmp;
754 }
755 tmp = 2.*M_PI * nx /boxlMpc / ny;
756 std::vector<double> kys(ny);
757 for( int j=0; j<ny; j++ ){
758 kys[j]=(j<ny/2)?double(j):double(j-ny);
759 kys[j] *= tmp;
760 }
761
762 //std::vector<double> extended_map( Nnx*Nny );
763 fftw_complex *fphi = new fftw_complex[ny*Nkx];
764 fftw_complex *fft= new fftw_complex[ny*(Nkx)];
765 std::vector<double> realsp(NN);
766
767 /*
768 fftw_plan plan_r2c,plan_c2r;
769 {
770 std::lock_guard<std::mutex> hold(mu);
771 plan_r2c = fftw_plan_dft_r2c_2d(ny,nx,&(surface_density[0])
772 ,fphi,FFTW_ESTIMATE);
773 plan_c2r = fftw_plan_dft_c2r_2d(ny,nx,fft,realsp.data(),FFTW_MEASURE);
774 }
775*/
776
777 fftw_execute_dft_r2c(plan.plan_r2c,&(surface_density[0]),fphi);
778
779 // fourier space
780 // std:: cout << " allocating fourier space maps " << std:: endl;
781
782 // build modes for each pixel in the fourier space
783 for( int i=0; i<Nkx; i++ ){
784 for( int j=0; j<ny; j++ ){
785
786 double k2 = kxs[i]*kxs[i] + kys[j]*kys[j];
787 size_t k = i+(Nkx)*j;
788
789 // null for k2 = 0 no divergence
790 if(k2 == 0){
791 fphi[k][0] = 0.;
792 fphi[k][1] = 0.;
793 }else{
794
795 //assert(k < Nny*Nkx);
796 //assert(!isnan(fphi[k][0]));
797
798 // fphi
799 //fphi[i+(Nkx)*j][0]= -2.*fNmap[i+(Nkx)*j][0]/k2;
800 //fphi[i+(Nkx)*j][1]= -2.*fNmap[i+(Nkx)*j][1]/k2;
801
802 // fphi
803 fphi[k][0] *= -2./k2;
804 fphi[k][1] *= -2./k2;
805
806 // apply window function
807 double w = Wphi_of_k(k2);
808 fphi[k][0] *= w;
809 fphi[k][1] *= w;
810 }
811 }
812 }
813
814 if(do_alpha){
815 // alpha1
816 {
817
818 // build modes for each pixel in the fourier space
819 for( int i=0; i<Nkx; i++ ){
820 for( int j=0; j<ny; j++ ){
821
822 size_t k = i + Nkx * j;
823 fft[k][0] = -kxs[i]*fphi[k][1];
824 fft[k][1] = kxs[i]*fphi[k][0];
825 //assert(!isnan(fft[k][0]));
826 }
827 }
828
829 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
830
831 alpha1_bar.resize(NN);
832 for( size_t i=0; i<NN; i++ ) alpha1_bar[i] = -1*float(realsp[i]/NN);
833 }
834
835 // alpha2
836 {
837
838 // build modes for each pixel in the fourier space
839 for( int j=0; j<ny; j++ ){
840 for( int i=0; i<Nkx; i++ ){
841 size_t k = i+(Nkx)*j;
842
843 // alpha
844 fft[k][0] = -kys[j]*fphi[k][1];
845 fft[k][1] = kys[j]*fphi[k][0];
846
847 }
848 }
849
850 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
851
852 alpha2_bar.resize(NN);
853 for( size_t i=0; i<NN; i++ ) alpha2_bar[i] = -1*float(realsp[i]/NN);
854 }
855
856 // phi
857 {
858
859 // build modes for each pixel in the fourier space
860 for( int i=0; i<Nkx; i++ ){
861 for( int j=0; j<ny; j++ ){
862 size_t k = i+(Nkx)*j;
863 fft[k][0] = fphi[k][0];
864 fft[k][1] = fphi[k][1];
865 }
866 }
867
868 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
869
870 phi_bar.resize(NN);
871 for( size_t i=0; i<NN; i++ ) phi_bar[i] = float(realsp[i]/NN);
872 }
873 }
874 // gamma1
875 {
876
877 // build modes for each pixel in the fourier space
878 for( int i=0; i<Nkx; i++ ){
879 for( int j=0; j<ny; j++ ){
880
881 size_t k = i+(Nkx)*j;
882 double tmp = 0.5*(kxs[i]*kxs[i] - kys[j]*kys[j]);
883 // gamma
884 fft[k][0] = tmp * fphi[k][0];
885 fft[k][1] = tmp * fphi[k][1];
886 }
887 }
888
889 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
890
891 gamma1_bar.resize(NN);
892 for(size_t i=0; i<NN; i++ ) gamma1_bar[i] = float( realsp[i]/NN);
893 }
894 // gamma2
895 {
896
897 // build modes for each pixel in the fourier space
898 for( int i=0; i<Nkx; i++ ){
899 for( int j=0; j<ny; j++ ){
900
901 size_t k = i+(Nkx)*j;
902 double tmp = kxs[i]*kys[j];
903 // gamma
904 fft[k][0] = tmp * fphi[k][0];
905 fft[k][1] = tmp * fphi[k][1];
906
907 }
908 }
909
910 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
911
912 gamma2_bar.resize(NN);
913 for( size_t i=0; i<NN; i++ ) gamma2_bar[i] = float(-realsp[i]/NN);
914 }
915
916 // kappa - this is done over because of the window in Fourier space
917 {
918
919 // build modes for each pixel in the fourier space
920 for( int i=0; i<Nkx; i++ ){
921 for( int j=0; j<ny; j++ ){
922
923 double k2 = -(kxs[i]*kxs[i] + kys[j]*kys[j])/2;
924
925 size_t k = i+(Nkx)*j;
926
927 // surface density
928 fft[k][0] = k2*fphi[k][0];
929 fft[k][1] = k2*fphi[k][1];
930 }
931 }
932
933 fftw_execute_dft_c2r(plan.plan_c2r,fft,realsp.data());
934
935 for( size_t i=0; i<NN; i++ ) surface_density[i] = (realsp[i]/NN);
936 }
937
938 // std:: cout << " remapping the map in the original size " << std:: endl;
939 delete[] fft;
940 delete[] fphi;
941}
942
943#endif
944/* MultiLENS_H_ */
945
946
947
COSMOLOGY
The cosmology and all the functions required to calculated quantities based on the cosmology.
Definition cosmo.h:52
CPFITS_READ
read only fits file interface
Definition cpfits.h:232
LensHalo
A base class for all types of lensing "halos" which are any mass distribution that cause lensing.
Definition lens_halos.h:56
LensHalo::getDist
PosType getDist() const
Definition lens_halos.h:118
LensHalo::operator=
LensHalo & operator=(const LensHalo &h)
Definition lens_halos.cpp:169
LensHaloMultiMap
A lens halo that calculates all lensing quantities on two grids - a low res long range grid and a hig...
Definition multimap.h:177
LensHaloMultiMap::getNy_sr
size_t getNy_sr() const
return number of pixels on a y-axis side in short range map (margins are hidden internally)
Definition multimap.h:276
LensHaloMultiMap::getCenter_lr
Point_2d getCenter_lr() const
center of short range map in physical Mpc
Definition multimap.h:251
LensHaloMultiMap::getNx_sr
size_t getNx_sr() const
return number of pixels on a x-axis side in short range map (margins are hidden internally)
Definition multimap.h:274
LensHaloMultiMap::plan_short_range
static my_fftw_plan plan_short_range
these plans will be shared between all instances of LensHaloMultiMap
Definition multimap.h:205
LensHaloMultiMap::resetsubmap
void resetsubmap(int i, const std::vector< long > &lower_left_pix)
Definition multimap.cpp:400
LensHaloMultiMap::getLowerLeft_sr
Point_2d getLowerLeft_sr(int i) const
lower left of short range map in physical Mpc
Definition multimap.h:254
LensHaloMultiMap::resetsubmapPhys
void resetsubmapPhys(int i, Point_2d ll)
Sets the least highres smaller map in physical coordinates relative to the center of the original map...
Definition multimap.cpp:386
LensHaloMultiMap::NumberOfShortRangeMaps
int NumberOfShortRangeMaps()
the number of short range grids
Definition multimap.h:261
LensHaloMultiMap::LensHaloMultiMap
LensHaloMultiMap(std::string fitsfile, std::string dir_data, double redshift, double mass_unit, double npane, int number_of_subfields, COSMOLOGY &c, bool write_subfields=false, std::string dir_scratch="", bool subtract_ave=true, double ffactor=5, double gfactor=5, double rsmooth2=0)
Definition multimap.cpp:26
LensHaloMultiMap::getNy_lr
size_t getNy_lr() const
return number of pixels on a y-axis side in original map
Definition multimap.h:271
LensHaloMultiMap::getNx_lr
size_t getNx_lr() const
return number of pixels on a x-axis side in original map
Definition multimap.h:269
LensHaloMultiMap::getLowerLeft_lr
Point_2d getLowerLeft_lr() const
lower left of long range map in physical Mpc
Definition multimap.h:247
LensHaloMultiMap::getCenter_sr
Point_2d getCenter_sr(int i) const
center of short range map in physical Mpc
Definition multimap.h:258
LensHaloMultiMap::getRangeMpc_sr
double getRangeMpc_sr() const
return range of long range map in physical Mpc
Definition multimap.h:266
LensHaloMultiMap::getUpperRight_lr
Point_2d getUpperRight_lr() const
upper right of short range map in physical Mpc
Definition multimap.h:249
LensHaloMultiMap::getNx
size_t getNx() const
return number of pixels on a x-axis side in original map
Definition multimap.h:279
LensHaloMultiMap::getUpperRight_sr
Point_2d getUpperRight_sr(int i) const
upper right of short range map in physical Mpc
Definition multimap.h:256
LensHaloMultiMap::getNy
size_t getNy() const
return number of pixels on a y-axis side in original map
Definition multimap.h:281
LensHaloMultiMap::getRangeMpc_lr
double getRangeMpc_lr() const
return range of long range map in physical Mpc
Definition multimap.h:264
LensingVariable
LensingVariable
output lensing variables
Definition standard.h:89
LensMap
The MOKA map structure, containing all quantities that define it.
Definition multimap.h:56
LensMap::read_sub
void read_sub(CPFITS_READ &cpfits, std::vector< long > &first, std::vector< long > &last, double Dist)
read a subsection of the fits map
Definition multimap.cpp:917
LensMap::ProcessFFTs
void ProcessFFTs(float zerosize, T Wphi_of_k, my_fftw_plan &plan_padded, bool do_alpha=true)
a thread safe version of ProcessFFTs
Definition multimap.h:430
LensMap::read
void read(std::string input_fits, double angDist)
read an entire map
Definition multimap.cpp:747
LensMap::Myread
void Myread(std::string fits_input_file)
read from a file that has been generated with LensMap::write()
Definition multimap.cpp:814
LensMap::surface_density
std::valarray< double > surface_density
values for the map
Definition multimap.h:65
LensMap::read_header
void read_header(std::string input_fits, double angDist)
read only header information
Definition multimap.cpp:711
LensMap::upperright
Point_2d upperright
boundery with centred grid
Definition multimap.h:77
LensMap::x_resolution
double x_resolution()
resolution in Mpc
Definition multimap.h:80
LensMap::write
void write(std::string filename)
write the fits file of the maps of all the lensing quantities.
Definition multimap.cpp:977
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition point.h:48
Generated on Sun Sep 8 2024 19:31:02 for GLAMERDOC++ by doxygen 1.12.0