GLAMERDOC++
Gravitational Lensing Code Library
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
19 struct 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 
56 struct 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 
176 class LensHaloMultiMap : public LensHalo
177 {
178 public:
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 
351 public:
352  LensMap long_range_map;
353  std::vector<LensMap> short_range_maps;
354 
355 private:
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 
416 template <typename T>
417 void 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 
429 template <typename T>
430 void 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 
726 template <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 
735 template <class T>
736 void 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 
LensHaloMultiMap::getNx_lr
size_t getNx_lr() const
return number of pixels on a x-axis side in original map
Definition: multimap.h:269
LensHalo::operator=
LensHalo & operator=(const LensHalo &h)
Definition: lens_halos.cpp:169
LensHaloMultiMap::getRangeMpc_lr
double getRangeMpc_lr() const
return range of long range map in physical Mpc
Definition: multimap.h:264
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::getUpperRight_sr
Point_2d getUpperRight_sr(int i) const
upper right of short range map in physical Mpc
Definition: multimap.h:256
LensMap::x_resolution
double x_resolution()
resolution in Mpc
Definition: multimap.h:80
LensHaloMultiMap::plan_short_range
static my_fftw_plan plan_short_range
these plans will be shared between all instances of LensHaloMultiMap
Definition: multimap.h:201
CPFITS_READ
read only fits file interface
Definition: cpfits.h:232
LensHaloMultiMap::getNy_lr
size_t getNy_lr() const
return number of pixels on a y-axis side in original map
Definition: multimap.h:271
LensMap::read_header
void read_header(std::string input_fits, double angDist)
read only header information
Definition: multimap.cpp:711
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
LensHalo
A base class for all types of lensing "halos" which are any mass distribution that cause lensing.
Definition: lens_halos.h:57
LensHaloMultiMap::getLowerLeft_sr
Point_2d getLowerLeft_sr(int i) const
lower left of short range map in physical Mpc
Definition: multimap.h:254
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition: point.h:48
LensHaloMultiMap::getNx
size_t getNx() const
return number of pixels on a x-axis side in original map
Definition: multimap.h:279
LensHaloMultiMap::getUpperRight_lr
Point_2d getUpperRight_lr() const
upper right of short range map in physical Mpc
Definition: multimap.h:249
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::upperright
Point_2d upperright
boundery with centred grid
Definition: multimap.h:77
LensHaloMultiMap::getNy
size_t getNy() const
return number of pixels on a y-axis side in original map
Definition: multimap.h:281
LensHaloMultiMap::getLowerLeft_lr
Point_2d getLowerLeft_lr() const
lower left of long range map in physical Mpc
Definition: multimap.h:247
standard.h
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
LensMap
The MOKA map structure, containing all quantities that define it.
Definition: multimap.h:56
LensHaloMultiMap::resetsubmap
void resetsubmap(int i, const std::vector< long > &lower_left_pix)
Definition: multimap.cpp:400
LensHaloMultiMap::getCenter_sr
Point_2d getCenter_sr(int i) const
center of short range map in physical Mpc
Definition: multimap.h:258
LensHaloMultiMap::getCenter_lr
Point_2d getCenter_lr() const
center of short range map in physical Mpc
Definition: multimap.h:251
LensMap::write
void write(std::string filename)
write the fits file of the maps of all the lensing quantities.
Definition: multimap.cpp:977
LensMap::surface_density
std::valarray< double > surface_density
values for the map
Definition: multimap.h:65
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::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
LensHalo::getDist
PosType getDist() const
Definition: lens_halos.h:122
InputParams.h
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
LensingVariable
LensingVariable
output lensing variables
Definition: standard.h:89
LensHaloMultiMap::getRangeMpc_sr
double getRangeMpc_sr() const
return range of long range map in physical Mpc
Definition: multimap.h:266
LensMap::read
void read(std::string input_fits, double angDist)
read an entire map
Definition: multimap.cpp:747
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
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18