source: Sophya/trunk/SophyaExt/FitsIOServer/fitsioserver.cc@ 471

Last change on this file since 471 was 471, checked in by ansari, 26 years ago

creation FitsioServer
File size: 24.2 KB
RevLine 
[471]1
2#include "nbmath.h"
3#include <math.h>
4#include <iostream.h>
5#include "fitsioserver.h"
6#include "dvlist.h"
7#include "strutil.h"
8
9 typedef map<string, MuTyV, less<string> > ValList;
10void FitsIoServer::load(TMatrix<double>& mat, char flnm[])
11{
12 int nbrows=0;
13 int nbcols=0;
14 FITS_tab_typ_ = TDOUBLE;
15 int status=0;
16 long naxis;
17 int n1, n2, n3;
18 DVList dvl;
19 // pointer to the FITS file, defined in fitsio.h
20 fitsfile *fptr;
21 fits_open_file(&fptr, flnm, READONLY, &status);
22 if( status ) printerror( status );
23 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
24 // close the file
25 fits_close_file(fptr, &status);
26 if( status ) printerror( status );
27
28 nbrows=n1;
29 nbcols=n2;
30 if (naxis == 1) nbcols=1;
31 // cout << " lecture du dvlist par load " << endl;
32 dvl.Print();
33 ValList::const_iterator it;
34 for(it = dvl.Begin(); it != dvl.End(); it++)
35 {
36 char datatype= (*it).second.typ;
37 char keyname[]="";
38 strcpy(keyname, (*it).first.substr(0,64).c_str());
39 int ival=0;
40 double dval=0.;
41 char strval[]="";
42 switch (datatype)
43 {
44 case 'I' :
45 ival=(*it).second.mtv.iv;
46 break;
47 case 'D' :
48 dval=(*it).second.mtv.dv;
49 break;
50 case 'S' :
51 strcpy(strval, (*it).second.mtv.strv);
52 break;
53 default :
54 cout << " probleme dans type mot cle optionnel" << endl;
55 break;
56 }
57 }
58
59 if (naxis > 2)
60 {
61 cout << " le fichier fits n'est pas une matrice, naxis= " << naxis << endl;
62 }
63 // number of components
64
65 if (mat.NRows() != nbrows || mat.NCols() != nbcols )
66 {
67 cout << " found " << nbrows << " rows ";
68 cout << " expected " << mat.NRows() << endl;
69 cout << " found " << nbcols << " columns " ;
70 cout << " expected " << mat.NCols() << endl;
71 mat.ReSize(nbrows,nbcols);
72 cout << " resize the vector to nbrows= " << nbrows << " nbcols= " << nbcols << endl;
73 }
74 int ij=0;
75 for (int j=0; j< nbcols; j++)
76 for (int i = 0; i < nbrows; i++) mat(i,j) = dtab_[ij++];
77 // cout << " chargement termine " << endl;
78}
79
80void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
81{
82 int npixels=0;
83 int nside=0;
84 FITS_tab_typ_ = TDOUBLE;
85 read_sphere(flnm, npixels, nside);
86 // number of pixels in the sphere
87 if (sph.NbPixels() != npixels)
88 {
89 cout << " found " << npixels << " pixels" << endl;
90 cout << " expected " << sph.NbPixels() << endl;
91 if (nside <= 0 )
92 {
93 cout << " FITSIOSERVER : no resolution parameter on fits file " << endl;
94 exit(0);
95 }
96 if (nside != sph.SizeIndex())
97 {
98 sph.Resize(nside);
99 cout << " resize the sphere to nside= " << nside << endl;
100 }
101 else
102 {
103 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
104 exit(0);
105 }
106 }
107 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
108 // cout << " chargement termine " << endl;
109}
110void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
111{
112 int npixels=0;
113 int nside=0;
114 FITS_tab_typ_ = TFLOAT;
115 read_sphere(flnm, npixels, nside);
116 // number of pixels in the sphere
117 if (sph.NbPixels() != npixels)
118 {
119 cout << " found " << npixels << " pixels" << endl;
120 cout << " expected " << sph.NbPixels() << endl;
121 if (nside <= 0 )
122 {
123 cout << " FITSIOSERVER : no resolution parameter on fits file " << endl;
124 exit(0);
125 }
126 if (nside != sph.SizeIndex())
127 {
128 sph.Resize(nside);
129 cout << " resize the sphere to nside= " << nside << endl;
130 }
131 else
132 {
133 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
134 exit(0);
135 }
136 }
137 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
138 // cout << " chargement termine " << endl;
139}
140
141void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
142{
143 int nside;
144 int nbrows=0;
145 int nbcols=0;
146 FITS_tab_typ_ = TDOUBLE;
147 int status=0;
148 long naxis;
149 int n1, n2, n3;
150 DVList dvl;
151 // pointer to the FITS file, defined in fitsio.h
152 fitsfile *fptr;
153 fits_open_file(&fptr, flnm, READONLY, &status);
154 if( status ) printerror( status );
155 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
156 // close the file
157 fits_close_file(fptr, &status);
158 if( status ) printerror( status );
159
160 nbrows=n1;
161 nbcols=n2;
162 if (naxis != 2)
163 {
164 cout << " le fichier fits n'est pas une localmap, naxis= " << naxis << endl;
165 }
166 // cout << " lecture du dvlist par load " << endl;
167 dvl.Print();
168 ValList::const_iterator it;
169 for(it = dvl.Begin(); it != dvl.End(); it++)
170 {
171 char datatype= (*it).second.typ;
172 char keyname[]="";
173 strcpy(keyname, (*it).first.substr(0,64).c_str());
174 int ival=0;
175 double dval=0.;
176 char strval[]="";
177 switch (datatype)
178 {
179 case 'I' :
180 ival=(*it).second.mtv.iv;
181 break;
182 case 'D' :
183 dval=(*it).second.mtv.dv;
184 break;
185 case 'S' :
186 strcpy(strval, (*it).second.mtv.strv);
187 break;
188 default :
189 cout << " probleme dans type mot cle optionnel" << endl;
190 break;
191 }
192 if (!strcmp(keyname, "NSIDE") )
193 {
194 // cout << " j'ai trouve NSIDE " << endl;
195 nside = ival;
196 }
197 //if (!strcmp(keyname, "ORDERING") )
198 // {
199 // cout << " j'ai trouve ORDERING " << endl;
200 // }
201 }
202 float theta0 = dvl.GetD("THETA0");
203 float phi0 = dvl.GetD("PHI0");
204 int x0 = dvl.GetI("X0");
205 int y0 = dvl.GetI("Y0");
206 int xo = dvl.GetI("XO");
207 int yo = dvl.GetI("YO");
208 float anglex=dvl.GetD("ANGLEX");
209 float angley=dvl.GetD("ANGLEY");
210 float angle = dvl.GetD("ANGLE");
211 //cout << " theta0= " << theta0 << endl;
212 //cout << " phi0= " << phi0 << endl;
213 //cout << " x0= " << x0 << endl;
214 //cout << " y0= " << y0 << endl;
215 //cout << " angle= " << angle << endl;
216 //cout << " anglex= " << anglex << endl;
217 //cout << " angley= " << angley << endl;
218
219 // number of components
220
221 if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
222 {
223 cout << " found " << nbrows << " x-pixels ";
224 cout << " expected " << lcm.Size_x() << endl;
225 cout << " found " << nbcols << " y-pixels " ;
226 cout << " expected " << lcm.Size_y() << endl;
227 lcm.ReSize(nbrows,nbcols);
228 cout << " resize the map to x-size= " << nbrows << " y-size= " << nbcols << endl;
229 }
230 lcm.SetSize(anglex, angley);
231 lcm.SetOrigin(theta0, phi0, x0, y0, angle);
232 int ij=0;
233 for (int j=0; j< nbcols; j++)
234 for (int i = 0; i < nbrows; i++) lcm(i,j) = dtab_[ij++];
235 //cout << " chargement termine " << endl;
236}
237
238void FitsIoServer::printerror(int status) const
239{
240 //Print out cfitsio error messages and exit program
241 if( status )
242 {
243 //print error report
244 fits_report_error(stderr, status);
245 //terminate the program, returning error status
246 exit( status );
247 }
248 return;
249}
250void FitsIoServer::save( TMatrix<double>& mat, char filename[])
251{
252 int nbrows = mat.NRows();
253 int nbcols = mat.NCols();
254 long naxis = nbcols > 1 ? 2 : 1;
255 //cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
256 FITS_tab_typ_ = TDOUBLE;
257 if (dtab_ != NULL ) delete[] dtab_;
258 dtab_=new double[nbrows*nbcols];
259
260
261 //pointer to the FITS file, defined in fitsio.h
262 fitsfile *fptr;
263
264
265 // initialize status before calling fitsio routines
266 int status = 0;
267
268 // delete old file if it already exists
269 remove(filename);
270
271 // create new FITS file
272 fits_create_file(&fptr, filename, &status);
273 if( status ) printerror( status );
274 int ij=0;
275 for (int j=0; j< nbcols; j++)
276 for (int i = 0; i < nbrows; i++) dtab_[ij++]= mat(i,j);
277
278 DVList dvl;
279
280 // cout << " save : debut ecriture " << endl;
281 planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
282 // cout << " save : fin ecriture " << endl;
283
284
285 // close the file
286 fits_close_file(fptr, &status);
287 if( status ) printerror( status );
288
289}
290
291//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
292void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
293{
294 int npixels = sph.NbPixels();
295 FITS_tab_typ_ = TDOUBLE;
296 if (dtab_ != NULL ) delete[] dtab_;
297 dtab_=new double[npixels];
298
299 //pointer to the FITS file, defined in fitsio.h
300 fitsfile *fptr;
301
302
303 // initialize status before calling fitsio routines
304 int status = 0;
305
306 // delete old file if it already exists
307 remove(filename);
308
309 // create new FITS file
310 fits_create_file(&fptr, filename, &status);
311 if( status ) printerror( status );
312
313 for (int j = 0; j < npixels; j++) dtab_[j]= sph(j);
314 DVList dvl;
315 dvl["NSIDE"] = sph.SizeIndex();
316 dvl["ORDERING"]=sph.TypeOfMap();
317
318 planck_write_img(fptr, 1, npixels, 0, 0, dvl);
319
320 // decider ulterieurement ce qu'on fait de ce qui suit, specifique
321 // pour l'instant, aux spheres gorski. Actuellement, on ne sauve pas
322 // les informations d'analyse harmonique
323 /*
324 // write supplementary keywords
325 fits_write_comment(fptr, " ", &status);
326 if( status ) printerror( status );
327
328 char comment[] = "Generated by synfast";
329 char svalue [] = "SYNFAST ";
330 strcat(svalue, version);
331 fits_write_key(fptr, TSTRING, "ORIG-MAP", svalue, comment, &status);
332 if( status ) printerror( status );
333
334 strcpy(comment,"HEALPIX Pixelisation");
335 strcpy(svalue, "HEALPIX");
336 fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
337 if( status ) printerror( status );
338
339
340 strcpy(comment,"pixel ordering scheme, either RING or NESTED");
341 strcpy(svalue, "RING");
342 fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
343 if( status ) printerror( status );
344
345 strcpy(comment,"Maximum multipole l used in map synthesis");
346 int lmax= sph.nlmax();
347 fits_write_key(fptr, TINT, "MAX-LPOL", &lmax, comment, &status);
348 if( status ) printerror( status );
349
350 strcpy(comment,"Random generator seed");
351 int iseed= sph.iseed();
352 fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
353 if( status ) printerror( status );
354
355 strcpy(comment,"FWHM in DEGREES");
356 float fwhm_deg= sph.fwhm()/60.0;
357 fits_write_key(fptr, TFLOAT, "FWHM", &fwhm_deg, comment, &status);
358 if( status ) printerror( status );
359
360 strcpy(comment,"--------------------------------------------");
361 fits_write_comment(fptr, comment, &status);
362 if( status ) printerror( status );
363
364 strcpy(comment," Above keywords are still likely to change !");
365 fits_write_comment(fptr, comment, &status);
366 if( status ) printerror( status );
367
368 strcpy(comment,"--------------------------------------------");
369 fits_write_comment(fptr, comment, &status);
370 if( status ) printerror( status );
371
372 // write information about the power spectrum
373 char card[]="";
374 sph.powfile(card);
375 strcpy(comment,"Input power spectrum in : ");
376 strcat(comment, card);
377 fits_write_comment(fptr, comment, &status);
378 if( status ) printerror( status );
379
380 strcpy(comment,"Quadrupole :");
381 fits_write_comment(fptr, comment, &status);
382 if( status ) printerror( status );
383
384 strcpy(comment,"C(2)= ");
385 sprintf(card,"%g",sph.quadrupole());
386 strcat(comment, card);
387 fits_write_comment(fptr, comment, &status);
388 if( status ) printerror( status );
389
390 strcpy(comment,"keywords relative to input power spectrum ");
391 fits_write_comment(fptr, comment, &status);
392 if( status ) printerror( status );
393 strcpy(comment,"have not yet been defined for PLANCK Surveyor Project");
394 fits_write_comment(fptr, comment, &status);
395 if( status ) printerror( status );
396 */
397
398 // close the file
399 fits_close_file(fptr, &status);
400 if( status ) printerror( status );
401
402}
403
404void FitsIoServer::save(SphericalMap<float>& sph, char filename[])
405{
406 int npixels = sph.NbPixels();
407 FITS_tab_typ_ = TFLOAT;
408 if (ftab_ != NULL ) delete[] ftab_;
409 ftab_=new float[npixels];
410
411
412 //pointer to the FITS file, defined in fitsio.h
413 fitsfile *fptr;
414
415
416 // initialize status before calling fitsio routines
417 int status = 0;
418
419 // delete old file if it already exists
420 remove(filename);
421
422 // create new FITS file
423 fits_create_file(&fptr, filename, &status);
424 if( status ) printerror( status );
425 for (int j = 0; j < npixels; j++) ftab_[j]= sph(j);
426 DVList dvl;
427 dvl["NSIDE"] = sph.SizeIndex();
428 dvl["ORDERING"]=sph.TypeOfMap();
429
430 planck_write_img(fptr, 1, npixels, 0, 0, dvl);
431
432
433 // close the file
434 fits_close_file(fptr, &status);
435 if( status ) printerror( status );
436
437}
438
439void FitsIoServer::save(LocalMap<double>& locm, char filename[])
440{
441 int nbrows = locm.Size_x();
442 int nbcols = locm.Size_y();
443 long naxis = 2;
444 cout << " nombre de pts en x : " << nbrows << " en y " << nbcols << endl;
445 FITS_tab_typ_ = TDOUBLE;
446 if (dtab_ != NULL ) delete[] dtab_;
447 dtab_=new double[nbrows*nbcols];
448
449 //pointer to the FITS file, defined in fitsio.h
450 fitsfile *fptr;
451
452
453 // initialize status before calling fitsio routines
454 int status = 0;
455
456 // delete old file if it already exists
457 remove(filename);
458
459 // create new FITS file
460 fits_create_file(&fptr, filename, &status);
461 if( status ) printerror( status );
462 int ij=0;
463 for (int j=0; j< nbcols; j++)
464 for (int i = 0; i < nbrows; i++) dtab_[ij++]= locm(i,j);
465
466 DVList dvl;
467 dvl["NSIDE"] = locm.SizeIndex();
468 dvl["ORDERING"]=locm.TypeOfMap();
469 float theta0;
470 float phi0;
471 int x0;
472 int y0;
473 float angle;
474 locm.Origin(theta0, phi0,x0, y0, angle);
475 float anglex;
476 float angley;
477 locm.Aperture(anglex, angley);
478 dvl["THETA0"] = theta0;
479 dvl["PHI0"] = phi0;
480 dvl["X0"] = x0;
481 dvl["Y0"] = y0;
482 dvl["ANGLE"] = angle;
483 dvl["ANGLEX"] = anglex;
484 dvl["ANGLEY"] = angley;
485 planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
486
487
488 // close the file
489 fits_close_file(fptr, &status);
490 if( status ) printerror( status );
491
492}
493
494
495void FitsIoServer::planck_write_img( fitsfile *fptr, int naxis,int n1, int n2, int n3, DVList& dvl)
496{
497
498 int status=0;
499 int bitpix=0;
500 if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
501 if ( FITS_tab_typ_ == TFLOAT) bitpix = FLOAT_IMG;
502 long naxes[3];
503 naxes[0] = n1;
504 naxes[1] = n2;
505 naxes[2] = n3;
506 if (n2 > 0 && naxis < 2) cout << " FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
507 if (n3 > 0 && naxis < 3) cout << " FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
508 fits_create_img(fptr, bitpix, naxis, naxes, &status);
509 if( status ) printerror( status );
510
511
512
513 long nelements= naxes[0];
514 if (naxis >=2) nelements*=naxes[1];
515 if (naxis == 3) nelements*=naxes[2];
516 switch (FITS_tab_typ_)
517 {
518 case TDOUBLE :
519 fits_write_img(fptr, TDOUBLE, 1, nelements, dtab_, &status);
520 if( status ) printerror( status );
521 break;
522 case TFLOAT :
523 fits_write_img(fptr, TFLOAT, 1, nelements, ftab_, &status);
524 if( status ) printerror( status );
525 break;
526 default :
527 cout << " fitsioserver : type de tableau non traite en ecriture " << endl;
528 break;
529 }
530
531 // write the current date
532 fits_write_date(fptr, &status);
533 if( status ) printerror( status );
534 // on convient d ecrire apres la date, les mots cles utilisateur.
535 // la date servira de repere pour relire ces mots cles.
536
537 dvl.Print();
538 ValList::const_iterator it;
539 for(it = dvl.Begin(); it != dvl.End(); it++)
540 {
541 int datatype= key_type_PL2FITS( (*it).second.typ);
542 char keyname[]="";
543 int flen_keyword = 9;
544 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
545 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
546 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
547 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
548 char comment[FLEN_COMMENT]="";
549 int ival=0;
550 double dval=0.;
551 char strval[]="";
552 switch (datatype)
553 {
554 case TINT :
555 ival=(*it).second.mtv.iv;
556 strcpy(comment,"I entier");
557 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
558 break;
559 case TDOUBLE :
560 dval=(*it).second.mtv.dv;
561 strcpy(comment,"D double");
562 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
563 break;
564 case TSTRING :
565 strcpy(strval, (*it).second.mtv.strv);
566 strcpy(comment,"S character string");
567 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
568 break;
569 default :
570 cout << " probleme dans type mot cle optionnel" << endl;
571 break;
572 }
573 if( status ) printerror( status );
574 }
575}
576void FitsIoServer::planck_read_img( fitsfile *fptr, long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
577{
578 int status=0;
579 long bitpix;
580 long naxes[3]={0,0,0};
581 char* comment=NULL;
582
583 fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
584 if( status ) printerror( status );
585 fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
586 if( status ) printerror( status );
587 int nfound;
588 int nkeys=(int)naxis;
589 fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
590 if( status ) printerror( status );
591
592 n1 = naxes[0] ;
593 n2 = naxes[1] ;
594 n3 = naxes[2] ;
595
596
597 long nelements= naxes[0];
598 if (naxis >=2) nelements*=naxes[1];
599 if (naxis == 3) nelements*=naxes[2];
600 int anynull;
601 double dnullval=0.;
602 float fnullval=0.;
603 // on laisse a fits le soin de convertir le type du tableau lu vers
604 // le type suppose par l'utilisateur de fitsioserver
605 //
606 switch ( FITS_tab_typ_)
607 {
608 case TDOUBLE :
609 if (dtab_ != NULL) delete [] dtab_;
610 dtab_=new double[nelements];
611 fits_read_img(fptr, TDOUBLE, 1, nelements, &dnullval, dtab_,
612 &anynull, &status);
613 if( status ) printerror( status );
614 break;
615 case TFLOAT :
616 if (ftab_ != NULL) delete [] ftab_;
617 ftab_=new float[nelements];
618 fits_read_img(fptr, TFLOAT, 1, nelements, &fnullval, ftab_,
619 &anynull, &status);
620 if( status ) printerror( status );
621 break;
622 default :
623 cout << " fitsio::read_img : type non traite: " << FITS_tab_typ_ << endl;
624 break;
625 }
626 status = 0;
627 char card[FLEN_CARD];
628 int num = 0;
629 char comment2[FLEN_COMMENT] = "x";
630 char keyname[]= "";
631 char datekey[]= "DATE";
632 char endkey[] = "END";
633 char typ='x';
634 int ival;
635 double dval;
636 char strval[]="";
637 // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
638 // on va jusqu'au mot cle DATE
639 int flen_keyword = 9;
640 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
641 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
642 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
643 while (status == 0 && strncmp(keyname, datekey,4) != 0 )
644 {
645 num++;
646 fits_read_record(fptr, num, card, &status);
647 strncpy(keyname,card,flen_keyword-1);
648 }
649 if (status != 0 )
650 {
651 cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
652 }
653 // on recupere la liste des mots-cles utilisateurs
654 while (status == 0)
655 {
656 num++;
657 // on lit un record pour recuperer le nom du mot-cle
658 fits_read_record(fptr, num, card, &status);
659 strncpy(keyname,card,flen_keyword-1);
660 char value[FLEN_VALUE];
661 // on recupere le premier caractere du commentaire, qui contient
662 // le code du type de la valeur
663 // (tant que ce n est pas le mot cle END)
664 fits_read_keyword(fptr, keyname, value, comment2, &status);
665 if ( strncmp(keyname, endkey,flen_keyword-1) != 0)
666 {
667 typ = comment2[0];
668 // quand le type est connu, on lit la valeur
669 switch (typ)
670 {
671 case 'I' :
672 fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
673 if( status ) printerror( status );
674 strip (keyname, 'B',' ');
675 dvl[keyname] = ival;
676 break;
677 case 'D' :
678 fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
679 if( status ) printerror( status );
680 strip (keyname, 'B',' ');
681 dvl[keyname] = dval;
682 break;
683 case 'S' :
684 fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
685 if( status ) printerror( status );
686 strip (keyname, 'B',' ');
687 strip(strval, 'B',' ');
688 dvl[keyname]=strval;
689 break;
690 default :
691 cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
692 break;
693 }
694 }
695 }
696
697}
698
699
700void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
701{
702
703 long naxes[2]={600, 300};
704 float map [ 600*300 ];
705 int npixels= naxes[0]*naxes[1];
706
707 cout << " image FITS en projection SINUS" << endl;
708 // table will have npixels rows
709 for(int j=0; j < npixels; j++) map[j]=0.;
710 for(int j=0; j<naxes[1]; j++)
711 {
712 double yd = (j+0.5)/naxes[1]-0.5;
713 double theta = (0.5-yd)*Pi;
714 double facteur=1./sin(theta);
715 for(int i=0; i<naxes[0]; i++)
716 {
717 double xa = (i+0.5)/naxes[0]-0.5;
718 double phi = 2.*Pi*xa*facteur+Pi;
719 float th=float(theta);
720 float ff=float(phi);
721 if (phi<2*Pi && phi>= 0)
722 {
723 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
724 }
725 }
726 }
727
728 write_picture(naxes, map, filename);
729
730}
731void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
732{
733 // le code de cete methode duplique celui de la precedente, la seule
734 //difference etant le type de sphere en entree. Ces methodes de projection
735 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
736 // ne me suis pas casse la tete, pour l instant
737
738 long naxes[2]={600, 300};
739 float map [ 600*300 ];
740 int npixels= naxes[0]*naxes[1];
741
742 cout << " image FITS en projection SINUS" << endl;
743 // table will have npixels rows
744 for(int j=0; j < npixels; j++) map[j]=0.;
745 for(int j=0; j<naxes[1]; j++)
746 {
747 double yd = (j+0.5)/naxes[1]-0.5;
748 double theta = (0.5-yd)*Pi;
749 double facteur=1./sin(theta);
750 for(int i=0; i<naxes[0]; i++)
751 {
752 double xa = (i+0.5)/naxes[0]-0.5;
753 double phi = 2.*Pi*xa*facteur+Pi;
754 float th=float(theta);
755 float ff=float(phi);
756 if (phi<2*Pi && phi>= 0)
757 {
758 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
759 }
760 }
761 }
762
763 write_picture(naxes, map, filename);
764
765
766}
767
768void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
769{
770
771 long naxes[2];
772 naxes[0] = lcm.Size_x();
773 naxes[1] = lcm.Size_y();
774 int npixels= naxes[0]*naxes[1];
775 float* map = new float[npixels];
776
777 // table will have npixels rows
778 for(int j=0; j < npixels; j++) map[j]=0.;
779 for(int j=0; j<naxes[1]; j++)
780 {
781 for(int i=0; i<naxes[0]; i++)
782 {
783 map[j*naxes[0]+i] = lcm(i, j);
784 }
785 }
786
787 write_picture(naxes, map, filename);
788 delete [] map;
789}
790
791void FitsIoServer::read_sphere(char flnm[], int& npixels, int& nside)
792{
793 int status=0;
794 long naxis;
795 nside=0;
796 int n1, n2, n3;
797 DVList dvl;
798 // pointer to the FITS file, defined in fitsio.h
799 fitsfile *fptr;
800 fits_open_file(&fptr, flnm, READONLY, &status);
801 if( status ) printerror( status );
802 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
803 // close the file
804 fits_close_file(fptr, &status);
805 if( status ) printerror( status );
806
807 npixels=n1;
808 dvl.Print();
809 ValList::const_iterator it;
810 for(it = dvl.Begin(); it != dvl.End(); it++)
811 {
812 char datatype= (*it).second.typ;
813 char keyname[]="";
814 strcpy(keyname, (*it).first.substr(0,64).c_str());
815 int ival=0;
816 double dval=0.;
817 char strval[]="";
818 switch (datatype)
819 {
820 case 'I' :
821 ival=(*it).second.mtv.iv;
822 break;
823 case 'D' :
824 dval=(*it).second.mtv.dv;
825 break;
826 case 'S' :
827 strcpy(strval, (*it).second.mtv.strv);
828 break;
829 default :
830 cout << " probleme dans type mot cle optionnel" << endl;
831 break;
832 }
833 if (!strcmp(keyname, "NSIDE") )
834 {
835 // cout << " j'ai trouve NSIDE " << endl;
836 nside = ival;
837 }
838 // if (!strcmp(keyname, "ORDERING") )
839 /// {
840 // cout << " j'ai trouve ORDERING " << endl;
841 // }
842 }
843
844 if (naxis != 1)
845 {
846 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
847 }
848}
849
850
851void FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
852{
853 //pointer to the FITS file, defined in fitsio.h
854 fitsfile *fptr;
855 int bitpix = FLOAT_IMG;
856 long naxis = 2;
857
858 // delete old file if it already exists
859 remove(filename);
860
861 // initialize status before calling fitsio routines
862 int status = 0;
863
864 // create new FITS file
865 fits_create_file(&fptr, filename, &status);
866 if( status ) printerror( status );
867
868 // write the required header keywords
869 fits_create_img(fptr, bitpix, naxis, naxes, &status);
870 if( status ) printerror( status );
871
872 // write the current date
873 fits_write_date(fptr, &status);
874 if( status ) printerror( status );
875
876
877 // first row in table to write
878 long firstrow = 1;
879 // first element in row
880 long firstelem = 1;
881 int colnum = 1;
882 int nelements=naxes[0]*naxes[1];
883 fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
884 if( status ) printerror( status );
885
886 // close the file
887 fits_close_file(fptr, &status);
888 if( status ) printerror( status );
889}
890
891
Note: See TracBrowser for help on using the repository browser.