Changeset 1218 in Sophya
- Timestamp:
- Oct 3, 2000, 2:14:14 PM (25 years ago)
- Location:
- trunk
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaExt/FitsIOServer/fitsfile.cc
r1209 r1218 8 8 9 9 10 void BnTblLine::setFormat(int dc, int fc, int ic, int cc, vector<string> names) 11 { 12 int nbcols = dc + fc + ic + cc; 13 int maxName = names.size(); 14 if (nbcols != maxName) 15 { 16 cout << " WARNING: BnTblLine:: length of vector of column names not equal to total number of columns" << endl; 17 maxName = nbcols < maxName ? nbcols : maxName; 18 } 19 ColName_ = vector<string>(nbcols); 20 for (int k=0; k < maxName; k++) ColName_[k] = names[k]; 21 if (dc >0) ddata_ = vector<double>(dc); 22 if (fc >0) fdata_ = vector<float>(fc); 23 if (ic >0) idata_ = vector<int>(fc); 24 if (cc >0) cdata_ = vector<string>(fc); 25 } 26 27 bool BnTblLine::sameFormat(const BnTblLine& btl) const 28 { 29 if (btl.ddata_.size() == ddata_.size() && btl.fdata_.size() == fdata_.size() && btl.idata_.size() == idata_.size() && btl.cdata_.size() == cdata_.size()) return true; 30 else return false; 31 } 32 33 void BnTblLine::Print() 34 { 35 int k; 36 cout << " ********* ligne ************* " << endl; 37 cout << " *** noms de variables " << endl; 38 for (k=0; k < ColName_.size(); k++) cout << ColName_[k] << " "; 39 cout << endl; 40 cout << " *** variables doubles " << endl; 41 for (k=0; k < ddata_.size(); k++) cout << ddata_[k] << " "; 42 cout << endl; 43 cout << " *** variables float " << endl; 44 for (k=0; k < fdata_.size(); k++) cout << fdata_[k] << " "; 45 cout << endl; 46 cout << " *** variables int " << endl; 47 for (k=0; k < idata_.size(); k++) cout << idata_[k] << " "; 48 cout << endl; 49 cout << " *** variables string " << endl; 50 for (k=0; k < cdata_.size(); k++) cout << cdata_[k] << " "; 51 cout << endl; 52 cout << " ***************************** " << endl; 53 } 54 55 56 57 /*! 58 \class SOPHYA::FitsIOHandler 59 The class structure is analogous to Sophya-PPersist system : 60 Each SOPHYA object XXX is associated with a object of class FITS_XXX 61 (inheriting from FitsFileHandler), to which input/output operations with FITS 62 files are delegated (through a class Hierarchy : FitsFile (virtual), 63 FitsInFile, FitsOutFile) . A typical example of use is the following : 64 65 \verbatim 66 int m=... ; 67 SphereHEALPix<r_8> sphere1(m); // definition of the SOPHYA object 68 .... fill the sphere .... 69 70 FITS_SphereHEALPix<r_8> fits_sph1(sphere1); 71 // delegated object 72 fits_sph.Write("myfile.fits"); // writing on FITS file 73 74 FITS_SphereHEALPix<r_8> fits_sph2("myfile.fits"); 75 // load a delegated object 76 // from FITS file 77 SphereHEALPix<r_8> sphere2=(SphereHEALPix<r_8>)fits_sph2; 78 // casting the delegated object 79 // into a SOPHYA object 80 \endverbatim 81 82 83 */ 84 85 /*! \fn void SOPHYA::FitsIOHandler::Read(char flnm[],int hdunum) 86 87 this method is called from inherited objects : 88 89 opens a file 'flnm' 90 91 gets parameters in extension-header (hdunum) 92 93 calls the method 'ReadFromFits' from the inherited object 94 */ 10 95 void FitsIOHandler::Read(char flnm[],int hdunum) 11 96 { … … 13 98 Read(ifts, hdunum); 14 99 } 100 101 /*! \fn void SOPHYA::FitsIOHandler::Read(FitsInFile& is, int hdunum) 102 Read the data on extension hdunum (or primary header, if hdunum=1) from FitsInFIle. If hdunum is not addressed, , one reads the next extension, with respect to the current position. 103 */ 15 104 void FitsIOHandler::Read(FitsInFile& is, int hdunum) 16 105 { … … 21 110 22 111 112 /*! \fn void SOPHYA::FitsIOHandler::Write(char flnm[]) 113 this method is called from inherited objects. 114 115 for writing a new object in a new fits-extension : 116 117 \warning By convention, primary header does not contain fits-image data : i.e. 118 all data are fits-extensions. The first relevant header will have hdunum=2. 119 For switching off this convention use the method : 120 121 firstImageOnPrimaryHeader() (see below) 122 123 In that case do not forget to precise hdunum=1 when reading data on primary header. 124 125 calls the method 'WriteToFits' from the inherited object 126 127 */ 23 128 void FitsIOHandler::Write(char flnm[]) 24 129 … … 33 138 } 34 139 140 141 /*! 142 \class SOPHYA::FitsIOHandler 143 Class (virtual) for managing FITS format files 144 */ 35 145 36 146 … … 87 197 } 88 198 199 /*! 200 \class SOPHYA::FitsInFile 201 202 class for saving SOPHYA objects on FITS Format Files (uses cfitsio lib) 203 */ 204 89 205 FitsInFile::FitsInFile() 90 206 { 91 207 InitNull(); 92 208 } 93 //FitsInFile::FitsInFile(char flnm[], int hdunum) 209 94 210 FitsInFile::FitsInFile(char flnm[]) 95 211 { … … 118 234 } 119 235 236 ////////////////////////////////////////////////////////// 237 // methods with general purpose 238 ///////////////////////////////////////////////////////// 120 239 121 240 int FitsInFile::NbBlocks(char flnm[]) … … 200 319 } 201 320 321 322 void FitsInFile::ReadFInit(int hdunum) 323 { 324 InitNull(); 325 // int status = 0; 326 327 // fits_open_file(&fptr_,flnm,READONLY,&status); 328 // if( status ) printerror( status ); 329 330 int status = 0; 331 332 if (hdunum <= 1) 333 { 334 hdunum_ = 1; 335 // presence of image ? 336 int naxis= 0; 337 fits_read_key(fptr_,TINT,"NAXIS",&naxis,NULL,&status); 338 if( status ) printerror( status ); 339 if (naxis > 0 ) // there is an image 340 { 341 hdutype_ = IMAGE_HDU; 342 GetImageParameters (fptr_, bitpix_, naxis_, naxisn_); 343 nbData_ = 1; 344 int k; 345 for (k=0; k<naxis_; k++) if (naxisn_[k] > 0) nbData_ *= naxisn_[k]; 346 KeywordsIntoDVList(fptr_, dvl_,hdunum_); 347 } 348 else 349 { 350 throw PException(" first header : no image, probably error in hdunum"); 351 } 352 // 353 } 354 else 355 { 356 hdunum_ = hdunum-1; 357 int hdutype; 358 fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status); 359 if( status ) printerror( status,":FitsFile::ReadF : erreur movabs"); 360 moveToFollowingHeader(); 361 } 362 } 363 364 202 365 void FitsInFile::GetImageParameters (fitsfile* fileptr,int& bitpix,int& naxis,vector<int>& naxisn) 203 366 { … … 230 393 231 394 232 void FitsInFile::ReadFInit(int hdunum) 233 { 234 InitNull(); 235 // int status = 0; 236 237 // fits_open_file(&fptr_,flnm,READONLY,&status); 238 // if( status ) printerror( status ); 239 240 int status = 0; 241 242 if (hdunum <= 1) 243 { 244 hdunum_ = 1; 245 // presence of image ? 246 int naxis= 0; 247 fits_read_key(fptr_,TINT,"NAXIS",&naxis,NULL,&status); 248 if( status ) printerror( status ); 249 if (naxis > 0 ) // there is an image 250 { 251 hdutype_ = IMAGE_HDU; 252 GetImageParameters (fptr_, bitpix_, naxis_, naxisn_); 253 nbData_ = 1; 254 int k; 255 for (k=0; k<naxis_; k++) if (naxisn_[k] > 0) nbData_ *= naxisn_[k]; 256 KeywordsIntoDVList(fptr_, dvl_,hdunum_); 257 } 258 else 259 { 260 throw PException(" first header : no image, probably error in hdunum"); 261 } 262 // 263 } 264 else 265 { 266 hdunum_ = hdunum-1; 267 int hdutype; 268 fits_movabs_hdu(fptr_,hdunum_,&hdutype,&status); 269 if( status ) printerror( status,":FitsFile::ReadF : erreur movabs"); 270 moveToFollowingHeader(); 271 } 272 } 273 395 396 397 /*! \fn DVList SOPHYA::FitsInFile::DVListFromPrimaryHeader() const 398 399 \return the keywords of primary header in a DVList 400 401 */ 274 402 DVList FitsInFile::DVListFromPrimaryHeader() const 275 403 { … … 308 436 309 437 438 439 440 441 /*! \fn int SOPHYA::FitsInFile::NbColsFromFits() const 442 \return number of columns (return 1 if IMAGE) 443 */ 310 444 int FitsInFile::NbColsFromFits() const 311 445 { … … 320 454 } 321 455 456 /*! \fn int SOPHYA::FitsInFile::NentriesFromFits(int nocol) const 457 \return number of data in the current IMAGE extension on FITS file, or number 458 of data of column number 'nocol' of the current BINTABLE extension 459 */ 322 460 int FitsInFile::NentriesFromFits(int nocol) const 323 461 { … … 334 472 } 335 473 474 /*! \fn char SOPHYA::FitsInFile::ColTypeFromFits(int nocol) const 475 476 return a character denoting data type of column number 'nocol' in a BINTABLE : 477 478 D : double 479 480 E : float 481 482 I : integer 483 484 S : character string 485 486 */ 487 336 488 char FitsInFile::ColTypeFromFits(int nocol) const 337 489 { … … 342 494 return types_[nocol]; 343 495 } 496 497 498 /*! \fn string SOPHYA::FitsInFile::ColNameFromFits(int nocol) const 499 500 \return name of the column number 'nocol' of the current BINTABLE extension 501 */ 502 344 503 string FitsInFile::ColNameFromFits(int nocol) const 345 504 { … … 350 509 return noms_[nocol]; 351 510 } 511 512 /*! \fn int DSOPHYA::FitsInFile::ColStringLengthFromFits(int nocol) const 513 514 \return number of characters of each data for the column number 'nocol' (if char* typed) of the current BINTABLE extension 515 */ 352 516 353 517 int FitsInFile::ColStringLengthFromFits(int nocol) const … … 365 529 return taille_des_chaines_[index]; 366 530 } 531 532 533 534 /*! \fn void SOPHYA::FitsInFile::GetBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char ** cdata) 535 536 Get the NoLine-th 'line' from the current BINTABLE extension on FITS file, 537 */ 538 367 539 void FitsInFile::GetBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char ** cdata) 368 540 { … … 405 577 } 406 578 579 /*! \fn void SOPHYA::FitsInFile::GetBinTabLine(long NoLine, BnTblLine& ligne) 580 Get the NoLine-th 'line' from the current BINTABLE extension on FITS file, 581 */ 407 582 void FitsInFile::GetBinTabLine(long NoLine, BnTblLine& ligne) 408 583 { … … 447 622 } 448 623 624 /*! \fn void SOPHYA::FitsInFile::GetBinTabLine(int NoLine, float* fdata) 625 626 Get the NoLine-th float 'line' from the current BINTABLE extension on FITS file, 627 */ 449 628 void FitsInFile::GetBinTabLine(int NoLine, float* fdata) 450 629 { … … 466 645 467 646 647 /*! \fn void SPOPHYA::FitsInFile::GetBinTabFCol(double* valeurs,int nentries, int NoCol) const 648 649 fill the array 'valeurs' with double data from the current BINTABLE extension on FITS file, from column number 'NoCol' 650 651 \param <nentries> number of data to be read 652 */ 468 653 void FitsInFile::GetBinTabFCol(double* valeurs,int nentries, int NoCol) const 469 654 { … … 499 684 } 500 685 686 /*! \fn void SOPHYA::FitsInFile::GetBinTabFCol(float* valeurs,int nentries, int NoCol) const 687 688 same as previous method with float data 689 */ 501 690 void FitsInFile::GetBinTabFCol(float* valeurs,int nentries, int NoCol) const 502 691 { … … 531 720 } 532 721 722 /*! \fn void SOPHYA::FitsInFile::GetBinTabFCol(int* valeurs,int nentries, int NoCol) const 723 724 same as previous method with int data 725 */ 726 533 727 void FitsInFile::GetBinTabFCol(int* valeurs,int nentries, int NoCol) const 534 728 { … … 561 755 } 562 756 757 /*! \fn void SOPHYA::FitsInFile::GetBinTabFCol(char** valeurs, int nentries, int NoCol) const 758 759 same as previous method with char* data 760 */ 761 563 762 void FitsInFile::GetBinTabFCol(char** valeurs, int nentries, int NoCol) const 564 763 { … … 587 786 } 588 787 788 /*! \fn void SOPHYA::FitsInFile::GetSingleColumn(double* map, int nentries) const 789 fill the array 'map' with double data from the current extension on FITS file. 790 If the extension is BINTABLE, the first column is provided. 791 792 \param <nentries> number of data to be read 793 */ 589 794 void FitsInFile::GetSingleColumn(double* map, int nentries) const 590 795 { … … 619 824 } 620 825 826 /*! \fn void SOPHYA::FitsInFile::GetSingleColumn(float* map, int nentries) const 827 same as above with float data 828 */ 621 829 void FitsInFile::GetSingleColumn(float* map, int nentries) const 622 830 { … … 649 857 } 650 858 859 /*! \fn void SOPHYA::FitsInFile::GetSingleColumn( int* map, int nentries) const 860 same as above with int data 861 */ 651 862 void FitsInFile::GetSingleColumn( int* map, int nentries) const 652 863 { … … 866 1077 } 867 1078 1079 1080 /*! 1081 \class SOPHYA::FitsOutFile 1082 Class for loading SOPHYA objects from FITS Format Files (uses cfitsio lib) 1083 */ 1084 868 1085 FitsOutFile::FitsOutFile() 869 1086 { … … 871 1088 } 872 1089 1090 /*! \fn SOPHYA::FitsOutFile::FitsOutFile(char flnm[], WriteMode wrm) 1091 1092 \param <WriteMode> enum , WriteMode = clear -> if alreadyy exists, the file will be overwritten (else created) ; WriteMode = append -> further objects will be appended to the file if it exists (else : file created). WriteMode = unknown -> file created if does not exist, else : exception. (the last situation is the default) 1093 1094 */ 873 1095 FitsOutFile::FitsOutFile(char flnm[], WriteMode wrm) 874 1096 { … … 929 1151 930 1152 1153 /*! \fn void SOPHYA::FitsOutFile::makeHeaderImageOnFits(char type, int nbdim, int* naxisn, DVList &dvl) 1154 1155 create an IMAGE header on FITS file. 1156 \param <type> type of data (see method ColTypeFromFits) 1157 \param <nbdim> number of dimensions : 1D, 2D, 3D etc. = NAXIS 1158 \param <naxisn> array containind sizes of the different dimensions 1159 */ 931 1160 void FitsOutFile::makeHeaderImageOnFits(char type, int nbdim, int* naxisn, DVList &dvl) 932 1161 { … … 968 1197 969 1198 } 1199 1200 1201 /*! \fn void SOPHYA::FitsOutFile::PutImageToFits(int nbData, double* map) const 1202 1203 write double data from array 'map'on an IMAGE extension 1204 \param <nbData> number of data to be written 1205 */ 970 1206 void FitsOutFile::PutImageToFits(int nbData, double* map) const 971 1207 { … … 977 1213 } 978 1214 1215 /*! \fn void SOPHYA::FitsOutFile::PutImageToFits(int nbData, float* map) const 1216 1217 same as previous method with float data 1218 */ 979 1219 void FitsOutFile::PutImageToFits(int nbData, float* map) const 980 1220 { … … 986 1226 987 1227 } 1228 1229 /*! \fn void SOPHYA::FitsOutFile::PutImageToFits( int nbData, int* map) const 1230 1231 same as previous method with int data */ 988 1232 void FitsOutFile::PutImageToFits( int nbData, int* map) const 989 1233 { … … 998 1242 999 1243 1244 /*! \fn void SOPHYA::FitsOutFile::makeHeaderBntblOnFits( string fieldType, vector<string> Noms, int nentries, int tfields, DVList &dvl, string extname, vector<int> taille_des_chaines) 1245 1246 create an BINTABLE header on FITS file. 1247 \param <fieldType> array conta 1248 ining characters denoting types of the different column (see method ColTypeFromFits) 1249 \param <Noms> array of the names of columns 1250 \param <nentries> number of data of each column 1251 \param <tfields> number of columns 1252 \param <dvl> a SOPHYA DVList containing keywords to be appended 1253 \param <extname> keyword EXTNAME for FITS file 1254 \param <taille_des_chaines> vector containing the number of characters of data for each char* typed column, with order of appearance in 'fieldType' 1255 */ 1000 1256 void FitsOutFile::makeHeaderBntblOnFits( string fieldType, vector<string> Noms, int nentries, int tfields, DVList &dvl, string extname, vector<int> taille_des_chaines) 1001 1257 { … … 1092 1348 } 1093 1349 1350 /*! \fn void SOPHYA::FitsOutFile::PutColToFits(int nocol, int nentries, double* donnees) const 1351 1352 write double data from array 'donnees ' on column number 'nocol' of a BINTABLE extension. 1353 \param <nentries> number of data to be written 1354 */ 1094 1355 void FitsOutFile::PutColToFits(int nocol, int nentries, double* donnees) const 1095 1356 { … … 1114 1375 if( status ) printerror( status,"erreur ecriture du fichier fits" ); 1115 1376 } 1377 1378 1379 1380 /*! \fn void SOPHYA::FitsOutFile::PutColToFits(int nocol, int nentries, float* donnees) const 1381 1382 same as previous method with float data 1383 */ 1116 1384 void FitsOutFile::PutColToFits(int nocol, int nentries, float* donnees) const 1117 1385 { … … 1140 1408 if( status ) printerror( status,"erreur ecriture du fichier fits" ); 1141 1409 } 1410 1411 1412 /*! \fn void FitsOutFile::PutColToFits(int nocol, int nentries, int* donnees) const 1413 1414 same as previous method with int data 1415 */ 1142 1416 void FitsOutFile::PutColToFits(int nocol, int nentries, int* donnees) const 1143 1417 { … … 1166 1440 if( status ) printerror( status," ecriture du fichier fits" ); 1167 1441 } 1442 1443 1444 /*! \fn void SOPHYA::FitsOutFile::PutColToFits(int nocol, int nentries, char** donnees) const 1445 same as previous method with char* data 1446 */ 1168 1447 void FitsOutFile::PutColToFits(int nocol, int nentries, char** donnees) const 1169 1448 { … … 1228 1507 1229 1508 1509 /* \fn void SOPHYA::FitsOutFile::DVListIntoPrimaryHeader(DVList& dvl) const 1510 1511 Put keywords from a DVList into the primary header of the fits-file 1512 */ 1230 1513 void FitsOutFile::DVListIntoPrimaryHeader(DVList& dvl) const 1231 1514 { -
trunk/SophyaExt/FitsIOServer/fitsfile.h
r1209 r1218 10 10 #define LEN_KEYWORD 9 11 11 12 // classes for saving/loading SOPHYA objects to/from FITS files... 13 // Guy le Meur (september 2000) 14 15 12 16 namespace SOPHYA { 13 17 14 struct BnTblLine 15 { 16 BnTblLine() {} 17 void setFormat(int dc, int fc, int ic, int cc, vector<string> names) 18 { 19 int nbcols = dc + fc + ic + cc; 20 int maxName = names.size(); 21 if (nbcols != maxName) 22 { 23 cout << " WARNING: BnTblLine:: length of vector of column names not equal to total number of columns" << endl; 24 maxName = nbcols < maxName ? nbcols : maxName; 25 } 26 ColName_ = vector<string>(nbcols); 27 for (int k=0; k < maxName; k++) ColName_[k] = names[k]; 28 if (dc >0) ddata_ = vector<double>(dc); 29 if (fc >0) fdata_ = vector<float>(fc); 30 if (ic >0) idata_ = vector<int>(fc); 31 if (cc >0) cdata_ = vector<string>(fc); 32 } 33 34 bool sameFormat(const BnTblLine& btl) const 35 { 36 if (btl.ddata_.size() == ddata_.size() && btl.fdata_.size() == fdata_.size() && btl.idata_.size() == idata_.size() && btl.cdata_.size() == cdata_.size()) return true; 37 else return false; 38 } 39 40 void Print() 41 { 42 int k; 43 cout << " ********* ligne ************* " << endl; 44 cout << " *** noms de variables " << endl; 45 for (k=0; k < ColName_.size(); k++) cout << ColName_[k] << " "; 46 cout << endl; 47 cout << " *** variables doubles " << endl; 48 for (k=0; k < ddata_.size(); k++) cout << ddata_[k] << " "; 49 cout << endl; 50 cout << " *** variables float " << endl; 51 for (k=0; k < fdata_.size(); k++) cout << fdata_[k] << " "; 52 cout << endl; 53 cout << " *** variables int " << endl; 54 for (k=0; k < idata_.size(); k++) cout << idata_[k] << " "; 55 cout << endl; 56 cout << " *** variables string " << endl; 57 for (k=0; k < cdata_.size(); k++) cout << cdata_[k] << " "; 58 cout << endl; 59 cout << " ***************************** " << endl; 60 } 61 vector<double> ddata_; 62 vector<float> fdata_; 63 vector<int> idata_; 64 vector<string> cdata_; 65 vector<string> ColName_; 66 67 }; 68 69 class FitsFile; 70 class FitsInFile; 71 class FitsOutFile; 72 enum WriteMode {append, clear, unknown}; 73 18 struct BnTblLine; 19 class FitsFile; 20 class FitsInFile; 21 class FitsOutFile; 22 enum WriteMode {append, clear, unknown}; 23 74 24 75 25 // 76 26 //! Class for managing Interface for SOPHYA objects to FITS Format Files (uses cfitsio lib) 77 27 78 /*! 79 The class structure is analogous to Sophya-PPersist system : 80 Each SOPHYA object XXX is associated with a object of class FITS_XXX 81 (inheriting from FitsFileHandler), to which input/output operations with FITS 82 files are delegated (through a class Hierarchy : FitsFile (virtual), 83 FitsInFile, FitsOutFile) . A typical example of use is the following : 84 85 \verbatim 86 int m=... ; 87 SphereHEALPix<r_8> sphere1(m); // definition of the SOPHYA object 88 .... fill the sphere .... 89 90 FITS_SphereHEALPix<r_8> fits_sph1(sphere1); 91 // delegated object 92 fits_sph.Write("myfile.fits"); // writing on FITS file 93 94 FITS_SphereHEALPix<r_8> fits_sph2("myfile.fits"); 95 // load a delegated object 96 // from FITS file 97 SphereHEALPix<r_8> sphere2=(SphereHEALPix<r_8>)fits_sph2; 98 // casting the delegated object 99 // into a SOPHYA object 100 \endverbatim 101 102 */ 103 104 class FitsIOHandler { 105 106 107 public: 108 109 virtual ~FitsIOHandler() {} 110 /*! 111 this method is called from inherited objects : 112 113 opens a file 'flnm' 114 115 gets parameters in extension-header (hdunum) 116 117 calls the method 'ReadFromFits' from the inherited object 118 119 120 */ 121 void Read(char flnm[],int hdunum= 0); 122 /*! 123 this method is called from inherited objects : 124 125 for writing a new object in a new fits-extension : 126 127 ??? 128 129 at the end of 130 131 the existing file (flnm), if OldFile=true. 132 133 If OldFile=false, an exception occurs 28 class FitsIOHandler { 29 30 31 public: 32 33 virtual ~FitsIOHandler() {} 34 void Read(char flnm[],int hdunum= 0); 35 void Write(char flnm[]) ; 36 void Read(FitsInFile& ifts, int hdunum=0); 37 void Write(FitsOutFile& ofts) ; 38 39 40 protected: 41 42 virtual void ReadFromFits(FitsInFile& is)=0; 43 virtual void WriteToFits(FitsOutFile& os) =0; 44 45 friend class FitsInFile; 46 friend class FitsOutFile; 47 }; 48 49 50 //! Class (virtual) for managing FITS format files 51 class FitsFile { 52 53 public: 54 55 FitsFile() { InitNull(); }; 56 virtual ~FitsFile(); 57 static string GetErrStatus(int status); 58 inline int statusF() const { return fits_status_;} 134 59 135 By convention, primary header does not contain fits-image data : i.e. 136 all data are fits-extensions. The first relevant header will have hdunum=2. 137 For switching off this convention use the method : 138 139 firstImageOnPrimaryHeader() (see below) 140 141 In that case do not forget to precise hdunum=1 when reading data on primary header. 142 143 calls the method 'WriteToFits' from the inherited object 144 145 */ 146 void Write(char flnm[]) ; 147 148 /*! 149 Read the data on extension hdunum (or primary header, if hdunum=1) from FitsInFIle. With default value for hdunum, one reads the next extension, with respect to the current position. 150 */ 151 void Read(FitsInFile& ifts, int hdunum=0); 152 void Write(FitsOutFile& ofts) ; 153 154 155 protected: 156 virtual void ReadFromFits(FitsInFile& is)=0; 157 virtual void WriteToFits(FitsOutFile& os) =0; 158 friend class FitsInFile; 159 friend class FitsOutFile; 160 private : 161 }; 162 163 164 165 class FitsFile 166 { 167 168 public: 169 170 FitsFile() 171 { 172 InitNull(); 173 }; 174 virtual ~FitsFile(); 175 176 static string GetErrStatus(int status); 177 178 179 180 181 inline int statusF() const { return fits_status_;} 182 183 184 protected: 185 186 void ResetStatus(int& status) ; 60 61 protected: 62 63 void ResetStatus(int& status) ; 187 64 static void printerror(int&) ; 188 65 static void printerror(int&,char* texte) ; 189 inline void InitNull() {fptr_ = NULL; hdutype_= 0; hdunum_ = 1;66 inline void InitNull() {fptr_ = NULL; hdutype_= 0; hdunum_ = 1; 190 67 fits_status_ = 0;} 191 68 192 //! pointer to the FITS file, defined in fitsio.h 193 fitsfile *fptr_; 194 195 //! image or bintable ? 196 int hdutype_; 197 198 //! index of header to be read/written 199 int hdunum_; 200 201 //! last status returned by fitsio library. updated only by several methods 202 int fits_status_; 203 204 205 }; 206 69 fitsfile *fptr_; /**< pointer to the FITS file, defined in fitsio.h */ 70 int hdutype_; /**< image or bintable ? */ 71 int hdunum_; /**< index of header to be read/written */ 72 int fits_status_; /**< last status returned by fitsio library. updated only by several methods */ 73 74 }; 75 76 //! Class for saving SOPHYA objects on FITS Format Files (uses cfitsio lib) 207 77 208 78 class FitsInFile : public FitsFile { … … 210 80 public: 211 81 FitsInFile(); 212 // FitsInFile(char flnm[], int hdunum=0);213 82 FitsInFile(char flnm[]); 214 83 ~FitsInFile() { ; }; 215 84 216 217 ////////////////////////////////////////////////////////// 218 // methods with general purpose 219 /////////////////////////////////////// 220 221 222 223 static int NbBlocks(char flnm[]); 224 static void GetBlockType(char flnm[], int hdunum, string& typeOfExtension, int& naxis, vector<int>& naxisn, string& dataType, DVList& dvl ); 225 226 227 228 // void ReadFInit(char flnm[],int hdunum=0); 229 void ReadFInit(int hdunum); 85 static int NbBlocks(char flnm[]); 86 static void GetBlockType(char flnm[], int hdunum, string& typeOfExtension, int& naxis, vector<int>& naxisn, string& dataType, DVList& dvl ); 87 void ReadFInit(int hdunum); 230 88 231 /*! return a reference on a DVList containing the keywords from FITS file 232 */ 233 inline const DVList& DVListFromFits() const { return dvl_;} 234 235 /* get the keywords of primary header in a DVList */ 236 DVList DVListFromPrimaryHeader() const; 237 238 void moveToFollowingHeader(); 89 /*! \return a reference on a DVList containing the keywords from FITS file */ 90 inline const DVList& DVListFromFits() const { return dvl_;} 91 92 DVList DVListFromPrimaryHeader() const; 93 void moveToFollowingHeader(); 94 95 239 96 240 97 … … 242 99 /////// methods for managing extensions //////////////// 243 100 ////////////////////////////////////////////////////////// 244 245 101 246 102 … … 250 106 251 107 252 /*! return true if the current header corresponds to a FITS image extension */ 108 109 /*! \return true if the current header corresponds to a FITS image extension */ 253 110 inline bool IsFitsImage() const { return (hdutype_ == IMAGE_HDU);} 254 111 255 112 256 113 257 /*! number of dimensions of an image extension : NAXIS parameter (in FITS notations) 258 */ 114 /*! \return number of dimensions of an image extension : NAXIS parameter (in FITS notations) */ 259 115 inline int nbDimOfImage() const {return naxis_;} 260 /*! a reference on a vector containing sizes of the NAXIS dimensions : NAXIS1, NAXIS2, NAXIS3 wtc. 261 */116 117 /*! \return a reference on a vector containing sizes of the NAXIS dimensions : NAXIS1, NAXIS2, NAXIS3 etc. */ 262 118 inline const vector<int>& dimOfImageAxes() const { return naxisn_;} 263 /*! 264 total number of data in the current IMAGE extension 265 */119 120 121 /*! \return total number of data in the current IMAGE extension */ 266 122 inline int nbOfImageData() const { return nbData_; } 267 123 … … 275 131 276 132 277 /*! return true if the current header corresponds to a FITS ASCII or BINTABLE extension */133 /*! \return true if the current header corresponds to a FITS ASCII or BINTABLE extension */ 278 134 inline bool IsFitsTable() const {return (hdutype_ == ASCII_TBL || hdutype_ == BINARY_TBL);} 279 135 280 136 281 static void GetBinTabParameters(fitsfile* fileptr, int& nbcols, int& nrows,137 static void GetBinTabParameters(fitsfile* fileptr, int& nbcols, int& nrows, 282 138 vector<int>& repeat, 283 139 vector<string>& noms, 284 140 vector<char>& types, 285 141 vector<int>& taille_des_chaines); 286 287 288 289 /*! return a character denoting data type of column number 'nocol' in a BINTABLE : 290 291 D : double 292 293 E : float 294 295 I : integer 296 297 S : character string 298 299 */ 300 char ColTypeFromFits(int nocol) const; 301 /*! name of the column number 'nocol' of the current BINTABLE extension 302 */ 303 string ColNameFromFits(int nocol) const; 304 305 /*! number of characters of each data for the column number 'nocol' (if char* typed) of the current BINTABLE extension 306 */ 307 int ColStringLengthFromFits(int nocol) const; 308 309 310 311 /*! 312 Get the NoLine-th 'line' from the current BINTABLE extension on FITS file, 313 */ 314 void GetBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char 142 char ColTypeFromFits(int nocol) const; 143 string ColNameFromFits(int nocol) const; 144 int ColStringLengthFromFits(int nocol) const; 145 void GetBinTabLine(int NoLine, double* ddata, float* fdata, int* idata, char 315 146 ** cdata) ; 316 /*! 317 Get the NoLine-th 'line' from the current BINTABLE extension on FITS file, 318 */ 319 void GetBinTabLine(long NoLine, BnTblLine& ligne) ; 320 321 /*! 322 Get the NoLine-th 'line' from the current BINTABLE extension on FITS file, 323 */ 324 void GetBinTabLine(int NoLine, float* fdata) ; 325 326 /*! 327 fill the array 'valeurs' with double data from the current BINTABLE extension on FITS file, from column number 'NoCol' 328 329 \param <nentries> number of data to be read 330 */ 331 void GetBinTabFCol(double* valeurs, int nentries, int NoCol) const; 332 333 /*! same as previous method with float data */ 334 void GetBinTabFCol(float* valeurs, int nentries, int NoCol) const; 335 /*! same as previous method with int data */ 336 void GetBinTabFCol(int* valeurs, int nentries, int NoCol) const; 337 /*! same as previous method with char* data */ 338 void GetBinTabFCol(char** valeurs,int nentries, int NoCol) const; 339 // Write elements into the FITS data array 147 void GetBinTabLine(long NoLine, BnTblLine& ligne) ; 148 void GetBinTabLine(int NoLine, float* fdata) ; 149 void GetBinTabFCol(double* valeurs, int nentries, int NoCol) const; 150 void GetBinTabFCol(float* valeurs, int nentries, int NoCol) const; 151 void GetBinTabFCol(int* valeurs, int nentries, int NoCol) const; 152 void GetBinTabFCol(char** valeurs,int nentries, int NoCol) const; 340 153 341 154 ///////////////////////////////////////////////////////////// … … 343 156 //////////////////////////////////////////////////////// 344 157 345 /*! return number of columns (return 1 if IMAGE) */346 158 int NbColsFromFits() const; 347 /*! number of data in the current IMAGE extension on FITS file, or number348 of data of column number 'nocol' of the current BINTABLE extension349 */350 159 int NentriesFromFits(int nocol) const; 351 160 352 161 353 /*!354 fill the array 'map' with double data from the current extension on FITS file.355 If the extension is BINTABLE, the first column is provided.356 357 \param <nentries> number of data to be read358 */359 162 void GetSingleColumn(double* map, int nentries) const; 360 163 361 /*! same as above with float data */362 164 void GetSingleColumn(float* map, int nentries) const; 363 165 364 /*! same as above with int data */365 166 void GetSingleColumn(int* map, int nentries) const; 366 167 … … 371 172 static void GetImageParameters (fitsfile* fileptr,int& bitpix,int& naxis,vector<int>& naxisn); 372 173 373 374 375 376 //! fits-Image parameter 377 int bitpix_; 378 379 //! fits-Image parameter 380 int naxis_; 381 382 //! fits-Image parameters : sizes of dimensions 383 vector<int> naxisn_; 384 385 //! fits-Image parameter: number of data 386 int nbData_; 387 388 //! Bintable parameter 389 int nrows_; 390 391 //! Bintable parameter 392 vector<int> repeat_; 393 394 //! Bintable parameter 395 int nbcols_; 396 397 //! Bintable parameter: column names 398 vector<string> noms_; 399 400 //! Bintable parameters: types of columns (D: double, E: float, I: integers, A: char*) 401 vector<char> types_; 402 403 //! Bintable parameters: length of the char* variables 404 vector<int> taille_des_chaines_; 405 406 //! DVList for transferring keywords 407 DVList dvl_; 408 409 410 }; 411 174 int bitpix_; /**< fits-Image parameter */ 175 int naxis_; /**< fits-Image parameter */ 176 vector<int> naxisn_; /**< fits-Image parameters : sizes of dimensions */ 177 int nbData_; /*< fits-Image parameter: number of data */ 178 int nrows_; /**< Bintable parameter */ 179 vector<int> repeat_; /**< Bintable parameter */ 180 int nbcols_; /**< Bintable parameter */ 181 vector<string> noms_; /**< Bintable parameter: column names */ 182 vector<char> types_; /**< Bintable parameters: types of columns (D: double, E: float, I: integers, A: char*) */ 183 DVList dvl_; /**< DVList for transferring keywords */ 184 vector<int> taille_des_chaines_; /**< Bintable parameters: length of the char* variables */ 185 186 }; 187 188 //! Class for loading SOPHYA objects from FITS Format Files (uses cfitsio lib) 412 189 413 190 class FitsOutFile : public FitsFile { … … 417 194 418 195 FitsOutFile(); 419 /*!420 \param <WriteMode> enum , WriteMode = clear -> if alreadyy exists, the file will be overwritten (else created) ; WriteMode = append -> further objects will be appended to the file if it exists (else : file created). WriteMode = unknown -> file created if does not exist, else : exception. (the last situation is the default)421 422 */423 196 FitsOutFile(char flnm[], WriteMode wrm = unknown ); 424 197 ~FitsOutFile() { ;}; … … 437 210 438 211 inline void firstImageOnPrimaryHeader() {imageOnPrimary_=true;} 439 440 /*! create an IMAGE header on FITS file. 441 \param <type> type of data (see method ColTypeFromFits) 442 \param <nbdim> number of dimensions : 1D, 2D, 3D etc. = NAXIS 443 \param <naxisn> array containind sizes of the different dimensions 444 */ 445 void makeHeaderImageOnFits(char type, int nbdim, int* naxisn, DVList &dvl) ; 446 447 448 /*! write double data from array 'map'on an IMAGE extension 449 \param <nbData> number of data to be written 450 451 */ 452 void PutImageToFits( int nbData, double* map) const; 453 454 /*! same as previous method with float data */ 455 void PutImageToFits(int nbData, float* map ) const; 456 457 /*! same as previous method with int data */ 458 void PutImageToFits(int nbData, int* map) const; 212 void makeHeaderImageOnFits(char type, int nbdim, int* naxisn, DVList &dvl) ; 213 void PutImageToFits( int nbData, double* map) const; 214 void PutImageToFits(int nbData, float* map ) const; 215 void PutImageToFits(int nbData, int* map) const; 459 216 460 217 … … 466 223 467 224 468 /*! create an BINTABLE header on FITS file.469 \param <fieldType> array conta470 ining characters denoting types of the different column (see method ColTypeFromFits)471 \param <Noms> array of the names of columns472 \param <nentries> number of data of each column473 \param <tfields> number of columns474 \param <dvl> a SOPHYA DVList containing keywords to be appended475 \param <extname> keyword EXTNAME for FITS file476 \param <taille_des_chaines> vector containing the number of characters of data for each char* typed column, with order of appearance in 'fieldType'477 */478 225 void makeHeaderBntblOnFits ( string fieldType, vector<string> Noms, int nentries, int tfields, DVList &dvl, string extname, vector<int> taille_des_chaines) ; 479 480 /*! write double data from array 'donnees ' on column number 'nocol' of a BINTABLE extension. 481 \param <nentries> number of data to be written 482 483 */ 484 void PutColToFits(int nocol, int nentries, double* donnees) const; 485 486 /*! same as previous method with float data */ 487 void PutColToFits(int nocol, int nentries, float* donnees) const; 488 489 /*! same as previous method with int data */ 490 void PutColToFits(int nocol, int nentries, int* donnees) const; 491 492 /*! same as previous method with char* data */ 493 void PutColToFits(int nocol, int nentries, char** donnees) const; 494 495 void PutBinTabLine(long NoLine, BnTblLine& ligne) const; 226 void PutColToFits(int nocol, int nentries, double* donnees) const; 227 void PutColToFits(int nocol, int nentries, float* donnees) const; 228 void PutColToFits(int nocol, int nentries, int* donnees) const; 229 void PutColToFits(int nocol, int nentries, char** donnees) const; 230 void PutBinTabLine(long NoLine, BnTblLine& ligne) const; 496 231 497 232 … … 501 236 502 237 503 /* Put keywords from a DVList into the primary header of the fits-file */504 238 void DVListIntoPrimaryHeader(DVList& dvl) const; 505 239 … … 515 249 }; 516 250 251 struct BnTblLine 252 { 253 BnTblLine() {} 254 void setFormat(int dc, int fc, int ic, int cc, vector<string> names); 255 bool sameFormat(const BnTblLine& btl) const; 256 257 void Print(); 258 259 vector<double> ddata_; 260 vector<float> fdata_; 261 vector<int> idata_; 262 vector<string> cdata_; 263 vector<string> ColName_; 264 265 }; 517 266 518 267 -
trunk/SophyaLib/Samba/lambdaBuilder.cc
r729 r1218 1 1 #include "lambdaBuilder.h" 2 2 #include "nbconst.h" 3 4 5 /*! 6 \class SOPHYA::Legendre 7 generate Legendre polynomials : use in two steps : 8 9 a) instanciate Legendre(\f$x\f$, \f$lmax\f$) ; \f$x\f$ is the value for wich Legendre polynomials will be required (usually equal to \f$\cos \theta\f$) and \f$lmax\f$ is the MAXIMUM value of the order of polynomials wich will be required in the following code (all polynomials, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula). 10 11 b) get the value of Legendre polynomial for a particular value of \f$l\f$ by calling the method getPl. 12 13 */ 3 14 4 15 … … 12 23 array_init(lmax); 13 24 } 25 26 /*! \fn void SOPHYA::Legendre::array_init(int_4 lmax) 27 28 compute all \f$P_l(x,l_{max})\f$ for \f$l=1,l_{max}\f$ 29 */ 14 30 void Legendre::array_init(int_4 lmax) 15 31 { … … 28 44 TVector<r_8>* LambdaLMBuilder::normal_l_ = NULL; 29 45 46 47 48 /*! \class SOPHYA::LambdaLMBuilder 49 50 51 This class generate the coefficients : 52 \f[ 53 \lambda_l^m=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}} 54 P_l^m(\cos{\theta}) 55 \f] 56 where \f$P_l^m\f$ are the associated Legendre polynomials. The above coefficients contain the theta-dependance of spheric harmonics : 57 \f[ 58 Y_{lm}(\cos{\theta})=\lambda_l^m(\cos{\theta}) e^{im\phi}. 59 \f] 60 61 Each object has a fixed theta (radians), and maximum l and m to be calculated 62 (lmax and mmax). 63 use the class in two steps : 64 a) instanciate LambdaLMBuilder(\f$\theta\f$, \f$lmax\f$, \f$mmax\f$) ; \f$lmax\f$ and \f$mmax\f$ are MAXIMUM values for which \f$\lambda_l^m\f$ will be required in the following code (all coefficients, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula). 65 b) get the values of coefficients for particular values of \f$l\f$ and \f$m\f$ by calling the method lamlm. 66 */ 67 68 30 69 LambdaLMBuilder::LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax) 31 70 { … … 84 123 85 124 125 /*! \fn void SOPHYA::LambdaLMBuilder::updateArrayRecurrence(int_4 lmax) 126 127 compute a static array of coefficients independant from theta (common to all instances of the LambdaBuilder Class 128 */ 86 129 void LambdaLMBuilder::updateArrayRecurrence(int_4 lmax) 87 130 { … … 98 141 } 99 142 100 143 /*! \fn void SOPHYA::LambdaLMBuilder::updateArrayLamNorm() 144 145 compute static arrays of coefficients independant from theta (common to all instances of the derived classes 146 */ 101 147 void LambdaLMBuilder::updateArrayLamNorm() 102 148 { … … 119 165 120 166 167 /*! \class SOPHYA::LambdaWXBuilder 168 169 This class generates the coefficients : 170 \f[ 171 _{w}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}} G^+_{lm} 172 \f] 173 \f[ 174 _{x}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}G^-_{lm} 175 \f] 176 where 177 \f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta}) 178 \f] 179 and 180 \f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right) 181 \f] 182 \f$P_l^m\f$ are the associated Legendre polynomials. 183 184 The coefficients express the theta-dependance of the \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ functions : 185 \f[W_{lm}(\cos{\theta}) = \sqrt{\frac{(l+2)!}{2(l-2)!}}_w\lambda_l^m(\cos{\theta})e^{im\phi} 186 \f] 187 \f[X_{lm}(\cos{\theta}) = -i\sqrt{\frac{(l+2)!}{2(l-2)!}}_x\lambda_l^m(\cos{\theta})e^{im\phi} 188 \f] 189 where \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ are defined as : 190 191 \f[ 192 W_{lm}(\cos{\theta})=-\frac{1}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left( 193 _{+2}Y_l^m(\cos{\theta})+_{-2}Y_l^m(\cos{\theta})\right) 194 \f] 195 \f[X_{lm}(\cos{\theta})=-\frac{i}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left( 196 _{+2}Y_l^m(\cos{\theta})-_{-2}Y_l^m(\cos{\theta})\right) 197 \f] 198 199 */ 121 200 122 201 … … 173 252 } 174 253 254 /*! \class SOPHYA::LambdaPMBuilder 255 256 This class generates the coefficients 257 \f[ 258 _{\pm}\lambda_l^m=2\sqrt{\frac{(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\left( G^+_{lm} \mp G^-_{lm}\right) 259 \f] 260 where 261 \f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta}) 262 \f] 263 and 264 \f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right) 265 \f] 266 and \f$P_l^m\f$ are the associated Legendre polynomials. 267 The coefficients express the theta-dependance of the spin-2 spherical harmonics : 268 \f[_{\pm2}Y_l^m(\cos{\theta})=_\pm\lambda_l^m(\cos{\theta})e^{im\phi} 269 \f] 270 */ 175 271 176 272 LambdaPMBuilder::LambdaPMBuilder(r_8 theta, int_4 lmax, int_4 mmax) : LambdaLMBuilder(theta, lmax, mmax) -
trunk/SophyaLib/Samba/lambdaBuilder.h
r864 r1218 10 10 namespace SOPHYA { 11 11 12 /*! 13 generate Legendre polynomials : use in two steps : 14 15 a) instanciate Legendre(\f$x\f$, \f$lmax\f$) ; \f$x\f$ is the value for wich Legendre polynomials will be required (usually equal to \f$\cos \theta\f$) and \f$lmax\f$ is the MAXIMUM value of the order of polynomials wich will be required in the following code (all polynomials, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula). 16 17 b) get the value of Legendre polynomial for a particular value of \f$l\f$ by calling the method getPl. 18 */ 12 /*! classe pour les polynomes de legendre*/ 19 13 class Legendre { 20 14 … … 32 26 33 27 private : 34 /*! compute all \f$P_l(x,l_{max})\f$ for \f$l=1,l_{max}\f$ */35 28 void array_init(int_4 lmax); 36 29 … … 42 35 43 36 44 /*!45 This class generate the coefficients :46 \f[47 \lambda_l^m=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}48 P_l^m(\cos{\theta})49 \f]50 where \f$P_l^m\f$ are the associated Legendre polynomials. The above coefficients contain the theta-dependance of spheric harmonics :51 \f[52 Y_{lm}(\cos{\theta})=\lambda_l^m(\cos{\theta}) e^{im\phi}.53 \f]54 55 Each object has a fixed theta (radians), and maximum l and m to be calculated56 (lmax and mmax).57 use the class in two steps :58 a) instanciate LambdaLMBuilder(\f$\theta\f$, \f$lmax\f$, \f$mmax\f$) ; \f$lmax\f$ and \f$mmax\f$ are MAXIMUM values for which \f$\lambda_l^m\f$ will be required in the following code (all coefficients, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula).59 b) get the values of coefficients for particular values of \f$l\f$ and \f$m\f$ by calling the method lamlm.60 */61 37 class LambdaLMBuilder { 62 38 … … 71 47 72 48 private: 73 /*! compute a static array of coefficients independant from theta (common to all instances of the LambdaBuilder Class */74 49 void updateArrayRecurrence(int_4 lmax); 75 50 … … 82 57 protected : 83 58 84 /*! compute static arrays of coefficients independant from theta (common to all instances of the derived classes */85 59 void updateArrayLamNorm(); 86 60 … … 96 70 97 71 98 /*!99 100 This class generates the coefficients :101 \f[102 _{w}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}} G^+_{lm}103 \f]104 \f[105 _{x}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}G^-_{lm}106 \f]107 where108 \f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})109 \f]110 and111 \f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)112 \f]113 \f$P_l^m\f$ are the associated Legendre polynomials.114 115 The coefficients express the theta-dependance of the \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ functions :116 \f[W_{lm}(\cos{\theta}) = \sqrt{\frac{(l+2)!}{2(l-2)!}}_w\lambda_l^m(\cos{\theta})e^{im\phi}117 \f]118 \f[X_{lm}(\cos{\theta}) = -i\sqrt{\frac{(l+2)!}{2(l-2)!}}_x\lambda_l^m(\cos{\theta})e^{im\phi}119 \f]120 where \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ are defined as :121 122 \f[123 W_{lm}(\cos{\theta})=-\frac{1}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(124 _{+2}Y_l^m(\cos{\theta})+_{-2}Y_l^m(\cos{\theta})\right)125 \f]126 \f[X_{lm}(\cos{\theta})=-\frac{i}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(127 _{+2}Y_l^m(\cos{\theta})-_{-2}Y_l^m(\cos{\theta})\right)128 \f]129 130 */131 72 class LambdaWXBuilder : public LambdaLMBuilder 132 73 { … … 156 97 }; 157 98 158 /*!159 160 This class generates the coefficients161 \f[162 _{\pm}\lambda_l^m=2\sqrt{\frac{(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\left( G^+_{lm} \mp G^-_{lm}\right)163 \f]164 where165 \f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})166 \f]167 and168 \f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)169 \f]170 and \f$P_l^m\f$ are the associated Legendre polynomials.171 The coefficients express the theta-dependance of the spin-2 spherical harmonics :172 \f[_{\pm2}Y_l^m(\cos{\theta})=_\pm\lambda_l^m(\cos{\theta})e^{im\phi}173 \f]174 */175 99 class LambdaPMBuilder : public LambdaLMBuilder 176 100 { -
trunk/SophyaLib/Samba/sphericaltransformserver.cc
r833 r1218 8 8 #include "nbmath.h" 9 9 10 10 /*! \class SOPHYA::SphericalTransformServer 11 12 Class for performing analysis and synthesis of sky maps using spin-0 or spin-2 spherical harmonics. 13 14 Maps must be SOPHYA SphericalMaps (SphereGorski or SphereThetaPhi). 15 16 Temperature and polarization (Stokes parameters) can be developped on spherical harmonics : 17 \f[ 18 \frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n}) 19 \f] 20 \f[ 21 Q(\hat{n})=\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EW_{lm}(\hat{n})+a_{lm}^BX_{lm}(\hat{n})\right) 22 \f] 23 \f[ 24 U(\hat{n})=-\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EX_{lm}(\hat{n})-a_{lm}^BW_{lm}(\hat{n})\right) 25 \f] 26 \f[ 27 \left(Q \pm iU\right)(\hat{n})=\sum_{lm}a_{\pm 2lm}\, _{\pm 2}Y_l^m(\hat{n}) 28 \f] 29 30 \f[ 31 Y_l^m(\hat{n})=\lambda_l^m(\theta)e^{im\phi} 32 \f] 33 \f[ 34 _{\pm}Y_l^m(\hat{n})=_{\pm}\lambda_l^m(\theta)e^{im\phi} 35 \f] 36 \f[ 37 W_{lm}(\hat{n})=\frac{1}{N_l}\,_{w}\lambda_l^m(\theta)e^{im\phi} 38 \f] 39 \f[ 40 X_{lm}(\hat{n})=\frac{-i}{N_l}\,_{x}\lambda_l^m(\theta)e^{im\phi} 41 \f] 42 43 (see LambdaLMBuilder, LambdaPMBuilder, LambdaWXBuilder classes) 44 45 power spectra : 46 47 \f[ 48 C_l^T=\frac{1}{2l+1}\sum_{m=0}^{+ \infty }\left|a_{lm}^T\right|^2=\langle\left|a_{lm}^T\right|^2\rangle 49 \f] 50 \f[ 51 C_l^E=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^E\right|^2=\langle\left|a_{lm}^E\right|^2\rangle 52 \f] 53 \f[ 54 C_l^B=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^B\right|^2=\langle\left|a_{lm}^B\right|^2\rangle 55 \f] 56 57 \arg 58 \b Synthesis : Get temperature and polarization maps from \f$a_{lm}\f$ coefficients or from power spectra, (methods GenerateFrom...). 59 60 \b Temperature: 61 \f[ 62 \frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n}) = \sum_{-\infty}^{+\infty}b_m(\theta)e^{im\phi} 63 \f] 64 65 with 66 \f[ 67 b_m(\theta)=\sum_{l=\left|m\right|}^{+\infty}a_{lm}^T\lambda_l^m(\theta) 68 \f] 69 70 \b Polarisation 71 \f[ 72 Q \pm iU = \sum_{-\infty}^{+\infty}b_m^{\pm}(\theta)e^{im\phi} 73 \f] 74 75 where : 76 \f[ 77 b_m^{\pm}(\theta) = \sum_{l=\left|m\right|}^{+\infty}a_{\pm 2lm}\,_{\pm}\lambda_l^m(\theta) 78 \f] 79 80 or : 81 \f[ 82 Q = \sum_{-\infty}^{+\infty}b_m^{Q}(\theta)e^{im\phi} 83 \f] 84 \f[ 85 U = \sum_{-\infty}^{+\infty}b_m^{U}(\theta)e^{im\phi} 86 \f] 87 88 where: 89 \f[ 90 b_m^{Q}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(a_{lm}^E\,_{w}\lambda_l^m(\theta)-ia_{lm}^B\,_{x}\lambda_l^m(\theta)\right) 91 \f] 92 \f[ 93 b_m^{U}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(ia_{lm}^E\,_{x}\lambda_l^m(\theta)+a_{lm}^B\,_{w}\lambda_l^m(\theta)\right) 94 \f] 95 96 Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ \f$\frac{\Delta T}{T}\f$, \f$Q\f$,\f$U\f$ can be computed by FFT. 97 98 99 \arg 100 \b Analysis : Get \f$a_{lm}\f$ coefficients or power spectra from temperature and polarization maps (methods DecomposeTo...). 101 102 \b Temperature: 103 \f[ 104 a_{lm}^T=\int\frac{\Delta T}{T}(\hat{n})Y_l^{m*}(\hat{n})d\hat{n} 105 \f] 106 107 approximated as : 108 \f[ 109 a_{lm}^T=\sum_{\theta_k}\omega_kC_m(\theta_k)\lambda_l^m(\theta_k) 110 \f] 111 where : 112 \f[ 113 C_m (\theta _k)=\sum_{\phi _{k\prime}}\frac{\Delta T}{T}(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}} 114 \f] 115 Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ (\f$\omega_k\f$ is the solid angle of each pixel of the slice \f$\theta_k\f$) \f$C_m\f$ can be computed by FFT. 116 117 \b polarisation: 118 119 \f[ 120 a_{\pm 2lm}=\sum_{\theta_k}\omega_kC_m^{\pm}(\theta_k)\,_{\pm}\lambda_l^m(\theta_k) 121 \f] 122 where : 123 \f[ 124 C_m^{\pm} (\theta _k)=\sum_{\phi _{k\prime}}\left(Q \pm iU\right)(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}} 125 \f] 126 or : 127 128 \f[ 129 a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(C_m^{Q}(\theta_k)\,_{w}\lambda_l^m(\theta_k)-iC_m^{U}(\theta_k)\,_{x}\lambda_l^m(\theta_k)\right) 130 \f] 131 \f[ 132 a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(iC_m^{Q}(\theta_k)\,_{x}\lambda_l^m(\theta_k)+C_m^{U}(\theta_k)\,_{w}\lambda_l^m(\theta_k)\right) 133 \f] 134 135 where : 136 \f[ 137 C_m^{Q} (\theta _k)=\sum_{\phi _{k\prime}}Q(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}} 138 \f] 139 \f[ 140 C_m^{U} (\theta _k)=\sum_{\phi _{k\prime}}U(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}} 141 \f] 142 143 */ 144 145 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const 146 147 synthesis of a temperature map from Alm coefficients 148 */ 11 149 template<class T> 12 150 void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const … … 125 263 126 264 265 /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::fourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const 266 267 \return a vector with nph elements which are sums :\f$\sum_{m=-mmax}^{mmax}b_m(\theta)e^{im\varphi}\f$ for nph values of \f$\varphi\f$ regularly distributed in \f$[0,\pi]\f$ ( calculated by FFT) 268 269 The object b_m (\f$b_m\f$) of the class Bm is a special vector which index goes from -mmax to mmax. 270 */ 127 271 template<class T> 128 272 TVector< complex<T> > SphericalTransformServer<T>::fourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const … … 205 349 206 350 //******************************************** 351 /*! \fn TVector<T> SOPHYA::SphericalTransformServer::RfourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const 352 353 same as fourierSynthesisFromB, but return a real vector, taking into account the fact that b(-m) is conjugate of b(m) */ 207 354 template<class T> 208 355 TVector<T> SphericalTransformServer<T>::RfourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const … … 285 432 //******************************************* 286 433 434 /*! \fn Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const 435 436 \return the Alm coefficients from analysis of a temperature map. 437 438 \param<nlmax> : maximum value of the l index 439 440 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut. 441 */ 287 442 template<class T> 288 443 Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const … … 346 501 return alm; 347 502 } 503 /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const 504 505 \return a vector with mmax elements which are sums : 506 \f$\sum_{k=0}^{nphi}datain(\theta,\varphi_k)e^{im\varphi_k}\f$ for (mmax+1) values of \f$m\f$ from 0 to mmax. 507 */ 348 508 template<class T> 349 509 TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const … … 384 544 385 545 //&&&&&&&&& nouvelle version 546 /* \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const 547 548 same as previous one, but with a "datain" which is real (not complex) */ 386 549 template<class T> 387 550 TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const … … 445 608 } 446 609 610 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm(SphericalMap<T>& mapq, 611 SphericalMap<T>& mapu, 612 int_4 pixelSizeIndex, 613 const Alm<T>& alme, 614 const Alm<T>& almb) const 615 616 synthesis of a polarization map from Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */ 447 617 template<class T> 448 618 void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq, … … 521 691 522 692 693 /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq, 694 const SphericalMap<T>& mapu, 695 Alm<T>& alme, 696 Alm<T>& almb, 697 int_4 nlmax, 698 r_8 cos_theta_cut) const 699 700 analysis of a polarization map into Alm coefficients. 701 702 The spheres \c mapq and \c mapu contain respectively the Stokes parameters. 703 704 \c a2lme and \c a2lmb will receive respectively electric and magnetic Alm's 705 nlmax : maximum value of the l index 706 707 \c cos_theta_cut : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut. 708 */ 523 709 template<class T> 524 710 void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq, … … 578 764 579 765 766 /*! \fn void SOPHYA::SphericalTransformServer::almFromWX(int_4 nlmax, int_4 nmmax, 767 r_8 phi0, r_8 domega, 768 r_8 theta, 769 const TVector<T>& dataq, 770 const TVector<T>& datau, 771 Alm<T>& alme, 772 Alm<T>& almb) const 773 774 Compute polarized Alm's as : 775 \f[ 776 a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(\,_{w}\lambda_l^m\tilde{Q}-i\,_{x}\lambda_l^m\tilde{U}\right)} 777 \f] 778 \f[ 779 a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(i\,_{x}\lambda_l^m\tilde{Q}+\,_{w}\lambda_l^m\tilde{U}\right)} 780 \f] 781 782 where \f$\tilde{Q}\f$ and \f$\tilde{U}\f$ are C-coefficients computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters. 783 784 \f$\omega_{pix}\f$ are solid angle of each pixel. 785 786 dataq, datau : Stokes parameters. 787 788 */ 580 789 template<class T> 581 790 void SphericalTransformServer<T>::almFromWX(int_4 nlmax, int_4 nmmax, … … 628 837 629 838 630 template<class T> 631 void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax, 839 /*! \fn void SOPHYA::SphericalTransformServer::almFromPM(int_4 nph, int_4 nlmax, 840 int_4 nmmax, 841 r_8 phi0, r_8 domega, 842 r_8 theta, 843 const TVector<T>& dataq, 844 const TVector<T>& datau, 845 Alm<T>& alme, 846 Alm<T>& almb) const 847 848 Compute polarized Alm's as : 849 \f[ 850 a_{lm}^E=-\frac{1}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}+\,_{-}\lambda_l^m\tilde{P^-}\right)} 851 \f] 852 \f[ 853 a_{lm}^B=\frac{i}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}-\,_{-}\lambda_l^m\tilde{P^-}\right)} 854 \f] 855 856 where \f$\tilde{P^{\pm}}=\tilde{Q}\pm\tilde{U}\f$ computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters,\f$Q\f$ and \f$U\f$ . 857 858 \f$\omega_{pix}\f$ are solid angle of each pixel. 859 860 dataq, datau : Stokes parameters. 861 862 */ 863 template<class T> 864 void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, 865 int_4 nmmax, 632 866 r_8 phi0, r_8 domega, 633 867 r_8 theta, … … 673 907 674 908 909 /*! \fn void SOPHYA::SphericalTransformServer::mapFromWX(int_4 nlmax, int_4 nmmax, 910 SphericalMap<T>& mapq, 911 SphericalMap<T>& mapu, 912 const Alm<T>& alme, 913 const Alm<T>& almb) const 914 915 synthesis of Stokes parameters following formulae : 916 917 \f[ 918 Q=\sum_{m=-mmax}^{mmax}b_m^qe^{im\varphi} 919 \f] 920 \f[ 921 U=\sum_{m=-mmax}^{mmax}b_m^ue^{im\varphi} 922 \f] 923 924 computed by FFT (method fourierSynthesisFromB called by the present one) 925 926 with : 927 928 \f[ 929 b_m^q=-\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(\,_{w}\lambda_l^ma_{lm}^E-i\,_{x}\lambda_l^ma_{lm}^B\right) } 930 \f] 931 \f[ 932 b_m^u=\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(i\,_{x}\lambda_l^ma_{lm}^E+\,_{w}\lambda_l^ma_{lm}^B\right) } 933 \f] 934 */ 675 935 template<class T> 676 936 void SphericalTransformServer<T>::mapFromWX(int_4 nlmax, int_4 nmmax, … … 744 1004 } 745 1005 } 1006 /*! \fn void SOPHYA::SphericalTransformServer::mapFromPM(int_4 nlmax, int_4 nmmax, 1007 SphericalMap<T>& mapq, 1008 SphericalMap<T>& mapu, 1009 const Alm<T>& alme, 1010 const Alm<T>& almb) const 1011 1012 synthesis of polarizations following formulae : 1013 1014 \f[ 1015 P^+ = \sum_{m=-mmax}^{mmax} {b_m^+e^{im\varphi} } 1016 \f] 1017 \f[ 1018 P^- = \sum_{m=-mmax}^{mmax} {b_m^-e^{im\varphi} } 1019 \f] 1020 1021 computed by FFT (method fourierSynthesisFromB called by the present one) 1022 1023 with : 1024 1025 \f[ 1026 b_m^+=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E+ia_{lm}^B \right) } 1027 \f] 1028 \f[ 1029 b_m^-=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E-ia_{lm}^B \right) } 1030 \f] 1031 */ 746 1032 template<class T> 747 1033 void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax, … … 809 1095 810 1096 1097 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sphq, 1098 SphericalMap<T>& sphu, 1099 int_4 pixelSizeIndex, 1100 const TVector<T>& Cle, 1101 const TVector<T>& Clb, 1102 const r_8 fwhm) const 1103 1104 synthesis of a polarization map from power spectra electric-Cl and magnetic-Cl (Alm's are generated randomly, following a gaussian distribution). 1105 \param fwhm FWHM in arcmin for random generation of Alm's (eg. 5) 1106 */ 811 1107 template<class T> 812 1108 void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq, … … 833 1129 GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb); 834 1130 } 1131 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sph, 1132 int_4 pixelSizeIndex, 1133 const TVector<T>& Cl, 1134 const r_8 fwhm) const 1135 1136 synthesis of a temperature map from power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */ 835 1137 template<class T> 836 1138 void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph, … … 846 1148 847 1149 1150 /*! \fn TVector<T> SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const 1151 1152 \return power spectrum from analysis of a temperature map. 1153 1154 \param<nlmax> : maximum value of the l index 1155 1156 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut. 1157 */ 848 1158 template <class T> 849 1159 TVector<T> SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const -
trunk/SophyaLib/Samba/sphericaltransformserver.h
r866 r1218 11 11 namespace SOPHYA { 12 12 13 //14 /*! Class for performing analysis and synthesis of sky maps using spin-0 or spin-2 spherical harmonics.15 13 16 Maps must be SOPHYA SphericalMaps (SphereGorski or SphereThetaPhi).17 18 Temperature and polarization (Stokes parameters) can be developped on spherical harmonics :19 \f[20 \frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n})21 \f]22 \f[23 Q(\hat{n})=\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EW_{lm}(\hat{n})+a_{lm}^BX_{lm}(\hat{n})\right)24 \f]25 \f[26 U(\hat{n})=-\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EX_{lm}(\hat{n})-a_{lm}^BW_{lm}(\hat{n})\right)27 \f]28 \f[29 \left(Q \pm iU\right)(\hat{n})=\sum_{lm}a_{\pm 2lm}\, _{\pm 2}Y_l^m(\hat{n})30 \f]31 32 \f[33 Y_l^m(\hat{n})=\lambda_l^m(\theta)e^{im\phi}34 \f]35 \f[36 _{\pm}Y_l^m(\hat{n})=_{\pm}\lambda_l^m(\theta)e^{im\phi}37 \f]38 \f[39 W_{lm}(\hat{n})=\frac{1}{N_l}\,_{w}\lambda_l^m(\theta)e^{im\phi}40 \f]41 \f[42 X_{lm}(\hat{n})=\frac{-i}{N_l}\,_{x}\lambda_l^m(\theta)e^{im\phi}43 \f]44 45 (see LambdaLMBuilder, LambdaPMBuilder, LambdaWXBuilder classes)46 47 power spectra :48 49 \f[50 C_l^T=\frac{1}{2l+1}\sum_{m=0}^{+ \infty }\left|a_{lm}^T\right|^2=\langle\left|a_{lm}^T\right|^2\rangle51 \f]52 \f[53 C_l^E=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^E\right|^2=\langle\left|a_{lm}^E\right|^2\rangle54 \f]55 \f[56 C_l^B=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^B\right|^2=\langle\left|a_{lm}^B\right|^2\rangle57 \f]58 59 \arg60 \b Synthesis : Get temperature and polarization maps from \f$a_{lm}\f$ coefficients or from power spectra, (methods GenerateFrom...).61 62 \b Temperature:63 \f[64 \frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n}) = \sum_{-\infty}^{+\infty}b_m(\theta)e^{im\phi}65 \f]66 67 with68 \f[69 b_m(\theta)=\sum_{l=\left|m\right|}^{+\infty}a_{lm}^T\lambda_l^m(\theta)70 \f]71 72 \b Polarisation73 \f[74 Q \pm iU = \sum_{-\infty}^{+\infty}b_m^{\pm}(\theta)e^{im\phi}75 \f]76 77 where :78 \f[79 b_m^{\pm}(\theta) = \sum_{l=\left|m\right|}^{+\infty}a_{\pm 2lm}\,_{\pm}\lambda_l^m(\theta)80 \f]81 82 or :83 \f[84 Q = \sum_{-\infty}^{+\infty}b_m^{Q}(\theta)e^{im\phi}85 \f]86 \f[87 U = \sum_{-\infty}^{+\infty}b_m^{U}(\theta)e^{im\phi}88 \f]89 90 where:91 \f[92 b_m^{Q}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(a_{lm}^E\,_{w}\lambda_l^m(\theta)-ia_{lm}^B\,_{x}\lambda_l^m(\theta)\right)93 \f]94 \f[95 b_m^{U}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(ia_{lm}^E\,_{x}\lambda_l^m(\theta)+a_{lm}^B\,_{w}\lambda_l^m(\theta)\right)96 \f]97 98 Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ \f$\frac{\Delta T}{T}\f$, \f$Q\f$,\f$U\f$ can be computed by FFT.99 100 101 \arg102 \b Analysis : Get \f$a_{lm}\f$ coefficients or power spectra from temperature and polarization maps (methods DecomposeTo...).103 104 \b Temperature:105 \f[106 a_{lm}^T=\int\frac{\Delta T}{T}(\hat{n})Y_l^{m*}(\hat{n})d\hat{n}107 \f]108 109 approximated as :110 \f[111 a_{lm}^T=\sum_{\theta_k}\omega_kC_m(\theta_k)\lambda_l^m(\theta_k)112 \f]113 where :114 \f[115 C_m (\theta _k)=\sum_{\phi _{k\prime}}\frac{\Delta T}{T}(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}116 \f]117 Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ (\f$\omega_k\f$ is the solid angle of each pixel of the slice \f$\theta_k\f$) \f$C_m\f$ can be computed by FFT.118 119 \b polarisation:120 121 \f[122 a_{\pm 2lm}=\sum_{\theta_k}\omega_kC_m^{\pm}(\theta_k)\,_{\pm}\lambda_l^m(\theta_k)123 \f]124 where :125 \f[126 C_m^{\pm} (\theta _k)=\sum_{\phi _{k\prime}}\left(Q \pm iU\right)(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}127 \f]128 or :129 130 \f[131 a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(C_m^{Q}(\theta_k)\,_{w}\lambda_l^m(\theta_k)-iC_m^{U}(\theta_k)\,_{x}\lambda_l^m(\theta_k)\right)132 \f]133 \f[134 a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(iC_m^{Q}(\theta_k)\,_{x}\lambda_l^m(\theta_k)+C_m^{U}(\theta_k)\,_{w}\lambda_l^m(\theta_k)\right)135 \f]136 137 where :138 \f[139 C_m^{Q} (\theta _k)=\sum_{\phi _{k\prime}}Q(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}140 \f]141 \f[142 C_m^{U} (\theta _k)=\sum_{\phi _{k\prime}}U(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}143 \f]144 145 */146 14 template <class T> 147 15 class SphericalTransformServer … … 150 18 public: 151 19 152 SphericalTransformServer() 153 { 154 fftIntfPtr_=new FFTPackServer; 155 fftIntfPtr_->setNormalize(false); 156 }; 157 ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;}; 158 /*! 20 SphericalTransformServer() 21 { 22 fftIntfPtr_=new FFTPackServer; 23 fftIntfPtr_->setNormalize(false); 24 }; 25 ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;}; 26 27 /*! 159 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. 160 29 */ 161 void SetFFTServer(FFTServerInterface* srv=NULL) 162 { 163 if (fftIntfPtr_!=NULL) delete fftIntfPtr_; 164 fftIntfPtr_=srv; 165 fftIntfPtr_->setNormalize(false); 166 } 167 /*! synthesis of a temperature map from Alm coefficients */ 168 void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; 169 /*! synthesis of a polarization map from Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */ 170 void GenerateFromAlm(SphericalMap<T>& mapq, SphericalMap<T>& mapu, int_4 pixelSizeIndex, const Alm<T>& alme, const Alm<T>& almb) const; 30 void SetFFTServer(FFTServerInterface* srv=NULL) 31 { 32 if (fftIntfPtr_!=NULL) delete fftIntfPtr_; 33 fftIntfPtr_=srv; 34 fftIntfPtr_->setNormalize(false); 35 } 171 36 172 /*! synthesis of a temperature map from power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */ 173 void GenerateFromCl(SphericalMap<T>& sph, int_4 pixelSizeIndex, 37 void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; 38 void GenerateFromAlm(SphericalMap<T>& mapq, SphericalMap<T>& mapu, int_4 pixelSizeIndex, const Alm<T>& alme, const Alm<T>& almb) const; 39 40 void GenerateFromCl(SphericalMap<T>& sph, int_4 pixelSizeIndex, 174 41 const TVector<T>& Cl, const r_8 fwhm) const; 175 /*! synthesis of a polarization map from power spectra electric-Cl and magnetic-Cl (Alm's are generated randomly, following a gaussian distribution). 176 \param fwhm FWHM in arcmin for random generation of Alm's (eg. 5) 177 178 */ 179 void GenerateFromCl(SphericalMap<T>& sphq, SphericalMap<T>& sphu, 42 void GenerateFromCl(SphericalMap<T>& sphq, SphericalMap<T>& sphu, 180 43 int_4 pixelSizeIndex, 181 44 const TVector<T>& Cle, const TVector<T>& Clb, 182 45 const r_8 fwhm) const; 183 /*!return the Alm coefficients from analysis of a temperature map.184 185 \param<nlmax> : maximum value of the l index186 187 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.188 */189 46 190 47 191 Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const; 192 /*!analysis of a polarization map into Alm coefficients. 48 Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const; 193 49 194 The spheres \c mapq and \c mapu contain respectively the Stokes parameters. 195 196 \c a2lme and \c a2lmb will receive respectively electric and magnetic Alm's 197 nlmax : maximum value of the l index 198 199 \c cos_theta_cut : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut. 200 */ 201 202 void DecomposeToAlm(const SphericalMap<T>& mapq, const SphericalMap<T>& mapu, 50 void DecomposeToAlm(const SphericalMap<T>& mapq, const SphericalMap<T>& mapu, 203 51 Alm<T>& a2lme, Alm<T>& a2lmb, 204 52 int_4 nlmax, r_8 cos_theta_cut) const; 205 53 206 /*!return power spectrum from analysis of a temperature map. 207 208 \param<nlmax> : maximum value of the l index 209 210 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut. 211 */ 212 TVector<T> DecomposeToCl(const SphericalMap<T>& sph, 54 TVector<T> DecomposeToCl(const SphericalMap<T>& sph, 213 55 int_4 nlmax, r_8 cos_theta_cut) const; 214 56 215 57 216 private: 217 /*! return a vector with nph elements which are sums :\f$\sum_{m=-mmax}^{mmax}b_m(\theta)e^{im\varphi}\f$ for nph values of \f$\varphi\f$ regularly distributed in \f$[0,\pi]\f$ ( calculated by FFT) 218 219 The object b_m (\f$b_m\f$) of the class Bm is a special vector which index goes from -mmax to mmax. 220 */ 221 TVector< complex<T> > fourierSynthesisFromB(const Bm<complex<T> >& b_m, 58 private: 59 TVector< complex<T> > fourierSynthesisFromB(const Bm<complex<T> >& b_m, 222 60 int_4 nph, r_8 phi0) const; 223 /*! same as fourierSynthesisFromB, but return a real vector, taking into account the fact that b(-m) is conjugate of b(m) */ 224 TVector<T> RfourierSynthesisFromB(const Bm<complex<T> >& b_m, 61 TVector<T> RfourierSynthesisFromB(const Bm<complex<T> >& b_m, 225 62 int_4 nph, r_8 phi0) const; 226 63 227 /*! return a vector with mmax elements which are sums : 228 \f$\sum_{k=0}^{nphi}datain(\theta,\varphi_k)e^{im\varphi_k}\f$ for (mmax+1) values of \f$m\f$ from 0 to mmax. 229 */ 230 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 64 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 231 65 const TVector<complex<T> > datain, 232 66 r_8 phi0) const; 233 /* same as previous one, but with a "datain" which is real (not complex) */ 234 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 67 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 235 68 const TVector<T> datain, 236 69 r_8 phi0) const; 237 70 238 /*! 239 Compute polarized Alm's as : 240 \f[ 241 a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(\,_{w}\lambda_l^m\tilde{Q}-i\,_{x}\lambda_l^m\tilde{U}\right)} 242 \f] 243 \f[ 244 a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(i\,_{x}\lambda_l^m\tilde{Q}+\,_{w}\lambda_l^m\tilde{U}\right)} 245 \f] 246 247 where \f$\tilde{Q}\f$ and \f$\tilde{U}\f$ are C-coefficients computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters. 248 249 \f$\omega_{pix}\f$ are solid angle of each pixel. 250 251 dataq, datau : Stokes parameters. 252 253 */ 254 void almFromWX(int_4 nlmax, int_4 nmmax, r_8 phi0, 71 void almFromWX(int_4 nlmax, int_4 nmmax, r_8 phi0, 255 72 r_8 domega, r_8 theta, 256 73 const TVector<T>& dataq, const TVector<T>& datau, 257 74 Alm<T>& alme, Alm<T>& almb) const; 258 /*! 259 Compute polarized Alm's as : 260 \f[ 261 a_{lm}^E=-\frac{1}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}+\,_{-}\lambda_l^m\tilde{P^-}\right)} 262 \f] 263 \f[ 264 a_{lm}^B=\frac{i}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}-\,_{-}\lambda_l^m\tilde{P^-}\right)} 265 \f] 266 267 where \f$\tilde{P^{\pm}}=\tilde{Q}\pm\tilde{U}\f$ computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters,\f$Q\f$ and \f$U\f$ . 268 269 \f$\omega_{pix}\f$ are solid angle of each pixel. 270 271 dataq, datau : Stokes parameters. 272 273 */ 274 void almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax, 75 void almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax, 275 76 r_8 phi0, r_8 domega, r_8 theta, 276 77 const TVector<T>& dataq, const TVector<T>& datau, 277 78 Alm<T>& alme, Alm<T>& almb) const; 278 79 279 /*! synthesis of Stokes parameters following formulae : 280 281 \f[ 282 Q=\sum_{m=-mmax}^{mmax}b_m^qe^{im\varphi} 283 \f] 284 \f[ 285 U=\sum_{m=-mmax}^{mmax}b_m^ue^{im\varphi} 286 \f] 287 288 computed by FFT (method fourierSynthesisFromB called by the present one) 289 290 with : 291 292 \f[ 293 b_m^q=-\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(\,_{w}\lambda_l^ma_{lm}^E-i\,_{x}\lambda_l^ma_{lm}^B\right) } 294 \f] 295 \f[ 296 b_m^u=\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(i\,_{x}\lambda_l^ma_{lm}^E+\,_{w}\lambda_l^ma_{lm}^B\right) } 297 \f] 298 */ 299 void mapFromWX(int_4 nlmax, int_4 nmmax, 80 void mapFromWX(int_4 nlmax, int_4 nmmax, 300 81 SphericalMap<T>& mapq, SphericalMap<T>& mapu, 301 82 const Alm<T>& alme, const Alm<T>& almb) const; 302 83 303 /*! synthesis of polarizations following formulae :304 84 305 \f[ 306 P^+ = \sum_{m=-mmax}^{mmax} {b_m^+e^{im\varphi} } 307 \f] 308 \f[ 309 P^- = \sum_{m=-mmax}^{mmax} {b_m^-e^{im\varphi} } 310 \f] 311 312 computed by FFT (method fourierSynthesisFromB called by the present one) 313 314 with : 315 316 \f[ 317 b_m^+=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E+ia_{lm}^B \right) } 318 \f] 319 \f[ 320 b_m^-=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E-ia_{lm}^B \right) } 321 \f] 322 */ 323 324 void mapFromPM(int_4 nlmax, int_4 nmmax, 85 void mapFromPM(int_4 nlmax, int_4 nmmax, 325 86 SphericalMap<T>& mapq, SphericalMap<T>& mapu, 326 87 const Alm<T>& alme, const Alm<T>& almb) const;
Note:
See TracChangeset
for help on using the changeset viewer.