Changeset 515 in Sophya for trunk/SophyaLib
- Timestamp:
- Oct 26, 1999, 2:44:44 PM (26 years ago)
- Location:
- trunk/SophyaLib
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/dvlist.cc
r514 r515 64 64 65 65 /* --Methode-- */ 66 DVList::DVList( DVList& dvl)66 DVList::DVList(const DVList& dvl) 67 67 { 68 68 Merge(dvl); -
trunk/SophyaLib/NTools/dvlist.h
r490 r515 68 68 69 69 DVList(); 70 DVList( DVList&);70 DVList(const DVList&); 71 71 DVList(char *flnm); 72 72 -
trunk/SophyaLib/NTools/fftserver.cc
r514 r515 92 92 checkint_cfft(l); 93 93 float* foo = new float[2*l]; 94 for (int i=0;i<l;i++){ 94 int i; 95 for (i=0;i<l;i++){ 95 96 foo[2*i]=inout[i].real(); 96 97 foo[2*i+1]=inout[i].imag(); … … 98 99 cfftf_(&l, foo, ws_cfft); 99 100 inout[0]=complex<float> (foo[0],foo[1]); 100 for (i nt i=1;i<l;i++) inout[l-i]= complex<float> (foo[2*i], foo[2*i+1]);101 for (i=1;i<l;i++) inout[l-i]= complex<float> (foo[2*i], foo[2*i+1]); 101 102 delete[] foo; 102 103 } … … 106 107 checkint_cdfft(l); 107 108 double* foo=new double[2*l]; 108 for (int i=0;i<l;i++){ 109 int i; 110 for (i=0;i<l;i++){ 109 111 foo[2*i]=inout[i].real(); 110 112 foo[2*i+1]=inout[i].imag(); … … 112 114 cdfftf_(&l, foo, ws_cdfft); 113 115 inout[0]=complex<double> (foo[0],foo[1]); 114 for (i nt i=1;i<l;i++) {116 for (i=1;i<l;i++) { 115 117 inout[l-i]= complex<double> (foo[2*i],foo[2*i+1]); 116 118 } … … 134 136 checkint_cfft(l); 135 137 float* foo = new float[2*l]; 136 for (int i=0;i<l;i++){ 138 int i; 139 for (i=0;i<l;i++){ 137 140 foo[2*i]=inout[i].real(); 138 141 foo[2*i+1]=inout[i].imag(); 139 142 } 140 143 cfftf_(&l, foo, ws_cfft); 141 for (i nt i=0;i<l;i++) inout[i]=complex<float> (foo[2*i],foo[2*i+1]);144 for (i=0;i<l;i++) inout[i]=complex<float> (foo[2*i],foo[2*i+1]); 142 145 delete[] foo; 143 146 } … … 147 150 checkint_cdfft(l); 148 151 double* foo = new double[2*l]; 149 for (int i=0;i<l;i++){ 152 int i; 153 for (i=0;i<l;i++){ 150 154 foo[2*i]=inout[i].real(); 151 155 foo[2*i+1]=inout[i].imag(); 152 156 } 153 157 cdfftf_(&l, foo, ws_cdfft); 154 for (i nt i=0;i<l;i++) inout[i]=complex<double> (foo[2*i],foo[2*i+1]);158 for (i=0;i<l;i++) inout[i]=complex<double> (foo[2*i],foo[2*i+1]); 155 159 delete[] foo; 156 160 } -
trunk/SophyaLib/Samba/scan.cc
r276 r515 116 116 PhiZero_ = s. PhiZero_; 117 117 sPix_=new r_8[ NmaxPts_]; 118 for (int_4 k=0; k<NmaxPts_; k++) sPix_[k]=s.sPix_[k]; 119 for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k]; 118 int_4 k; 119 for (k=0; k<NmaxPts_; k++) sPix_[k]=s.sPix_[k]; 120 for (k=0; k<9; k++) Rota_[k]=s. Rota_[k]; 120 121 } 121 122 //++ -
trunk/SophyaLib/Samba/spheregorski.cc
r487 r515 1 1 #include "spheregorski.h" 2 2 #include "strutil.h" 3 #include <math.h> 3 4 #include <complex> 4 5 #include "piocmplx.h" … … 886 887 iy_hi = iy/128; 887 888 ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low)); 888 ipf = ipf / pow(ns_max/nside,2);// ! in {0, nside**2 - 1} 889 return ( ipf + face_num*pow(nside,2));// ! in {0, 12*nside**2 - 1} 889 // ipf = ipf / pow(ns_max/nside,2.);// ! in {0, nside**2 - 1} 890 // return ( ipf + face_num*pow(nside,2));// ! in {0, 12*nside**2 - 1} 891 // $CHECK$ Reza 25/10/99 , pow remplace par * 892 // ipf = ipf / ((ns_max/nside)*(ns_max/nside)); // ! in {0, nside**2 - 1} 893 // return ( ipf + face_num*(nside*nside);// ! in {0, 12*nside**2 - 1} 890 894 } 891 895 -
trunk/SophyaLib/Samba/spherethetaphi.cc
r487 r515 4 4 #include "piocmplx.h" 5 5 #include <iostream.h> 6 7 8 //*************************************************************** 9 //++ 10 // Class SphereThetaPhi 11 // 12 // 13 // include spherethetaphi.h nbmath.h 14 // 15 // Découpage de la sphère selon theta et phi, chaque 16 // hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du 17 // beurre), chacune des m bandes de theta ainsi définies étant découpée par 18 // des méridiens équirepartis, ce découpage étant fait de sorte que tous 19 // les pixels aient la même surface et soient le plus carré possible. 20 // On commence par découper l'hémisphère de z positif en partant du pôle et 21 // en allant vers l'équateur. Le premier pixel est la calotte polaire, 22 // il est circulaire et centré sur theta=0. 23 //-- 24 //++ 25 // 26 // Links Parents 27 // 28 // SphericalMap 29 //-- 30 //++ 31 // 32 // Links Descendants 33 // 34 // 35 //-- 36 37 /* --Methode-- */ 38 //++ 39 // Titre Constructeurs 40 //-- 41 //++ 42 43 template <class T> 44 SphereThetaPhi<T>::SphereThetaPhi() 45 46 //-- 47 { 48 InitNul(); 49 } 50 51 52 /* --Methode-- */ 53 54 //++ 55 template <class T> 56 SphereThetaPhi<T>::SphereThetaPhi(int m) 57 58 // Constructeur : m est le nombre de bandes en theta sur un hémisphère 59 // (la calotte constituant la premiere bande). 60 // pet est le nombre de pixels (pétales) de la bande en contact avec la 61 // calotte polaire. Pour l'instant pet est inopérant! 62 //-- 63 { 64 InitNul(); 65 Pixelize(m); 66 } 67 68 template <class T> 69 SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s) 70 { 71 if(s.mInfo_) mInfo_= new DVList(*s.mInfo_); 72 NTheta_= s.NTheta_; 73 NPix_ = s.NPix_; 74 NPhi_ = new int[NTheta_]; 75 Theta_ = new double[NTheta_+1]; 76 TNphi_ = new int[NTheta_+1]; 77 for(int k = 0; k < NTheta_; k++) 78 { 79 NPhi_[k] = s.NPhi_[k]; 80 Theta_[k]= s.Theta_[k]; 81 TNphi_[k]= s.TNphi_[k]; 82 } 83 Theta_[NTheta_]= s.Theta_[NTheta_]; 84 TNphi_[NTheta_]= s.TNphi_[NTheta_]; 85 Omega_ = s.Omega_; 86 pixels_= s.pixels_; 87 } 88 89 //++ 90 // Titre Destructeur 91 //-- 92 //++ 93 template <class T> 94 SphereThetaPhi<T>::~SphereThetaPhi() 95 96 //-- 97 { 98 Clear(); 99 } 100 101 //++ 102 // Titre Méthodes 103 //-- 104 template <class T> 105 void SphereThetaPhi<T>::InitNul() 106 // 107 // initialise à zéro les variables de classe pointeurs 108 { 109 NTheta_= 0; 110 NPix_ = 0; 111 Theta_ = NULL; 112 NPhi_ = NULL; 113 TNphi_ = NULL; 114 pixels_.Reset(); 115 } 116 117 /* --Methode-- */ 118 template <class T> 119 void SphereThetaPhi<T>::Clear() 120 121 { 122 if(Theta_) delete[] Theta_; 123 if(NPhi_ ) delete[] NPhi_; 124 if(TNphi_) delete[] TNphi_; 125 InitNul(); 126 } 127 128 //++ 129 template <class T> 130 void SphereThetaPhi<T>::Resize(int m) 131 // re-pixelize the sphere 132 //-- 133 { 134 Clear(); 135 Pixelize(m); 136 } 137 138 /* --Methode-- */ 139 //++ 140 template <class T> 141 int SphereThetaPhi<T>::NbPixels() const 142 143 // Retourne le nombre de pixels du découpage 144 //-- 145 { 146 return(NPix_); 147 } 148 149 /* --Methode-- */ 150 //++ 151 template <class T> 152 T& SphereThetaPhi<T>::PixVal(int k) 153 154 // Retourne la valeur du contenu du pixel d'indice k 155 //-- 156 { 157 if((k < 0) || (k >= NPix_)) 158 { 159 //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 160 cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl; 161 THROW(rangeCheckErr); 162 } 163 return pixels_(k); 164 } 165 166 //++ 167 template <class T> 168 T const& SphereThetaPhi<T>::PixVal(int k) const 169 170 // Retourne la valeur du contenu du pixel d'indice k 171 //-- 172 { 173 if((k < 0) || (k >= NPix_)) 174 { 175 cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl; 176 //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 177 throw "SphereThetaPhi::PIxVal Pixel index out of range"; 178 } 179 return *(pixels_.Data()+k); 180 } 181 182 /* --Methode-- */ 183 //++ 184 template <class T> 185 int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const 186 187 // Retourne l'indice du pixel vers lequel pointe une direction définie par 188 // ses coordonnées sphériques 189 //-- 190 { 191 double dphi; 192 int i,j,k; 193 bool fgzn = false; 194 195 if((theta > Pi) || (theta < 0.)) return(-1); 196 if((phi < 0.) || (phi > DeuxPi)) return(-1); 197 if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;} 198 199 // La bande d'indice kt est limitée par les valeurs de theta contenues dans 200 // Theta_[kt] et Theta_[kt+1] 201 for( i=1; i< NTheta_; i++ ) 202 if( theta < Theta_[i] ) break; 203 204 dphi= DeuxPi/(double)NPhi_[i-1]; 205 206 if (fgzn) k= NPix_-TNphi_[i]+(int)(phi/dphi); 207 else k= TNphi_[i-1]+(int)(phi/dphi); 208 return(k); 209 } 210 211 /* --Methode-- */ 212 //++ 213 template <class T> 214 void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const 215 216 // Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k 217 //-- 218 { 219 int i; 220 bool fgzn = false; 221 222 if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; } 223 if( k >= NPix_/2) {fgzn = true; k = NPix_-1-k;} 224 225 // recupère l'indice i de la tranche qui contient le pixel k 226 for( i=0; i< NTheta_; i++ ) 227 if( k < TNphi_[i+1] ) break; 228 229 // angle theta 230 theta= 0.5*(Theta_[i]+Theta_[i+1]); 231 if (fgzn) theta= Pi-theta; 232 233 // angle phi 234 k -= TNphi_[i]; 235 phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5); 236 if (fgzn) phi= DeuxPi-phi; 237 } 238 239 //++ 240 template <class T> 241 double SphereThetaPhi<T>::PixSolAngle(int dummy) const 242 243 // Pixel Solid angle (steradians) 244 // All the pixels have the same solid angle. The dummy argument is 245 // for compatibility with eventual pixelizations which would not 246 // fulfil this requirement. 247 //-- 248 { 249 return Omega_; 250 } 251 252 /* --Methode-- */ 253 //++ 254 template <class T> 255 void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax) 256 257 // Retourne les valeurs de theta et phi limitant le pixel d'indice k 258 //-- 259 { 260 int j; 261 double dphi; 262 bool fgzn= false; 263 264 if((k < 0) || (k >= NPix_)) { 265 tetMin= -99999.; 266 phiMin= -99999.; 267 tetMax= -99999.; 268 phiMax= -99999.; 269 return; 270 } 271 272 // si on se trouve dans l'hémisphère Sud 273 if(k >= NPix_/2) { 274 fgzn= true; 275 k= NPix_-1-k; 276 } 277 278 // on recupere l'indice i de la tranche qui contient le pixel k 279 int i; 280 for( i=0; i< NTheta_; i++ ) 281 if(k < TNphi_[i+1]) break; 282 283 // valeurs limites de theta dans l'hemisphere Nord 284 tetMin= Theta_[i]; 285 tetMax= Theta_[i+1]; 286 // valeurs limites de theta dans l'hemisphere Sud 287 if (fgzn) { 288 tetMin= Pi-Theta_[i+1]; 289 tetMax= Pi-Theta_[i]; 290 } 291 292 // pixel correspondant dans l'hemisphere Nord 293 if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1; 294 295 // indice j de discretisation ( phi= j*dphi ) 296 j= k-TNphi_[i]; 297 dphi= DeuxPi/(double)NPhi_[i]; 298 299 // valeurs limites de phi 300 phiMin= dphi*(double)(j); 301 phiMax= dphi*(double)(j+1); 302 return; 303 } 304 305 /* --Methode-- */ 306 //++ 307 template <class T> 308 int SphereThetaPhi<T>::NbThetaSlices() const 309 310 // Retourne le nombre de tranches en theta sur la sphere 311 //-- 312 { 313 int nbslices; 314 nbslices= 2*NTheta_; 315 return(nbslices); 316 } 317 318 /* --Methode-- */ 319 //++ 320 template <class T> 321 int SphereThetaPhi<T>::NPhi(int kt) const 322 323 // Retourne le nombre de pixels en phi de la tranche kt 324 //-- 325 { 326 int nbpix; 327 // verification 328 if((kt < 0) || (kt >= 2*NTheta_)) return(-1); 329 330 // si on se trouve dans l'hemisphere Sud 331 if(kt >= NTheta_) { 332 kt= 2*NTheta_-1-kt; 333 } 334 335 // nombre de pixels 336 nbpix= NPhi_[kt]; 337 return(nbpix); 338 } 339 340 341 /* --Methode-- */ 342 //++ 343 template <class T> 344 void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax) 345 346 // Retourne les valeurs de theta limitant la tranche kt 347 //-- 348 { 349 bool fgzn= false; 350 // verification 351 if( (kt< 0) || (kt>= 2*NTheta_) ) { 352 tetMin= -99999.; 353 tetMax= -99999.; 354 return; 355 } 356 357 // si on se trouve dans l'hemisphere Sud 358 if( kt >= NTheta_ ) { 359 fgzn= true; 360 kt= 2*NTheta_-1-kt; 361 } 362 363 // valeurs limites de theta dans l'hemisphere Nord 364 tetMin= Theta_[kt]; 365 tetMax= Theta_[kt+1]; 366 // valeurs limites de theta dans l'hemisphere Sud 367 if (fgzn) { 368 tetMin= Pi-Theta_[kt+1]; 369 tetMax= Pi-Theta_[kt]; 370 } 371 } 372 373 /* --Methode-- */ 374 //++ 375 template <class T> 376 void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax) 377 378 // Retourne les valeurs de phi limitant le pixel jp de la tranche kt 379 //-- 380 { 381 // verification 382 if((kt < 0) || (kt >= 2*NTheta_)) { 383 phiMin= -99999.; 384 phiMax= -99999.; 385 return; 386 } 387 388 // si on se trouve dans l'hemisphere Sud 389 if(kt >= NTheta_) kt= 2*NTheta_-1-kt; 390 391 // verifie si la tranche kt contient au moins jp pixels 392 if( (jp< 0) || (jp >= NPhi_[kt]) ) { 393 phiMin= -88888.; 394 phiMax= -88888.; 395 return; 396 } 397 398 double dphi= DeuxPi/(double)NPhi_[kt]; 399 phiMin= dphi*(double)(jp); 400 phiMax= dphi*(double)(jp+1); 401 return; 402 } 403 404 /* --Methode-- */ 405 //++ 406 template <class T> 407 int SphereThetaPhi<T>::Index(int kt,int jp) const 408 409 // Retourne l'indice du pixel d'indice jp dans la tranche kt 410 //-- 411 { 412 int k; 413 bool fgzn= false; 414 415 // si on se trouve dans l'hemisphere Sud 416 if(kt >= NTheta_) { 417 fgzn= true; 418 kt= 2*NTheta_-1-kt; 419 } 420 421 // si la tranche kt contient au moins i pixels 422 if( (jp>=0) && (jp<NPhi_[kt]) ) 423 { 424 // dans l'hemisphere Sud 425 if (fgzn) k= NPix_-TNphi_[kt+1]+jp; 426 // dans l'hemisphere Nord 427 else k= TNphi_[kt]+jp; 428 } 429 else 430 { 431 k= 9999; 432 printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp); 433 } 434 return(k); 435 } 436 437 /* --Methode-- */ 438 //++ 439 template <class T> 440 void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp) 441 442 // Retourne les indices kt et jp du pixel d'indice k 443 //-- 444 { 445 bool fgzn= false; 446 // si on se trouve dans l'hemisphere Sud 447 if(k >= NPix_/2) 448 { 449 fgzn= true; 450 k= NPix_-1-k; 451 } 452 453 // on recupere l'indice kt de la tranche qui contient le pixel k 454 int i; 455 for(i = 0; i < NTheta_; i++) 456 if(k < TNphi_[i+1]) break; 457 458 // indice kt de tranche 459 if (fgzn) kt= 2*NTheta_-1-i; 460 else kt= i; 461 462 // indice jp de pixel 463 if (fgzn) jp= TNphi_[i+1]-k-1; 464 else jp= k-TNphi_[i]; 465 } 466 //++ 467 template <class T> 468 void SphereThetaPhi<T>::Pixelize(int m) 469 470 // effectue le découpage en pixels (m et pet ont la même signification 471 // que pour le constructeur) 472 // 473 // Chaque bande de theta sera découpée en partant de phi=0 ... 474 // L'autre hémisphère est parcourue dans le même sens en phi et de 475 // l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus 476 // proche de l'équateur a z>0 est celui de plus petit phi de la bande la 477 // plus proche de l'equateur a z<0). 478 //-- 479 { 480 int ntotpix,i,j; 481 482 // Decodage et controle des arguments d'appel : 483 // au moins 2 et au plus 16384 bandes d'un hemisphere en theta 484 if (m < 2) m = 2; 485 if (m > 16384) m = 16384; 486 487 // On memorise les arguments d'appel 488 NTheta_ = m; 489 490 // On commence par decouper l'hemisphere z>0. 491 // Creation des vecteurs contenant : 492 // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...) 493 Theta_= new double[m+1]; 494 495 // Le nombre de pixels en phi de chacune des bandes en theta 496 NPhi_ = new int[m]; 497 498 // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une 499 // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere) 500 TNphi_= new int[m+1]; 501 502 // Calcul du nombre total de pixels dans chaque bande pour optimiser 503 // le rapport largeur/hauteur des pixels 504 505 //calotte polaire 506 TNphi_[0]= 0; 507 NPhi_[0] = 1; 508 509 //bandes jusqu'a l'equateur 510 for(j = 1; j < m; j++) 511 { 512 TNphi_[j]= TNphi_[j-1]+NPhi_[j-1]; 513 NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ; 514 } 515 516 // Nombre total de pixels sur l'hemisphere 517 ntotpix = TNphi_[m-1]+NPhi_[m-1]; 518 TNphi_[m]= ntotpix; 519 // et sur la sphere entiere 520 NPix_= 2*ntotpix; 521 522 // Creation et initialisation du vecteur des contenus des pixels 523 pixels_.ReSize(NPix_); 524 pixels_.Reset(); 525 526 // Determination des limites des bandes en theta : 527 // omeg est l'angle solide couvert par chaque pixel, 528 // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg 529 // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide 530 //de la 531 // calotte allant du pole a la limite haute de la bande kt vaut 532 // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg... 533 534 double omeg2pi= 1./(double)ntotpix; 535 Omega_ = omeg2pi*DeuxPi; 536 537 for(j=0; j <= m; j++) 538 { 539 Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi); 540 } 541 } 542 543 //++ 544 template <class T> 545 void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const 546 547 // Retourne, pour la tranche en theta d'indice 'index' le theta 548 // correspondant, un vecteur (Peida) contenant les phi des pixels de 549 // la tranche, un vecteur (Peida) contenant les valeurs de pixel 550 // correspondantes 551 //-- 552 553 { 554 cout << "entree GetThetaSlice, couche no " << index << endl; 555 556 if(index < 0 || index > NbThetaSlices()) 557 { 558 // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 559 cout << " SphereThetaPhi::GetThetaSlice : exceptions a mettre en place" <<endl; 560 THROW(rangeCheckErr); 561 } 562 563 int iring= Index(index,0); 564 int bid = this->NPhi(index); 565 int lring = bid; 566 cout << " iring= " << iring << " lring= " << lring << endl; 567 568 phi.ReSize(lring); 569 value.ReSize(lring); 570 double Te= 0.; 571 double Fi= 0.; 572 for(int kk = 0; kk < lring; kk++) 573 { 574 PixThetaPhi(kk+iring,Te,Fi); 575 phi(kk)= Fi; 576 value(kk)= PixVal(kk+iring); 577 } 578 theta= Te; 579 } 580 581 template <class T> 582 void SphereThetaPhi<T>::setmNPhi(int* array, int m) 583 //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta 584 //-- 585 { 586 NPhi_= new int[m]; 587 for(int k = 0; k < m; k++) NPhi_[k]= array[k]; 588 } 589 590 template <class T> 591 void SphereThetaPhi<T>::setmTNphi(int* array, int m) 592 //remplit le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes 593 //-- 594 { 595 TNphi_= new int[m]; 596 for(int k = 0; k < m; k++) TNphi_[k]= array[k]; 597 } 598 599 template <class T> 600 void SphereThetaPhi<T>::setmTheta(double* array, int m) 601 //remplit le tableau contenant les valeurs limites de theta 602 //-- 603 { 604 Theta_= new double[m]; 605 for(int k = 0; k < m; k++) Theta_[k]= array[k]; 606 } 607 608 template <class T> 609 void SphereThetaPhi<T>::setDataBlock(T* data, int m) 610 // remplit le vecteur des contenus des pixels 611 { 612 pixels_.FillFrom(m,data); 613 } 614 615 template <class T> 616 void SphereThetaPhi<T>::print(ostream& os) const 617 { 618 if(mInfo_) os << " DVList Info= " << *mInfo_ << endl; 619 // 620 os << " NTheta_= " << NTheta_ << endl; 621 os << " NPix_ = " << NPix_ << endl; 622 os << " Omega_ = " << Omega_ << endl; 623 624 os << " contenu de NPhi_ : "; 625 int i; 626 for(i=0; i < NTheta_; i++) 627 { 628 if(i%5 == 0) os << endl; 629 os << NPhi_[i] <<", "; 630 } 631 os << endl; 632 633 os << " contenu de Theta_ : "; 634 for(i=0; i < NTheta_+1; i++) 635 { 636 if(i%5 == 0) os << endl; 637 os << Theta_[i] <<", "; 638 } 639 os << endl; 640 641 os << " contenu de TNphi_ : "; 642 for(i=0; i < NTheta_+1; i++) 643 { 644 if(i%5 == 0) os << endl; 645 os << TNphi_[i] <<", "; 646 } 647 os << endl; 648 649 os << " contenu de pixels : "; 650 for(i=0; i < NPix_; i++) 651 { 652 if(i%5 == 0) os << endl; 653 os << pixels_(i) <<", "; 654 } 655 os << endl; 656 } 657 658 /////////////////////////////////////////////////////////// 659 // -------------------------------------------------------- 660 // Les objets delegues pour la gestion de persistance 661 // -------------------------------------------------------- 662 ////////////////////////////////////////////////////////// 663 664 template <class T> 665 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi() 666 { 667 dobj= new SphereThetaPhi<T>; 668 ownobj= true; 669 } 670 671 template <class T> 672 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename) 673 { 674 dobj= new SphereThetaPhi<T>; 675 ownobj= true; 676 Read(filename); 677 } 678 679 template <class T> 680 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj) 681 { 682 dobj= new SphereThetaPhi<T>(obj); 683 ownobj= true; 684 } 685 686 template <class T> 687 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj) 688 { 689 dobj= obj; 690 ownobj= false; 691 } 692 693 template <class T> 694 FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi() 695 { 696 if (ownobj && dobj) delete dobj; 697 } 698 699 template <class T> 700 AnyDataObj* FIO_SphereThetaPhi<T>::DataObj() 701 { 702 return(dobj); 703 } 704 705 template <class T> 706 void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is) 707 { 708 cout << " FIO_SphereThetaPhi:: ReadSelf " << endl; 709 710 if(dobj == NULL) 711 { 712 dobj= new SphereThetaPhi<T>; 713 } 714 715 // Pour savoir s'il y avait un DVList Info associe 716 char strg[256]; 717 is.GetLine(strg, 255); 718 bool hadinfo= false; 719 if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo= true; 720 if(hadinfo) 721 { // Lecture eventuelle du DVList Info 722 is >> dobj->Info(); 723 } 724 725 int mNTheta; 726 is.GetI4(mNTheta); 727 dobj->setSizeIndex(mNTheta); 728 729 int mNPix; 730 is.GetI4(mNPix); 731 dobj->setNbPixels(mNPix); 732 733 double mOmeg; 734 is.GetR8(mOmeg); 735 dobj->setPixSolAngle(mOmeg); 736 737 int* mNphi= new int[mNTheta]; 738 is.GetI4s(mNphi, mNTheta); 739 dobj->setmNPhi(mNphi, mNTheta); 740 delete [] mNphi; 741 742 int* mTNphi= new int[mNTheta+1]; 743 is.GetI4s(mTNphi, mNTheta+1); 744 dobj->setmTNphi(mTNphi, mNTheta+1); 745 delete [] mTNphi; 746 747 double* mTheta= new double[mNTheta+1]; 748 is.GetR8s(mTheta, mNTheta+1); 749 dobj->setmTheta(mTheta, mNTheta+1); 750 delete [] mTheta; 751 752 T* mPixels= new T[mNPix]; 753 PIOSReadArray(is, mPixels, mNPix); 754 dobj->setDataBlock(mPixels, mNPix); 755 delete [] mPixels; 756 } 757 758 template <class T> 759 void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const 760 { 761 cout << " FIO_SphereThetaPhi:: WriteSelf " << endl; 762 763 if(dobj == NULL) 764 { 765 cout << " WriteSelf:: dobj= null " << endl; 766 return; 767 } 768 769 char strg[256]; 770 int mNTheta= dobj->SizeIndex(); 771 int mNPix = dobj->NbPixels(); 772 773 if(dobj->ptrInfo()) 774 { 775 sprintf(strg,"SphereThetaPhi: NSlices=%6d NPix=%9d HasInfo",mNTheta,mNPix); 776 os.PutLine(strg); 777 os << dobj->Info(); 778 } 779 else 780 { 781 sprintf(strg,"SphereThetaPhi: NSlices=%6d NPix=%9d ",mNTheta,mNPix); 782 os.PutLine(strg); 783 } 784 785 os.PutI4(mNTheta); 786 os.PutI4(mNPix); 787 os.PutR8(dobj->PixSolAngle(0)); 788 os.PutI4s(dobj->getmNPhi() , mNTheta); 789 os.PutI4s(dobj->getmTNphi(), mNTheta+1); 790 os.PutR8s(dobj->getmTheta(), mNTheta+1); 791 //os.Put((dobj->getDataBlock())->Data(), mNPix); 792 PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix); 793 } 6 794 7 795 #ifdef __CXX_PRAGMA_TEMPLATES__ … … 25 813 template class FIO_SphereThetaPhi< complex<double> >; 26 814 #endif 27 28 //***************************************************************29 //++30 // Class SphereThetaPhi31 //32 //33 // include spherethetaphi.h nbmath.h34 //35 // Découpage de la sphère selon theta et phi, chaque36 // hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du37 // beurre), chacune des m bandes de theta ainsi définies étant découpée par38 // des méridiens équirepartis, ce découpage étant fait de sorte que tous39 // les pixels aient la même surface et soient le plus carré possible.40 // On commence par découper l'hémisphère de z positif en partant du pôle et41 // en allant vers l'équateur. Le premier pixel est la calotte polaire,42 // il est circulaire et centré sur theta=0.43 //--44 //++45 //46 // Links Parents47 //48 // SphericalMap49 //--50 //++51 //52 // Links Descendants53 //54 //55 //--56 57 /* --Methode-- */58 //++59 // Titre Constructeurs60 //--61 //++62 63 template <class T>64 SphereThetaPhi<T>::SphereThetaPhi()65 66 //--67 {68 InitNul();69 }70 71 72 /* --Methode-- */73 74 //++75 template <class T>76 SphereThetaPhi<T>::SphereThetaPhi(int m)77 78 // Constructeur : m est le nombre de bandes en theta sur un hémisphère79 // (la calotte constituant la premiere bande).80 // pet est le nombre de pixels (pétales) de la bande en contact avec la81 // calotte polaire. Pour l'instant pet est inopérant!82 //--83 {84 InitNul();85 Pixelize(m);86 }87 88 template <class T>89 SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)90 {91 if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);92 NTheta_= s.NTheta_;93 NPix_ = s.NPix_;94 NPhi_ = new int[NTheta_];95 Theta_ = new double[NTheta_+1];96 TNphi_ = new int[NTheta_+1];97 for(int k = 0; k < NTheta_; k++)98 {99 NPhi_[k] = s.NPhi_[k];100 Theta_[k]= s.Theta_[k];101 TNphi_[k]= s.TNphi_[k];102 }103 Theta_[NTheta_]= s.Theta_[NTheta_];104 TNphi_[NTheta_]= s.TNphi_[NTheta_];105 Omega_ = s.Omega_;106 pixels_= s.pixels_;107 }108 109 //++110 // Titre Destructeur111 //--112 //++113 template <class T>114 SphereThetaPhi<T>::~SphereThetaPhi()115 116 //--117 {118 Clear();119 }120 121 //++122 // Titre Méthodes123 //--124 template <class T>125 void SphereThetaPhi<T>::InitNul()126 //127 // initialise à zéro les variables de classe pointeurs128 {129 NTheta_= 0;130 NPix_ = 0;131 Theta_ = NULL;132 NPhi_ = NULL;133 TNphi_ = NULL;134 pixels_.Reset();135 }136 137 /* --Methode-- */138 template <class T>139 void SphereThetaPhi<T>::Clear()140 141 {142 if(Theta_) delete[] Theta_;143 if(NPhi_ ) delete[] NPhi_;144 if(TNphi_) delete[] TNphi_;145 InitNul();146 }147 148 //++149 template <class T>150 void SphereThetaPhi<T>::Resize(int m)151 // re-pixelize the sphere152 //--153 {154 Clear();155 Pixelize(m);156 }157 158 /* --Methode-- */159 //++160 template <class T>161 int SphereThetaPhi<T>::NbPixels() const162 163 // Retourne le nombre de pixels du découpage164 //--165 {166 return(NPix_);167 }168 169 /* --Methode-- */170 //++171 template <class T>172 T& SphereThetaPhi<T>::PixVal(int k)173 174 // Retourne la valeur du contenu du pixel d'indice k175 //--176 {177 if((k < 0) || (k >= NPix_))178 {179 //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));180 cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;181 THROW(rangeCheckErr);182 }183 return pixels_(k);184 }185 186 //++187 template <class T>188 T const& SphereThetaPhi<T>::PixVal(int k) const189 190 // Retourne la valeur du contenu du pixel d'indice k191 //--192 {193 if((k < 0) || (k >= NPix_))194 {195 cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;196 //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));197 throw "SphereThetaPhi::PIxVal Pixel index out of range";198 }199 return *(pixels_.Data()+k);200 }201 202 /* --Methode-- */203 //++204 template <class T>205 int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const206 207 // Retourne l'indice du pixel vers lequel pointe une direction définie par208 // ses coordonnées sphériques209 //--210 {211 double dphi;212 int i,j,k;213 bool fgzn = false;214 215 if((theta > Pi) || (theta < 0.)) return(-1);216 if((phi < 0.) || (phi > DeuxPi)) return(-1);217 if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}218 219 // La bande d'indice kt est limitée par les valeurs de theta contenues dans220 // Theta_[kt] et Theta_[kt+1]221 for( i=1; i< NTheta_; i++ )222 if( theta < Theta_[i] ) break;223 224 dphi= DeuxPi/(double)NPhi_[i-1];225 226 if (fgzn) k= NPix_-TNphi_[i]+(int)(phi/dphi);227 else k= TNphi_[i-1]+(int)(phi/dphi);228 return(k);229 }230 231 /* --Methode-- */232 //++233 template <class T>234 void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const235 236 // Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k237 //--238 {239 int i;240 bool fgzn = false;241 242 if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }243 if( k >= NPix_/2) {fgzn = true; k = NPix_-1-k;}244 245 // recupère l'indice i de la tranche qui contient le pixel k246 for( i=0; i< NTheta_; i++ )247 if( k < TNphi_[i+1] ) break;248 249 // angle theta250 theta= 0.5*(Theta_[i]+Theta_[i+1]);251 if (fgzn) theta= Pi-theta;252 253 // angle phi254 k -= TNphi_[i];255 phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);256 if (fgzn) phi= DeuxPi-phi;257 }258 259 //++260 template <class T>261 double SphereThetaPhi<T>::PixSolAngle(int dummy) const262 263 // Pixel Solid angle (steradians)264 // All the pixels have the same solid angle. The dummy argument is265 // for compatibility with eventual pixelizations which would not266 // fulfil this requirement.267 //--268 {269 return Omega_;270 }271 272 /* --Methode-- */273 //++274 template <class T>275 void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax)276 277 // Retourne les valeurs de theta et phi limitant le pixel d'indice k278 //--279 {280 int j;281 double dphi;282 bool fgzn= false;283 284 if((k < 0) || (k >= NPix_)) {285 tetMin= -99999.;286 phiMin= -99999.;287 tetMax= -99999.;288 phiMax= -99999.;289 return;290 }291 292 // si on se trouve dans l'hémisphère Sud293 if(k >= NPix_/2) {294 fgzn= true;295 k= NPix_-1-k;296 }297 298 // on recupere l'indice i de la tranche qui contient le pixel k299 int i;300 for( i=0; i< NTheta_; i++ )301 if(k < TNphi_[i+1]) break;302 303 // valeurs limites de theta dans l'hemisphere Nord304 tetMin= Theta_[i];305 tetMax= Theta_[i+1];306 // valeurs limites de theta dans l'hemisphere Sud307 if (fgzn) {308 tetMin= Pi-Theta_[i+1];309 tetMax= Pi-Theta_[i];310 }311 312 // pixel correspondant dans l'hemisphere Nord313 if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;314 315 // indice j de discretisation ( phi= j*dphi )316 j= k-TNphi_[i];317 dphi= DeuxPi/(double)NPhi_[i];318 319 // valeurs limites de phi320 phiMin= dphi*(double)(j);321 phiMax= dphi*(double)(j+1);322 return;323 }324 325 /* --Methode-- */326 //++327 template <class T>328 int SphereThetaPhi<T>::NbThetaSlices() const329 330 // Retourne le nombre de tranches en theta sur la sphere331 //--332 {333 int nbslices;334 nbslices= 2*NTheta_;335 return(nbslices);336 }337 338 /* --Methode-- */339 //++340 template <class T>341 int SphereThetaPhi<T>::NPhi(int kt) const342 343 // Retourne le nombre de pixels en phi de la tranche kt344 //--345 {346 int nbpix;347 // verification348 if((kt < 0) || (kt >= 2*NTheta_)) return(-1);349 350 // si on se trouve dans l'hemisphere Sud351 if(kt >= NTheta_) {352 kt= 2*NTheta_-1-kt;353 }354 355 // nombre de pixels356 nbpix= NPhi_[kt];357 return(nbpix);358 }359 360 361 /* --Methode-- */362 //++363 template <class T>364 void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax)365 366 // Retourne les valeurs de theta limitant la tranche kt367 //--368 {369 bool fgzn= false;370 // verification371 if( (kt< 0) || (kt>= 2*NTheta_) ) {372 tetMin= -99999.;373 tetMax= -99999.;374 return;375 }376 377 // si on se trouve dans l'hemisphere Sud378 if( kt >= NTheta_ ) {379 fgzn= true;380 kt= 2*NTheta_-1-kt;381 }382 383 // valeurs limites de theta dans l'hemisphere Nord384 tetMin= Theta_[kt];385 tetMax= Theta_[kt+1];386 // valeurs limites de theta dans l'hemisphere Sud387 if (fgzn) {388 tetMin= Pi-Theta_[kt+1];389 tetMax= Pi-Theta_[kt];390 }391 }392 393 /* --Methode-- */394 //++395 template <class T>396 void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax)397 398 // Retourne les valeurs de phi limitant le pixel jp de la tranche kt399 //--400 {401 // verification402 if((kt < 0) || (kt >= 2*NTheta_)) {403 phiMin= -99999.;404 phiMax= -99999.;405 return;406 }407 408 // si on se trouve dans l'hemisphere Sud409 if(kt >= NTheta_) kt= 2*NTheta_-1-kt;410 411 // verifie si la tranche kt contient au moins jp pixels412 if( (jp< 0) || (jp >= NPhi_[kt]) ) {413 phiMin= -88888.;414 phiMax= -88888.;415 return;416 }417 418 double dphi= DeuxPi/(double)NPhi_[kt];419 phiMin= dphi*(double)(jp);420 phiMax= dphi*(double)(jp+1);421 return;422 }423 424 /* --Methode-- */425 //++426 template <class T>427 int SphereThetaPhi<T>::Index(int kt,int jp) const428 429 // Retourne l'indice du pixel d'indice jp dans la tranche kt430 //--431 {432 int k;433 bool fgzn= false;434 435 // si on se trouve dans l'hemisphere Sud436 if(kt >= NTheta_) {437 fgzn= true;438 kt= 2*NTheta_-1-kt;439 }440 441 // si la tranche kt contient au moins i pixels442 if( (jp>=0) && (jp<NPhi_[kt]) )443 {444 // dans l'hemisphere Sud445 if (fgzn) k= NPix_-TNphi_[kt+1]+jp;446 // dans l'hemisphere Nord447 else k= TNphi_[kt]+jp;448 }449 else450 {451 k= 9999;452 printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);453 }454 return(k);455 }456 457 /* --Methode-- */458 //++459 template <class T>460 void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp)461 462 // Retourne les indices kt et jp du pixel d'indice k463 //--464 {465 bool fgzn= false;466 // si on se trouve dans l'hemisphere Sud467 if(k >= NPix_/2)468 {469 fgzn= true;470 k= NPix_-1-k;471 }472 473 // on recupere l'indice kt de la tranche qui contient le pixel k474 int i;475 for(i = 0; i < NTheta_; i++)476 if(k < TNphi_[i+1]) break;477 478 // indice kt de tranche479 if (fgzn) kt= 2*NTheta_-1-i;480 else kt= i;481 482 // indice jp de pixel483 if (fgzn) jp= TNphi_[i+1]-k-1;484 else jp= k-TNphi_[i];485 }486 //++487 template <class T>488 void SphereThetaPhi<T>::Pixelize(int m)489 490 // effectue le découpage en pixels (m et pet ont la même signification491 // que pour le constructeur)492 //493 // Chaque bande de theta sera découpée en partant de phi=0 ...494 // L'autre hémisphère est parcourue dans le même sens en phi et de495 // l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus496 // proche de l'équateur a z>0 est celui de plus petit phi de la bande la497 // plus proche de l'equateur a z<0).498 //--499 {500 int ntotpix,i,j;501 502 // Decodage et controle des arguments d'appel :503 // au moins 2 et au plus 16384 bandes d'un hemisphere en theta504 if (m < 2) m = 2;505 if (m > 16384) m = 16384;506 507 // On memorise les arguments d'appel508 NTheta_ = m;509 510 // On commence par decouper l'hemisphere z>0.511 // Creation des vecteurs contenant :512 // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)513 Theta_= new double[m+1];514 515 // Le nombre de pixels en phi de chacune des bandes en theta516 NPhi_ = new int[m];517 518 // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une519 // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)520 TNphi_= new int[m+1];521 522 // Calcul du nombre total de pixels dans chaque bande pour optimiser523 // le rapport largeur/hauteur des pixels524 525 //calotte polaire526 TNphi_[0]= 0;527 NPhi_[0] = 1;528 529 //bandes jusqu'a l'equateur530 for(j = 1; j < m; j++)531 {532 TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];533 NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;534 }535 536 // Nombre total de pixels sur l'hemisphere537 ntotpix = TNphi_[m-1]+NPhi_[m-1];538 TNphi_[m]= ntotpix;539 // et sur la sphere entiere540 NPix_= 2*ntotpix;541 542 // Creation et initialisation du vecteur des contenus des pixels543 pixels_.ReSize(NPix_);544 pixels_.Reset();545 546 // Determination des limites des bandes en theta :547 // omeg est l'angle solide couvert par chaque pixel,548 // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg549 // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide550 //de la551 // calotte allant du pole a la limite haute de la bande kt vaut552 // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...553 554 double omeg2pi= 1./(double)ntotpix;555 Omega_ = omeg2pi*DeuxPi;556 557 for(j=0; j <= m; j++)558 {559 Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);560 }561 }562 563 //++564 template <class T>565 void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const566 567 // Retourne, pour la tranche en theta d'indice 'index' le theta568 // correspondant, un vecteur (Peida) contenant les phi des pixels de569 // la tranche, un vecteur (Peida) contenant les valeurs de pixel570 // correspondantes571 //--572 573 {574 cout << "entree GetThetaSlice, couche no " << index << endl;575 576 if(index < 0 || index > NbThetaSlices())577 {578 // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));579 cout << " SphereThetaPhi::GetThetaSlice : exceptions a mettre en place" <<endl;580 THROW(rangeCheckErr);581 }582 583 int iring= Index(index,0);584 int bid = this->NPhi(index);585 int lring = bid;586 cout << " iring= " << iring << " lring= " << lring << endl;587 588 phi.ReSize(lring);589 value.ReSize(lring);590 double Te= 0.;591 double Fi= 0.;592 for(int kk = 0; kk < lring; kk++)593 {594 PixThetaPhi(kk+iring,Te,Fi);595 phi(kk)= Fi;596 value(kk)= PixVal(kk+iring);597 }598 theta= Te;599 }600 601 template <class T>602 void SphereThetaPhi<T>::setmNPhi(int* array, int m)603 //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta604 //--605 {606 NPhi_= new int[m];607 for(int k = 0; k < m; k++) NPhi_[k]= array[k];608 }609 610 template <class T>611 void SphereThetaPhi<T>::setmTNphi(int* array, int m)612 //remplit le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes613 //--614 {615 TNphi_= new int[m];616 for(int k = 0; k < m; k++) TNphi_[k]= array[k];617 }618 619 template <class T>620 void SphereThetaPhi<T>::setmTheta(double* array, int m)621 //remplit le tableau contenant les valeurs limites de theta622 //--623 {624 Theta_= new double[m];625 for(int k = 0; k < m; k++) Theta_[k]= array[k];626 }627 628 template <class T>629 void SphereThetaPhi<T>::setDataBlock(T* data, int m)630 // remplit le vecteur des contenus des pixels631 {632 pixels_.FillFrom(m,data);633 }634 635 template <class T>636 void SphereThetaPhi<T>::print(ostream& os) const637 {638 if(mInfo_) os << " DVList Info= " << *mInfo_ << endl;639 //640 os << " NTheta_= " << NTheta_ << endl;641 os << " NPix_ = " << NPix_ << endl;642 os << " Omega_ = " << Omega_ << endl;643 644 os << " contenu de NPhi_ : ";645 for(int i=0; i < NTheta_; i++)646 {647 if(i%5 == 0) os << endl;648 os << NPhi_[i] <<", ";649 }650 os << endl;651 652 os << " contenu de Theta_ : ";653 for(int i=0; i < NTheta_+1; i++)654 {655 if(i%5 == 0) os << endl;656 os << Theta_[i] <<", ";657 }658 os << endl;659 660 os << " contenu de TNphi_ : ";661 for(int i=0; i < NTheta_+1; i++)662 {663 if(i%5 == 0) os << endl;664 os << TNphi_[i] <<", ";665 }666 os << endl;667 668 os << " contenu de pixels : ";669 for(int i=0; i < NPix_; i++)670 {671 if(i%5 == 0) os << endl;672 os << pixels_(i) <<", ";673 }674 os << endl;675 }676 677 ///////////////////////////////////////////////////////////678 // --------------------------------------------------------679 // Les objets delegues pour la gestion de persistance680 // --------------------------------------------------------681 //////////////////////////////////////////////////////////682 683 template <class T>684 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()685 {686 dobj= new SphereThetaPhi<T>;687 ownobj= true;688 }689 690 template <class T>691 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)692 {693 dobj= new SphereThetaPhi<T>;694 ownobj= true;695 Read(filename);696 }697 698 template <class T>699 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)700 {701 dobj= new SphereThetaPhi<T>(obj);702 ownobj= true;703 }704 705 template <class T>706 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)707 {708 dobj= obj;709 ownobj= false;710 }711 712 template <class T>713 FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()714 {715 if (ownobj && dobj) delete dobj;716 }717 718 template <class T>719 AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()720 {721 return(dobj);722 }723 724 template <class T>725 void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)726 {727 cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;728 729 if(dobj == NULL)730 {731 dobj= new SphereThetaPhi<T>;732 }733 734 // Pour savoir s'il y avait un DVList Info associe735 char strg[256];736 is.GetLine(strg, 255);737 bool hadinfo= false;738 if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo= true;739 if(hadinfo)740 { // Lecture eventuelle du DVList Info741 is >> dobj->Info();742 }743 744 int mNTheta;745 is.GetI4(mNTheta);746 dobj->setSizeIndex(mNTheta);747 748 int mNPix;749 is.GetI4(mNPix);750 dobj->setNbPixels(mNPix);751 752 double mOmeg;753 is.GetR8(mOmeg);754 dobj->setPixSolAngle(mOmeg);755 756 int* mNphi= new int[mNTheta];757 is.GetI4s(mNphi, mNTheta);758 dobj->setmNPhi(mNphi, mNTheta);759 delete [] mNphi;760 761 int* mTNphi= new int[mNTheta+1];762 is.GetI4s(mTNphi, mNTheta+1);763 dobj->setmTNphi(mTNphi, mNTheta+1);764 delete [] mTNphi;765 766 double* mTheta= new double[mNTheta+1];767 is.GetR8s(mTheta, mNTheta+1);768 dobj->setmTheta(mTheta, mNTheta+1);769 delete [] mTheta;770 771 T* mPixels= new T[mNPix];772 PIOSReadArray(is, mPixels, mNPix);773 dobj->setDataBlock(mPixels, mNPix);774 delete [] mPixels;775 }776 777 template <class T>778 void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const779 {780 cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;781 782 if(dobj == NULL)783 {784 cout << " WriteSelf:: dobj= null " << endl;785 return;786 }787 788 char strg[256];789 int mNTheta= dobj->SizeIndex();790 int mNPix = dobj->NbPixels();791 792 if(dobj->ptrInfo())793 {794 sprintf(strg,"SphereThetaPhi: NSlices=%6d NPix=%9d HasInfo",mNTheta,mNPix);795 os.PutLine(strg);796 os << dobj->Info();797 }798 else799 {800 sprintf(strg,"SphereThetaPhi: NSlices=%6d NPix=%9d ",mNTheta,mNPix);801 os.PutLine(strg);802 }803 804 os.PutI4(mNTheta);805 os.PutI4(mNPix);806 os.PutR8(dobj->PixSolAngle(0));807 os.PutI4s(dobj->getmNPhi() , mNTheta);808 os.PutI4s(dobj->getmTNphi(), mNTheta+1);809 os.PutR8s(dobj->getmTheta(), mNTheta+1);810 //os.Put((dobj->getDataBlock())->Data(), mNPix);811 PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);812 }
Note:
See TracChangeset
for help on using the changeset viewer.