lagSphericTransform.cc 7.11 KB
Newer Older
1 2
#include <stdlib.h>
#include <iostream>
3 4
#include <algorithm> //copy, sort
#include <iterator> // for ostream_iterator
5 6 7

#include "lagSphericTransform.h"
#include "lagsht_exceptions.h"
8

9 10
#include "walltimer.h"    //timing

11

12

13 14 15 16
#define DEBUG 0

namespace LagSHT {

17

18 19 20
  void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& flmn, 
					     std::vector<r_8>& fijk,
					     std::vector< std::complex<r_8> >& flmk) {
21 22
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
23
  if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
24 25 26 27
    
  fijk.resize(Ntot3DPix);

  //Perform the Inverse Laguerre Transform first
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
28 29

  flmk.resize(NtotCoeff);
30 31 32 33 34

#if DEBUG >= 1
  cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
  cout << "Start   Inverse Laguerre...." << endl;
#endif
35 36
  
  tstack_push("Laguerre MultiSynthesis");
37
  lagTrans_->MultiSynthesis(flmn, flmk, Nalm_);
38
  tstack_pop("Laguerre MultiSynthesis");
39 40 41 42 43 44 45


#if DEBUG >= 1
  cout << "End   Inverse Laguerre...." << endl;
  cout << "Start Inverse H.Spherical Transform " << endl;
#endif

46
  tstack_push("SHT Synthesis");
47 48 49 50 51 52 53
  int off7n, off7k;
  //perform the synthesis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 2
54
    std::cout << "Synthesis... at shell num: " << n << std::endl;
55
#endif
56
    sphere_->alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
57 58

  }//loop on shell
59 60 61

    tstack_pop("SHT Synthesis");

62
#if DEBUG >= 1
63
    std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
64 65 66 67 68
#endif
}//Synthesis



69 70 71
  void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk, 
					    std::vector< std::complex<r_8> >& flmn,
					    std::vector< std::complex<r_8> >& almk
72
					  ) {
73 74
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
75 76 77 78
  if(fijk.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error: fijk");
  if(almk.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Analysis size error: almk");

  //done in Laguerre routine  flmn.resize(NtotCoeff);
79 80

#if DEBUG >= 1
81
  std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
82 83 84 85
#endif
  //Perform the Spherical Hamonic analysis first
    
  int off7n, off7k;
86 87 88

  tstack_push("SHT Analysis");

89 90 91 92 93 94
  //Ylm analysis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 1
95
    std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
96 97
#endif

98
    sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
99
    
100 101 102
  }//loop on shell

#if DEBUG >= 1
103 104
  std::cout << "End   Harmonic Spheric..." << std::endl;
  std::cout << "Start Laguerre..." << std::endl;
105 106
#endif

107 108
  tstack_pop("SHT Analysis");
  tstack_push("Laguerre MultiAnalysis");
109 110

  //Laguerre Transform
111
  lagTrans_->MultiAnalysis(almk, flmn, Nalm_);
112 113

#if DEBUG >= 1
114
  std::cout << "End Laguerre..." << std::endl;
115 116
#endif

117 118 119
  tstack_pop("Laguerre MultiAnalysis");


120 121
#if DEBUG >= 2
  for (int i=0;i<NtotCoeff; i++) {
122
    std::cout << "flmn("<<i<<"): " <<flmn[i] <<  std::endl;
123 124 125
  }
#endif
#if DEBUG >= 1
126
  std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
127 128 129 130 131
#endif
}//Analysis



132 133 134 135 136 137
  void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& Elmn, 
					     const std::vector< std::complex<r_8> >& Blmn,
					     int spin,
					     std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im,
					     std::vector< std::complex<r_8> >& Elmk, 
					     std::vector< std::complex<r_8> >& Blmk
138 139 140 141 142 143 144 145 146 147 148 149 150 151
					   )   {
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(Elmn.size() != (size_t)NtotCoeff || Blmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
    
  fijk_re.resize(Ntot3DPix);
  fijk_im.resize(Ntot3DPix);

  //Perform the Inverse Laguerre Transform first

  Elmk.resize(NtotCoeff);
  Blmk.resize(NtotCoeff);

#if DEBUG >= 1
152 153
  std::cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << std::endl;
  std::cout << "Start   Inverse Laguerre...." << std::endl;
154 155 156 157 158 159 160 161
#endif
  
  tstack_push("Laguerre MultiSynthesis");
  lagTrans_->MultiSynthesis(Elmn, Elmk, Nalm_);
  lagTrans_->MultiSynthesis(Blmn, Blmk, Nalm_);
  tstack_pop("Laguerre MultiSynthesis");

#if DEBUG >= 1
162 163
  std::cout << "End   Inverse Laguerre...." << std::endl;
  std::cout << "Start Inverse H.Spherical Transform " << std::endl;
164 165 166 167 168 169 170 171 172 173
#endif

  tstack_push("SHT Synthesis");
  int off7n, off7k;
  //perform the synthesis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 2
174
    std::cout << "Synthesis... at shell num: " << n << std::endl;
175 176 177 178 179 180 181 182 183 184 185 186 187
#endif
    sphere_->alm2map_spin(reinterpret_cast<const r_8*>(&(Elmk.data()[off7n])),
			  reinterpret_cast<const r_8*>(&(Blmk.data()[off7n])),
			  &fijk_re.data()[off7k],
			  &fijk_im.data()[off7k],
			  spin,
			  false); //the false is to forbidd accumulation

  }//loop on shell

    tstack_pop("SHT Synthesis");

#if DEBUG >= 1
188
    std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
189 190 191 192
#endif
}//Synthesis


