![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
Changeset 3510 in Sophya for trunk/SophyaLib/Samba
- Timestamp:
- Aug 8, 2008, 3:11:45 PM (17 years ago)
- Location:
- trunk/SophyaLib/Samba
- Files:
-
- 4 edited
-
alm.cc (modified) (6 diffs)
-
alm.h (modified) (2 diffs)
-
sphericaltransformserver.cc (modified) (4 diffs)
-
sphericaltransformserver.h (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/Samba/alm.cc
r2885 r3510 1 1 #include "sopnamsp.h" 2 2 #include "alm.h" 3 4 5 /*! 6 \class Alm 7 \ingroup Samba 8 Class for manipulating the coefficients \f$a_{lm}\f$ of the development 9 in spherical harmonics of a function efined on a sphere. 10 */ 11 12 /*! 13 fwhm specifies the gaussian beam half witdh in arc.minutes 14 */ 3 15 template <class T> 4 16 Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm) 5 17 6 18 { 19 int_4 nlmax= clin.NElts()-1; 20 21 //alm.ReSize(nlmax); 22 this->ReSizeRow(nlmax+1); 23 RandomGenerator rg(1, false); 24 GenFromCl(clin, fwhm, rg); 25 } 26 27 /*! 28 fwhm specifies the gaussian beam half witdh in arc.minutes 29 */ 30 template <class T> 31 Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm, RandomGenerator & rg) 32 { 33 int_4 nlmax= clin.NElts()-1; 34 35 //alm.ReSize(nlmax); 36 this->ReSizeRow(nlmax+1); 37 GenFromCl(clin, fwhm, rg); 38 } 7 39 8 40 41 42 template <class T> 43 void Alm<T>::GenFromCl(const TVector<T> & clin, const r_8 fwhm, RandomGenerator & rg) 44 { 9 45 /*======================================================================= 10 46 creates the a_lm from the power spectrum, … … 21 57 int_4 nlmax= clin.NElts()-1; 22 58 23 //alm.ReSize(nlmax);24 this->ReSizeRow(nlmax+1);25 26 59 r_8 sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI; 27 60 int_4 n_l = nlmax+1; … … 29 62 30 63 // --- smoothes the initial power spectrum --- 31 TVector<T> cl =clin;64 TVector<T> cl(clin, false); 32 65 int l; 33 66 for (l=0;l<n_l;l++) … … 45 78 T rms=sqrt(cl(l)); 46 79 // ------ m = 0 ------ 47 complex<T> zeta1( NorRand());80 complex<T> zeta1((T)rg.Gaussian() ); 48 81 (*this)(l,0) = zeta1 * rms; 49 82 … … 52 85 { 53 86 complex<T> aux1(hsqrt2); 54 complex<T> aux2( NorRand(),NorRand());87 complex<T> aux2((T)rg.Gaussian() , (T)rg.Gaussian() ); 55 88 zeta1=aux1*aux2; 56 89 (*this)(l,m)=rms*zeta1; … … 79 112 return powsp; 80 113 } 114 115 /*! 116 \class Bm 117 \ingroup Samba 118 Class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ 119 (then the size of the vector will be actually \f$2m_{max}+1)\f$. 120 This class is used by the spherical harmonics transform server 121 */ 122 81 123 82 124 #ifdef __CXX_PRAGMA_TEMPLATES__ -
trunk/SophyaLib/Samba/alm.h
r3077 r3510 6 6 #define ALM_SEEN 7 7 8 #include "s randgen.h"8 #include "stsrand.h" 9 9 #include "nbmath.h" 10 10 #include "triangmtx.h" … … 13 13 namespace SOPHYA { 14 14 15 / *! class for the coefficients \f$a_{lm}\f$ of the development of a function defined on a sphere, in spherical harmonics */15 // ----- Classe Alm 16 16 template <class T> 17 class Alm : public TriangularMatrix<complex<T> > 18 { 19 public : 17 class Alm : public TriangularMatrix< complex<T> > { 18 public : 20 19 21 Alm() : TriangularMatrix<complex<T> >() {;}; 22 /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */ 23 Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;} 24 Alm(const Alm<T>& a, bool share=false) : TriangularMatrix<complex<T> >(a, share) {;} 20 Alm() : TriangularMatrix<complex<T> >() {} 21 //! instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s 22 Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {} 23 //! copy constructor 24 Alm(const Alm<T>& a, bool share=false) : TriangularMatrix<complex<T> >(a, share) {} 25 //! Constructor, generates alm from a cl vector 26 Alm(const TVector<T> & clin, const r_8 fwhm) ; 27 //! Constructor, generates alm from a cl vector, using the specified random generator object 28 Alm(const TVector<T> & clin, const r_8 fwhm, RandomGenerator & rg); 25 29 26 Alm(const TVector<T>& clin, const r_8 fwhm) ; 27 /*! resize with a new lmax */ 28 inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);} 29 inline sa_size_t Lmax() const {return this->rowNumber()-1;} 30 TVector<T> powerSpectrum() const; 30 //! Resize with a new lmax 31 inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);} 32 inline sa_size_t Lmax() const {return this->rowNumber()-1;} 31 33 32 inline Alm<T>& SetComplexT(complex<T> a) 33 { 34 this->SetT(a); 35 return *this; 36 } 37 inline Alm<T>& operator = (complex<T> a) 38 { 39 return SetComplexT(a); 40 } 34 //! Computes the corresponding pwer spectrum (Cl vector) 35 TVector<T> powerSpectrum() const; 41 36 37 inline Alm<T>& SetComplexT(complex<T> a) 38 { this->SetT(a); return *this; } 39 inline Alm<T>& operator = (complex<T> a) 40 { return SetComplexT(a); } 41 42 private: 43 void GenFromCl(const TVector<T> & clin, const r_8 fwhm, RandomGenerator& rg); 42 44 }; 43 /*! class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ (then the size of the vector will be actually \f$2m_{max}+1)\f$ */ 45 46 47 //------------ Classe Bm , vecteur avec index -mMax <= index <= +mMax 44 48 template <class T> 45 class Bm 46 { 47 public : 48 Bm(): nmmax_(0) {;}; 49 /* instanciate from the maximum absolute value of the index */ 50 Bm(sa_size_t mmax) : nmmax_(mmax) { bm_.ReSize( (2*mmax+1) );} 51 Bm(const Bm<T>& b, bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;} 52 /*! resize with a new value of mmax */ 53 inline void ReSizeToMmax(sa_size_t mmax) 54 { 55 nmmax_= mmax; 56 bm_.ReSize(2* nmmax_+1); 57 } 58 //inline sa_size_t Size() const {return 2*nmmax_+1;}; 59 inline T& operator()(sa_size_t m) {return bm_( adr_i(m));}; 60 inline T const& operator()(sa_size_t m) const {return *(bm_.Begin()+ adr_i(m));}; 61 /*! return the current value of the maximum absolute value of the index */ 62 inline sa_size_t Mmax() const 63 { 64 return nmmax_; 65 } 66 inline Bm<T>& operator = (const Bm<T>& b) 67 { 68 if (this != &b) 69 { 49 class Bm { 50 public : 51 Bm(): nmmax_(0) {;}; 52 //! Constructor with the specification of the maximum absolute value of the index : |mMax| 53 //! Copy contrsuctor 54 Bm(sa_size_t mmax) : nmmax_(mmax) { bm_.ReSize( (2*mmax+1) );} 55 Bm(const Bm<T>& b, bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) { } 56 //! resize with a new value of mmax 57 inline void ReSizeToMmax(sa_size_t mmax) 58 { nmmax_= mmax; bm_.ReSize(2* nmmax_+1); } 59 60 //! Acces to element with index m ( -mMax <= m <= +mMax ) 61 inline T operator()(sa_size_t m) const { return bm_( m+nmmax_); } 62 //! Acces to element with index m ( -mMax <= m <= +mMax ) 63 inline T& operator()(sa_size_t m) { return bm_( m+nmmax_); } 64 65 //! Return the current value of the maximum absolute value of the index 66 inline sa_size_t Mmax() const { return nmmax_; } 67 68 inline Bm<T>& operator = (const Bm<T>& b) 69 { 70 if (this != &b) { 70 71 nmmax_= b.nmmax_; 71 72 bm_= b.bm_; 72 } 73 return *this; 74 } 75 76 private: 77 /*! return the address in the array representing the vector */ 78 inline sa_size_t adr_i(sa_size_t i) const 79 { 80 return (i+nmmax_); 81 } 73 } 74 return *this; 75 } 82 76 83 77 84 85 sa_size_t nmmax_; 86 NDataBlock<T> bm_; 87 78 private: 79 sa_size_t nmmax_; 80 NDataBlock<T> bm_; 88 81 }; 89 82 -
trunk/SophyaLib/Samba/sphericaltransformserver.cc
r3508 r3510 6 6 #include "sphericaltransformserver.h" 7 7 #include "tvector.h" 8 #include "s randgen.h"8 #include "stsrand.h" 9 9 #include "nbmath.h" 10 10 #include "timing.h" … … 152 152 153 153 */ 154 155 //! Default constructor - Creates a non thread-safe RandomGenerator to be used by GenerateFromCl 156 template<class T> 157 SphericalTransformServer<T>::SphericalTransformServer() 158 : rg_(1, false) 159 { 160 fftIntfPtr_=new FFTPackServer(true); // preserveinput = true 161 fftIntfPtr_->setNormalize(false); 162 } 163 164 //! Constructor with the specification of a RandomGenerator object to be used by GenerateFromCl 165 template<class T> 166 SphericalTransformServer<T>::SphericalTransformServer(RandomGenerator const & rg) 167 : rg_(rg) 168 { 169 fftIntfPtr_=new FFTPackServer(true); // preserveinput = true 170 fftIntfPtr_->setNormalize(false); 171 } 172 173 template<class T> 174 SphericalTransformServer<T>::~SphericalTransformServer() 175 { 176 if (fftIntfPtr_!=NULL) delete fftIntfPtr_; 177 } 178 179 /*! 180 Set a fft server. The constructor sets a default fft server (fft-pack). 181 So it is not necessary to call this method for a standard use. 182 \warning The FFTServerInterface object should NOT overwrite the input arrays 183 */ 184 template<class T> 185 void SphericalTransformServer<T>::SetFFTServer(FFTServerInterface* srv) 186 { 187 if (fftIntfPtr_!=NULL) delete fftIntfPtr_; 188 fftIntfPtr_=srv; 189 fftIntfPtr_->setNormalize(false); 190 } 191 154 192 155 193 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const … … 1275 1313 // Alm<T> a2lme = almFromCl(Cle, fwhm); 1276 1314 // Alm<T> a2lmb = almFromCl(Clb, fwhm); 1277 Alm<T> a2lme(Cle, fwhm );1278 Alm<T> a2lmb(Clb, fwhm );1315 Alm<T> a2lme(Cle, fwhm, rg_); 1316 Alm<T> a2lmb(Clb, fwhm, rg_); 1279 1317 1280 1318 GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb); … … 1293 1331 { 1294 1332 1295 Alm<T> alm(Cl, fwhm );1333 Alm<T> alm(Cl, fwhm, rg_); 1296 1334 GenerateFromAlm(sph,pixelSizeIndex, alm ); 1297 1335 } -
trunk/SophyaLib/Samba/sphericaltransformserver.h
r3003 r3510 5 5 #include "fftservintf.h" 6 6 #include "fftpserver.h" 7 #include "stsrand.h" 7 8 #include "alm.h" 8 9 #include "lambdaBuilder.h" … … 18 19 public: 19 20 20 SphericalTransformServer() 21 { 22 fftIntfPtr_=new FFTPackServer(true); // preserveinput = true 23 fftIntfPtr_->setNormalize(false); 24 }; 25 ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;}; 21 SphericalTransformServer(); 22 SphericalTransformServer(RandomGenerator const & rg); 23 virtual ~SphericalTransformServer(); 24 void SetFFTServer(FFTServerInterface* srv=NULL); 26 25 27 /*!28 Set a fft server. The constructor sets a default fft server (fft-pack). So it is not necessary to call this method for a standard use.29 \warning The FFTServerInterface object should NOT overwrite the input arrays30 */31 void SetFFTServer(FFTServerInterface* srv=NULL)32 {33 if (fftIntfPtr_!=NULL) delete fftIntfPtr_;34 fftIntfPtr_=srv;35 fftIntfPtr_->setNormalize(false);36 }37 26 38 27 void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; … … 44 33 int_4 pixelSizeIndex, 45 34 const TVector<T>& Cle, const TVector<T>& Clb, 46 const r_8 fwhm) const;35 const r_8 fwhm) const; 47 36 48 37 … … 118 107 119 108 FFTServerInterface* fftIntfPtr_; 109 mutable RandomGenerator rg_; 120 110 121 111 };
Note:
See TracChangeset
for help on using the changeset viewer.
