lagSphericTransform.cc 3.04 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
//#include "lagsht_healpixhelper.h"
9

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

12

13

14 15 16 17
#define DEBUG 0

namespace LagSHT {

18

19 20 21 22
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn, 
					   vector<r_8>& fijk) {
  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 28 29 30 31 32 33
    
  fijk.resize(Ntot3DPix);

  //Perform the Inverse Laguerre Transform first
  vector< complex<r_8> > flmk(NtotCoeff); 

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


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

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

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

  }//loop on shell
58 59 60

    tstack_pop("SHT Synthesis");

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



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

#if DEBUG >= 1
  cout << "SphLagTransform::Analysis.... Go..." << endl;
#endif
  //Perform the Spherical Hamonic analysis first
    
82
  //  vector< complex<r_8> > flmk(NtotCoeff);
83 84

  int off7n, off7k;
85 86 87

  tstack_push("SHT Analysis");

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

#if DEBUG >= 1
    cout << "Start Harmonic Spheric... at n = "<< n << endl;
#endif

97
    sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
98 99


100 101 102 103 104 105 106
  }//loop on shell

#if DEBUG >= 1
  cout << "End   Harmonic Spheric..." << endl;
  cout << "Start Laguerre..." << endl;
#endif

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

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

#if DEBUG >= 1
  cout << "End Laguerre..." << endl;
#endif

117 118 119
  tstack_pop("Laguerre MultiAnalysis");


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



}// Fin de namespace