193 194 195 196 197 198 199 200
  void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk_re, 
					    const std::vector<r_8>& fijk_im,
					    int spin,
					    std::vector< std::complex<r_8> >& Elmn, 
					    std::vector< std::complex<r_8> >& Blmn,
					    std::vector< std::complex<r_8> >& Elmk, 
					    std::vector< std::complex<r_8> >& Blmk
					    ) {
201 202 203
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(fijk_re.size() != (size_t)Ntot3DPix || fijk_im.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
204 205
//done in Laguerre routine   Elmn.resize(NtotCoeff);
//   Blmn.resize(NtotCoeff);
206 207 208 209 210 211 212 213

  //Perform the Inverse Laguerre Transform first

  Elmk.resize(NtotCoeff);
  Blmk.resize(NtotCoeff);


#if DEBUG >= 1
214
  std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
215 216 217 218 219 220 221 222 223 224 225 226 227
#endif
  //Perform the Spherical Hamonic analysis first
    
  int off7n, off7k;

  tstack_push("SHT Analysis");

  //Ylm analysis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 1
228
    std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
229 230 231 232 233 234 235 236 237 238 239 240 241
#endif

    sphere_->map2alm_spin(&(fijk_re.data()[off7k]), 
			  &(fijk_im.data()[off7k]),
			  reinterpret_cast<r_8*>(&(Elmk.data()[off7n])),
			  reinterpret_cast<r_8*>(&(Blmk.data()[off7n])),
			  spin,
			  false); //the false is to forbidd accumulation


  }//loop on shell

#if DEBUG >= 1
242 243
  std::cout << "End   Harmonic Spheric..." << std::endl;
  std::cout << "Start Laguerre..." << std::endl;
244 245 246 247 248 249 250 251 252 253
#endif

  tstack_pop("SHT Analysis");
  tstack_push("Laguerre MultiAnalysis");

  //Laguerre Transform
  lagTrans_->MultiAnalysis(Elmk, Elmn, Nalm_);
  lagTrans_->MultiAnalysis(Blmk, Blmn, Nalm_);

#if DEBUG >= 1
254
  std::cout << "End Laguerre..." << std::endl;
255 256 257 258 259 260
#endif

  tstack_pop("Laguerre MultiAnalysis");

#if DEBUG >= 2
  for (int i=0;i<NtotCoeff; i++) {
261
    std::cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] <<  std::endl;
262 263 264
  }
#endif
#if DEBUG >= 1
265
  std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
266 267 268
#endif
}//Analysis

269
}// Fin de namespace