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

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

ajout save sphereGorski
File size: 52.4 KB
Line 
1//************************************************************************
2// Class for loadind and saving from FITS-formatted file to DPC objects
3// (G. Le Meur ; Francois Touze) OCT. 99
4//
5// methods 'load(X& x, char f[])' get from FITS file "f" a DPC object x
6// from DPC class X.
7// methods 'save(X& x, char f[])' save a DPC object x from DPC // class X
8// onto a FITS file "f" .
9
10//************************************************************************
11
12#include <iostream.h>
13#include <list>
14
15#include "fitsioserver.h"
16#include "strutil.h"
17
18void FitsIoServer::load(TMatrix<double>& mat,char flnm[])
19{
20 int nbrows=0;
21 int nbcols=0;
22 FITS_tab_typ_ = TDOUBLE;
23 long naxis;
24 int n1, n2, n3;
25 DVList dvl;
26 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
27
28 nbrows=n1;
29 nbcols=n2;
30 if (naxis == 1) nbcols=1;
31 if (naxis > 2)
32 {
33 cout<<" FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
34 }
35
36 // number of components
37 if (mat.NRows() != nbrows || mat.NCols() != nbcols )
38 {
39 cout << " found " << nbrows << " rows ";
40 cout << " expected " << mat.NRows() << endl;
41 cout << " found " << nbcols << " columns " ;
42 cout << " expected " << mat.NCols() << endl;
43 mat.ReSize(nbrows,nbcols);
44 cout << " resize the vector to nbrows= " << nbrows << " nbcols= " << nbcols << endl;
45 }
46 int ij=0;
47 for (int j=0; j< nbcols; j++)
48 for (int i = 0; i < nbrows; i++) mat(i,j) = (double)r_8tab_[ij++];
49}
50
51void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum)
52
53 //********************************************************/
54 //* move to the HDU which has the specified number hdunum*/
55 //* in the FITS file, perform read operations and write */
56 //* the elements in an NTuple. */
57 //********************************************************/
58{
59
60 // pointer to the FITS file, defined in fitsio.h
61 fitsfile *fptr;
62 int status = 0;
63 if( fits_open_file(&fptr,flnm,READONLY,&status) )
64 printerror( status );
65
66 // move to the HDU
67 int hdutype;
68 if( fits_movabs_hdu(fptr,hdunum,&hdutype,&status) )
69 printerror( status );
70
71 // get number of keywords
72 int nkeys,keypos;
73 if( fits_get_hdrpos(fptr,&nkeys,&keypos,&status) )
74 printerror( status );
75
76
77 if (hdutype == BINARY_TBL)
78 printf("\nReading binary table in HDU %d:\n",hdunum);
79 else
80 {
81 printf("Error:: this HDU is not a binary table\n");
82 exit ( status );
83 }
84
85 // number of columns
86 int tfields;
87 if( fits_get_num_cols(fptr,&tfields,&status) )
88 printerror( status );
89
90 // to get table size
91 long naxis[2];
92 int nfound;
93 if( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxis, &nfound, &status) )
94 printerror( status );
95 int nrows= naxis[1];
96
97 //Information about each column
98 char **ttype, **tform;
99 ttype= new char*[tfields];
100 tform= new char*[tfields];
101 for( int ii = 0; ii < tfields; ii++)
102 {
103 ttype[ii]= new char[FLEN_VALUE];
104 tform[ii]= new char[FLEN_VALUE];
105 }
106
107 // number of TTYPEn,TFORMn and TUNITn keywords found
108 int num= 0;
109
110 // read the column names from the TTYPEn keywords
111 fits_read_keys_str(fptr,"TTYPE",1,tfields,ttype,&nfound,&status);
112 num += nfound;
113 // read the column names from the TFORMn keywords
114 fits_read_keys_str(fptr,"TFORM",1,tfields,tform,&nfound,&status);
115 num += nfound;
116
117 //for(int ii = 0; ii < tfields; ii++)
118 // printf("\nColumn name & Format %8s %8s", ttype[ii],tform[ii]);
119
120 // select the columns with float data values
121 int typecode;
122 long repeat,width;
123 list<int> column;
124 for( int ii = 0; ii < tfields; ii++)
125 {
126 fits_binary_tform(tform[ii],&typecode, &repeat, &width, &status);
127 //printf("\n%4s %3d %2ld %2ld", tform[ii], typecode, repeat, width);
128 if(typecode == TFLOAT)
129 column.push_back(ii+1);
130 }
131
132 // get input elements to create the NTuple
133 char **clname;
134 clname= new char*[column.size()];
135 for(int ii = 0; ii < column.size(); ii++)
136 clname[ii]= new char[FLEN_VALUE];
137
138 list<int>::iterator itr;
139 int index= 0;
140 for(itr= column.begin(); itr != column.end(); itr++)
141 strcpy(clname[index++],ttype[*itr-1]);
142
143 for( int ii = 0; ii < tfields; ii++)
144 {
145 delete [] ttype[ii];
146 delete [] tform[ii];
147 }
148 delete [] ttype;
149 delete [] tform;
150
151
152 // check if the specified keyword BLK exists
153 int blk= 512;
154 if(check_keyword(fptr,nkeys,"BLK"))
155 {
156 if( fits_read_key(fptr,TINT,"BLK",&blk,NULL,&status) )
157 printerror( status );
158 }
159 // create a NTuple
160 NTuple nt0(column.size(),clname,blk);
161
162 for(int ii = 0; ii < column.size(); ii++)
163 delete [] clname[ii];
164 delete [] clname;
165
166 float value[1];
167 long felem = 1;
168 char strnull[10];
169 strcpy(strnull, " ");
170 int anynull= 0;
171 long longnull = 0;
172
173 // value to represent undefined array elements
174 float floatnull = FLOATNULLVALUE;
175
176 // read elements from columns and fill NTuple
177 for (int k = 0; k < nrows; k++)
178 {
179 int j= 0;
180 float *xnt= new float[column.size()];
181 list<int>::iterator itr;
182 for(itr= column.begin(); itr != column.end(); itr++)
183 {
184 fits_read_col(fptr,TFLOAT,*itr,k+1,felem,1,&floatnull,value,
185 &anynull, &status);
186 xnt[j++]= value[0];
187 }
188 nt0.Fill(xnt);
189 delete[] xnt;
190 }
191
192 // the TUNITn keywords are optional, if they exist they are put
193 // in the DVLIst object
194 char keyname[LEN_KEYWORD]= "";
195 char strval[FLEN_VALUE];
196 for(int ii = 0; ii < tfields; ii++)
197 {
198 fits_make_keyn("TUNIT",ii+1,keyname,&status);
199 if(check_keyword(fptr,nkeys,keyname))
200 {
201 num++;
202 if( fits_read_key_str(fptr,keyname,strval,NULL,&status) )
203 printerror(status);
204 strip (keyname, 'B',' ');
205 strip(strval, 'B',' ');
206 nt0.Info()[keyname]= strval;
207 }
208 }
209
210 // add the number of mandatory keywords of a binary table
211 num += 8;
212
213 // put names and values of other reserved keywords in a DVList object
214 char comment[FLEN_COMMENT];
215 char dtype;
216 for(int j = num+1; j <= nkeys; j++)
217 {
218 fits_read_keyn(fptr,j,keyname,strval,comment,&status);
219 fits_get_keytype(strval,&dtype,&status);
220 strip (keyname, 'B',' ');
221 strip(strval, 'B',' ');
222
223 switch( dtype )
224 {
225 case 'C':
226 nt0.Info()[keyname]= strval;
227 break;
228 case 'I':
229 int ival;
230 ctoi(strval,&ival);
231 nt0.Info()[keyname]= (int_4)ival;
232 break;
233 case 'F':
234 double dval;
235 ctof(strval,&dval);
236 nt0.Info()[keyname]= dval;
237 break;
238 }
239 }
240
241 // copy in the input NTuple
242 ntpl= nt0;
243
244 if( fits_close_file(fptr, &status) )
245 printerror(status);
246
247 printf("\n");
248 return;
249}
250
251
252void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
253{
254 int npixels=0;
255 int nside=0;
256 long naxis;
257 int n1, n2, n3;
258
259 FITS_tab_typ_ = TDOUBLE;
260
261 DVList dvl;
262 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
263 if (naxis != 1)
264 {
265 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
266 }
267 npixels=n1;
268 nside= dvl.GetI("NSIDE");
269
270 // number of pixels in the sphere
271 if (sph.NbPixels() != npixels)
272 {
273 cout << " found " << npixels << " pixels" << endl;
274 cout << " expected " << sph.NbPixels() << endl;
275 if (nside <= 0 )
276 {
277 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
278 exit(0);
279 }
280 if (nside != sph.SizeIndex())
281 {
282 sph.Resize(nside);
283 cout << " resizing the sphere to nside= " << nside << endl;
284 }
285 else
286 {
287 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
288 exit(0);
289 }
290 }
291 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
292}
293
294
295void FitsIoServer::load(SphereGorski<float>& sph, char flnm[], int hdunum)
296{
297 int npixels=0;
298 int nside=0;
299
300 FITS_tab_typ_ = TFLOAT;
301
302 DVList dvl;
303 planck_read_bntbl(flnm, hdunum, npixels, dvl);
304 //dvl.Print();
305 nside= dvl.GetI("NSIDE");
306 const char* ordering= dvl.GetS("ORDERING").c_str();
307 char* ring = "'RING";
308 if ( strncmp(ordering, ring,3) != 0)
309
310 // if (ordering!="RING ")
311 {
312 cout << " numerotation non RING" << endl;
313 }
314 // number of pixels in the sphere
315 if (sph.NbPixels() != npixels)
316 {
317 cout << " found " << npixels << " pixels" << endl;
318 cout << " expected " << sph.NbPixels() << endl;
319 if (nside <= 0 )
320 {
321 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
322 exit(0);
323 }
324 if (nside != sph.SizeIndex())
325 {
326 sph.Resize(nside);
327 cout << " resizing the sphere to nside= " << nside << endl;
328 }
329 else
330 {
331 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
332 exit(0);
333 }
334 }
335 // cout << " fin relecture debut transfert" << endl;
336 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
337}
338void FitsIoServer::load(SphereGorski<double>& sph, char flnm[], int hdunum)
339{
340 int npixels=0;
341 int nside=0;
342
343 FITS_tab_typ_ = TDOUBLE;
344
345 DVList dvl;
346 planck_read_bntbl(flnm, hdunum, npixels, dvl);
347 //dvl.Print();
348 nside= dvl.GetI("NSIDE");
349 const char* ordering= dvl.GetS("ORDERING").c_str();
350 char* ring = "'RING";
351 if ( strncmp(ordering, ring,3) != 0)
352
353 // if (ordering!="RING ")
354 {
355 cout << " numerotation non RING" << endl;
356 }
357 // number of pixels in the sphere
358 if (sph.NbPixels() != npixels)
359 {
360 cout << " found " << npixels << " pixels" << endl;
361 cout << " expected " << sph.NbPixels() << endl;
362 if (nside <= 0 )
363 {
364 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
365 exit(0);
366 }
367 if (nside != sph.SizeIndex())
368 {
369 sph.Resize(nside);
370 cout << " resizing the sphere to nside= " << nside << endl;
371 }
372 else
373 {
374 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
375 exit(0);
376 }
377 }
378 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
379}
380
381void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
382{
383 int npixels=0;
384 int nside=0;
385 long naxis;
386 int n1, n2, n3;
387 DVList dvl;
388
389 FITS_tab_typ_ = TFLOAT;
390
391 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
392 if (naxis != 1)
393 {
394 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
395 }
396 npixels=n1;
397 nside= dvl.GetI("NSIDE");
398
399 // number of pixels in the sphere
400 if (sph.NbPixels() != npixels)
401 {
402 cout << " found " << npixels << " pixels" << endl;
403 cout << " expected " << sph.NbPixels() << endl;
404 if (nside <= 0 )
405 {
406 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
407 exit(0);
408 }
409 if (nside != sph.SizeIndex())
410 {
411 sph.Resize(nside);
412 cout << " resizing the sphere to nside= " << nside << endl;
413 }
414 else
415 {
416 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
417 exit(0);
418 }
419 }
420 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
421}
422
423void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
424{
425 int nbrows=0;
426 int nbcols=0;
427 FITS_tab_typ_ = TDOUBLE;
428 long naxis;
429 int n1, n2, n3;
430 DVList dvl;
431 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
432
433 nbrows=n1;
434 nbcols=n2;
435 if (naxis != 2)
436 {
437 cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
438 }
439 float theta0 = dvl.GetD("THETA0");
440 float phi0 = dvl.GetD("PHI0");
441 int x0 = dvl.GetI("X0");
442 int y0 = dvl.GetI("Y0");
443 int xo = dvl.GetI("XO");
444 int yo = dvl.GetI("YO");
445 float anglex=dvl.GetD("ANGLEX");
446 float angley=dvl.GetD("ANGLEY");
447 float angle = dvl.GetD("ANGLE");
448
449 // number of components
450
451 if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
452 {
453 cout << " found " << nbrows << " x-pixels ";
454 cout << " expected " << lcm.Size_x() << endl;
455 cout << " found " << nbcols << " y-pixels " ;
456 cout << " expected " << lcm.Size_y() << endl;
457 lcm.ReSize(nbrows,nbcols);
458 cout << " resizing the map to x-size= " << nbrows << " y-size= " << nbcols << endl;
459 }
460 lcm.SetSize(anglex, angley);
461 lcm.SetOrigin(theta0, phi0, x0, y0, angle);
462 int ij=0;
463 for (int j=0; j< nbcols; j++)
464 for (int i = 0; i < nbrows; i++) lcm(i,j) = (double)r_8tab_[ij++];
465}
466
467void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
468{
469 FITS_tab_typ_ = TFLOAT;
470 long naxis;
471 int siz_x;
472 int siz_y;
473 int n3;
474 DVList dvl;
475
476 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
477
478
479 if (naxis != 2)
480 {
481 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
482 }
483
484 DpcImg.Allocate(siz_x, siz_y, 0., 0);
485 float* pixelPtr=DpcImg.ImagePtr();
486 // verifications de type
487 PBaseDataTypes dataT=DataType((r_4)0);
488 int TypeDonnees=dvl.GetI("DATATYPE");
489 if (int(dataT) != TypeDonnees)
490 {
491 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not float " << endl;
492 cout << " eventual conversion to float is achieved by cfitsio lib " << endl;
493 }
494
495 memcpy(pixelPtr, r_4tab_, siz_x*siz_y*DataSize(dataT));
496 const char* nom=dvl.GetS("NAME").c_str();
497 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
498 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
499 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
500 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
501 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
502 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
503
504}
505
506
507void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
508{
509 FITS_tab_typ_ = TINT;
510 long naxis;
511 int siz_x;
512 int siz_y;
513 int n3;
514 DVList dvl;
515
516 // pointer to the FITS file, defined in fitsio.h
517 //fitsfile *fptr;
518 // initialize status before calling fitsio routines
519 //int status = 0;
520
521 //fits_open_file(&fptr, flnm, READONLY, &status);
522 //if( status ) printerror( status );
523 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
524 // close the file
525 //fits_close_file(fptr, &status);
526 //if( status ) printerror( status );
527
528 if (naxis != 2)
529 {
530 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
531 }
532
533 DpcImg.Allocate(siz_x, siz_y, 0., 0);
534 int_4* pixelPtr=DpcImg.ImagePtr();
535
536
537 // verifications de type
538 PBaseDataTypes dataT=DataType((int_4)0);
539 int TypeDonnees=dvl.GetI("DATATYPE");
540 if (int(dataT) != TypeDonnees)
541 {
542 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not int_4 " << endl;
543 cout << " eventual conversion to int_4 is achieved by cfitsio lib " << endl;
544 }
545 memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
546 const char* nom=dvl.GetS("NAME").c_str();
547 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
548 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
549 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
550 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
551 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
552 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
553}
554
555
556void FitsIoServer::save( TMatrix<double>& mat, char filename[])
557{
558 int nbrows = mat.NRows();
559 int nbcols = mat.NCols();
560 long naxis = nbcols > 1 ? 2 : 1;
561 cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
562 FITS_tab_typ_ = TDOUBLE;
563 if (r_8tab_ != NULL ) delete[] r_8tab_;
564 r_8tab_=new r_8[nbrows*nbcols];
565
566
567 int ij=0;
568 for (int j=0; j< nbcols; j++)
569 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)mat(i,j);
570
571 DVList dvl;
572
573 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
574 delete[] r_8tab_;
575}
576
577void FitsIoServer::save(NTuple& ntpl,char flnm[])
578
579 //****************************************************/
580 //* read the elements of the NTuple ntpl, and create */
581 //* a FITS file with a binary table extension */
582 //****************************************************/
583{
584
585 // delete old file if it already exists
586 remove(flnm);
587
588 // pointer to the FITS file, defined in fitsio.h
589 fitsfile *fptr;
590 int status = 0;
591 if( fits_create_file(&fptr, flnm, &status) )
592 printerror( status );
593
594 // table will have tfields columns
595 int tfields= ntpl.NVar();
596
597 // table will have nrows rows
598 int nrows= ntpl.NEntry();
599
600 // extension name
601 char extname[] = "NTuple_Binary_tbl";
602
603 // define the name, and the datatype for the tfields columns
604 char **ttype, **tform;
605 ttype= new char*[tfields];
606 tform= new char*[tfields];
607 for(int i = 0; i < tfields; i++)
608 {
609 ttype[i]= new char[FLEN_VALUE];
610 strcpy(ttype[i], ntpl.NomIndex(i));
611 tform[i]= new char[FLEN_VALUE];
612 strcpy(tform[i], "1E");
613 }
614
615 // create a new empty binary table onto the FITS file
616 // physical units if they exist, are defined in the DVList object
617 // so the null pointer is given for the tunit parameters.
618 if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
619 NULL,extname,&status) )
620 printerror( status );
621
622 for( int ii = 0; ii < tfields; ii++)
623 {
624 delete [] ttype[ii];
625 delete [] tform[ii];
626 }
627 delete [] ttype;
628 delete [] tform;
629
630
631 // first row in table to write
632 int firstrow = 1;
633
634 // first element in row (ignored in ASCII tables)
635 int firstelem = 1;
636
637 for(int i = 0; i < tfields; i++)
638 {
639 float *dens= new float[nrows];
640 for(int j = 0; j < nrows; j++)
641 {
642 dens[j]= ntpl.GetVal(j,i);
643 }
644 fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
645 delete []dens;
646 }
647
648 // number of blocks per event
649 int blk= ntpl.BLock();
650 fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
651
652 // get names and values from the join DVList object
653 DVList dvl= ntpl.Info();
654 dvl.Print();
655 DVList::ValList::const_iterator it;
656 for(it = dvl.Begin(); it != dvl.End(); it++)
657 {
658 char keytype= (*it).second.typ;
659 char keyname[9]="";
660 strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
661 char comment[FLEN_COMMENT]="";
662
663 switch (keytype)
664 {
665 case 'I' :
666 {
667 int ival=(*it).second.mtv.iv;
668 strcpy(comment,"I entier");
669 fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
670 break;
671 }
672 case 'D' :
673 {
674 double dval=(*it).second.mtv.dv;
675 strcpy(comment,"D double");
676 fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
677 break;
678 }
679 case 'S' :
680 {
681 char strval[]="";
682 strcpy(strval,(*it).second.mtv.strv);
683 strcpy(comment,"S character string");
684 fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
685 break;
686 }
687 }
688 }
689
690 //close the FITS file
691 if ( fits_close_file(fptr, &status) )
692 printerror( status );
693 return;
694}
695
696
697
698
699//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
700void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
701{
702 int npixels = sph.NbPixels();
703 FITS_tab_typ_ = TDOUBLE;
704 if (r_8tab_ != NULL ) delete[] r_8tab_;
705 r_8tab_=new r_8[npixels];
706
707
708 for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
709 DVList dvl;
710 dvl["NSIDE"] = (int_4) sph.SizeIndex();
711 dvl["ORDERING"]=sph.TypeOfMap();
712
713 planck_write_img(filename, 1, npixels, 0, 0, dvl);
714
715 // decider ulterieurement ce qu'on fait de ce qui suit, specifique
716 // pour l'instant, aux spheres gorski.
717
718 /*
719 // write supplementary keywords
720 fits_write_comment(fptr, " ", &status);
721 if( status ) printerror( status );
722
723
724 strcpy(comment,"HEALPIX Pixelisation");
725 strcpy(svalue, "HEALPIX");
726 fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
727 if( status ) printerror( status );
728
729
730 strcpy(comment,"pixel ordering scheme, either RING or NESTED");
731 strcpy(svalue, "RING");
732 fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
733 if( status ) printerror( status );
734
735
736 strcpy(comment,"Random generator seed");
737 int iseed= sph.iseed();
738 fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
739 if( status ) printerror( status );
740
741
742 strcpy(comment,"--------------------------------------------");
743 fits_write_comment(fptr, comment, &status);
744 if( status ) printerror( status );
745
746 strcpy(comment," Above keywords are still likely to change !");
747 fits_write_comment(fptr, comment, &status);
748 if( status ) printerror( status );
749
750 strcpy(comment,"--------------------------------------------");
751 fits_write_comment(fptr, comment, &status);
752 if( status ) printerror( status );
753
754 */
755
756 delete[] r_8tab_;
757}
758
759void FitsIoServer::save(SphericalMap<float>& sph, char filename[])
760{
761 int npixels = sph.NbPixels();
762 FITS_tab_typ_ = TFLOAT;
763 if (r_4tab_ != NULL ) delete[] r_4tab_;
764 r_4tab_=new r_4[npixels];
765
766
767 for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
768 DVList dvl;
769 dvl["NSIDE"] = (int_4)sph.SizeIndex();
770 dvl["ORDERING"]=sph.TypeOfMap();
771
772 planck_write_img(filename, 1, npixels, 0, 0, dvl);
773 delete[] r_4tab_;
774
775}
776void FitsIoServer::save(SphereGorski<float>& sph, char filename[])
777{
778 int npixels = sph.NbPixels();
779 FITS_tab_typ_ = TFLOAT;
780 if (r_4tab_ != NULL ) delete[] r_4tab_;
781 r_4tab_=new r_4[npixels];
782
783
784 for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
785 DVList dvl;
786 dvl.SetS("PIXTYPE","HEALPIX ");
787 dvl["ORDERING"]=sph.TypeOfMap();
788 dvl["NSIDE"] = (int_4)sph.SizeIndex();
789 dvl["FIRSTPIX"]=0;
790 dvl["LASTPIX"]=npixels-1;
791 char* typeOfContent="TEMPERATURE";
792 char* extname="SIMULATION";
793 char* comment1=" Sky Map Pixelisation Specific Keywords";
794 planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
795 delete[] r_4tab_;
796
797}
798void FitsIoServer::save(SphereGorski<double>& sph, char filename[])
799{
800 int npixels = sph.NbPixels();
801 FITS_tab_typ_ = TDOUBLE;
802 if (r_8tab_ != NULL ) delete[] r_8tab_;
803 r_8tab_=new r_8[npixels];
804
805
806 for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
807 DVList dvl;
808 dvl.SetS("PIXTYPE","HEALPIX ");
809 dvl["ORDERING"]=sph.TypeOfMap();
810 dvl["NSIDE"] = (int_4)sph.SizeIndex();
811 dvl["FIRSTPIX"]=0;
812 dvl["LASTPIX"]=npixels-1;
813 char* typeOfContent="TEMPERATURE";
814 char* extname="SIMULATION";
815 char* comment1=" Sky Map Pixelisation Specific Keywords";
816 planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
817 delete[] r_8tab_;
818
819}
820
821void FitsIoServer::save(LocalMap<double>& locm, char filename[])
822{
823 int nbrows = locm.Size_x();
824 int nbcols = locm.Size_y();
825 long naxis = 2;
826 cout << " nombre de pts en x : " << nbrows << " en y " << nbcols << endl;
827 FITS_tab_typ_ = TDOUBLE;
828 if (r_8tab_ != NULL ) delete[] r_8tab_;
829 r_8tab_=new r_8[nbrows*nbcols];
830
831 int ij=0;
832 for (int j=0; j< nbcols; j++)
833 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)locm(i,j);
834
835 DVList dvl;
836 dvl["NSIDE"] = (int_4) locm.SizeIndex();
837 dvl["ORDERING"]=locm.TypeOfMap();
838 double theta0;
839 double phi0;
840 int x0;
841 int y0;
842 double angle;
843 locm.Origin(theta0,phi0,x0,y0,angle);
844 double anglex;
845 double angley;
846 locm.Aperture(anglex,angley);
847 dvl["THETA0"] = theta0;
848 dvl["PHI0"] = phi0;
849 dvl["X0"] = (int_4)x0;
850 dvl["Y0"] = (int_4)y0;
851 dvl["ANGLE"] = angle;
852 dvl["ANGLEX"] = anglex;
853 dvl["ANGLEY"] = angley;
854 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
855 delete[] r_8tab_;}
856
857
858void FitsIoServer::save(const ImageR4& DpcImg,char flnm[])
859{
860 long naxis=2;
861 int siz_x = DpcImg.XSize();
862 int siz_y = DpcImg.YSize();
863 int nbpixels = siz_x*siz_y;
864 FITS_tab_typ_ = TFLOAT;
865
866
867
868 // write FITS image
869 DVList dvl;
870
871 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
872 dvl["ORG_X"] = DpcImg.XOrg();
873 dvl["ORG_Y"] = DpcImg.YOrg();
874 dvl["PIXSZ_X"] = DpcImg.XPxSize();
875 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
876 dvl["IDENT"] = DpcImg.Ident();
877
878 //dvl["NAME"] = DpcImg.Nom();
879
880 // j utilise la methode SetS parce que ses parametres sont const et
881 // que l'argument DpcImg est const dans la presente method.
882 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
883 // est non const et provoque un warning sur mac (CodeWarrior)
884 dvl.SetS("NAME",DpcImg.Nom());
885
886 dvl["NBSAT"] = DpcImg.nbSat;
887 dvl["NBNUL"] = DpcImg.nbNul;
888 dvl["MINPIX"] = DpcImg.minPix;
889 dvl["MAXPIX"] = DpcImg.maxPix;
890 dvl["MOYPIX"] = DpcImg.moyPix;
891 dvl["SIGPIX"] = DpcImg.sigPix;
892 dvl["FOND"] = DpcImg.fond;
893 dvl["SIGFON"] = DpcImg.sigmaFond;
894
895 // get the values of the DpcImage
896 if (r_4tab_ != NULL) delete [] r_4tab_;
897 r_4tab_=new r_4[siz_x*siz_y];
898 PBaseDataTypes dataT=DataType((r_4)0);
899 memcpy( r_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
900 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
901
902
903 delete [] r_4tab_;
904}
905
906
907void FitsIoServer::save(const ImageI4& DpcImg,char flnm[])
908{
909 long naxis=2;
910 int siz_x = DpcImg.XSize();
911 int siz_y = DpcImg.YSize();
912 int nbpixels = siz_x*siz_y;
913 FITS_tab_typ_ = TINT;
914
915 // pointer to the FITS file, defined in fitsio.h
916 //fitsfile *fptr;
917
918 // initialize status before calling fitsio routines
919 //int status = 0;
920
921 // delete old file if it already exists
922 //remove(flnm);
923
924 // create new FITS file
925 //fits_create_file(&fptr, flnm, &status);
926 //if( status ) printerror( status );
927
928
929 // write FITS image
930 DVList dvl;
931
932 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
933 dvl["ORG_X"] = DpcImg.XOrg();
934 dvl["ORG_Y"] = DpcImg.YOrg();
935 dvl["PIXSZ_X"] = DpcImg.XPxSize();
936 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
937 dvl["IDENT"] = DpcImg.Ident();
938
939 //dvl["NAME"] = DpcImg.Nom();
940 // j utilise la methode SetS parce que ses parametres sont const et
941 // que l'argument DpcImg est const dans la presente method.
942 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
943 // est non const et provoque un warning sur mac (CodeWarrior)
944 dvl.SetS("NAME",DpcImg.Nom());
945
946 dvl["NBSAT"] = DpcImg.nbSat;
947 dvl["NBNUL"] = DpcImg.nbNul;
948 dvl["MINPIX"] = DpcImg.minPix;
949 dvl["MAXPIX"] = DpcImg.maxPix;
950 dvl["MOYPIX"] = DpcImg.moyPix;
951 dvl["SIGPIX"] = DpcImg.sigPix;
952 dvl["FOND"] = DpcImg.fond;
953 dvl["SIGFON"] = DpcImg.sigmaFond;
954
955 // get the values of the DpcImage
956 if (i_4tab_ != NULL) delete [] i_4tab_;
957 i_4tab_=new int_4[siz_x*siz_y];
958 PBaseDataTypes dataT=DataType((int_4)0);
959 memcpy( i_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
960
961
962 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
963
964 // close the file
965 //fits_close_file(fptr, &status);
966 //if( status ) printerror( status );
967
968 delete [] i_4tab_;
969}
970
971
972
973
974void FitsIoServer::planck_write_img(char flnm[], int naxis,int n1, int n2, int n3, DVList& dvl)
975{
976
977
978 // pointer to the FITS file, defined in fitsio.h
979 fitsfile *fptr;
980
981 // initialize status before calling fitsio routines
982 int status = 0;
983
984 // delete old file if it already exists
985 remove(flnm);
986
987 // create new FITS file
988 fits_create_file(&fptr, flnm, &status);
989 if( status ) printerror( status );
990
991 int bitpix=0;
992 if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
993 if ( FITS_tab_typ_ == TFLOAT) bitpix = FLOAT_IMG;
994 if ( FITS_tab_typ_ == TINT) bitpix = LONG_IMG;
995 long naxes[3];
996 naxes[0] = n1;
997 naxes[1] = n2;
998 naxes[2] = n3;
999 if (n2 > 0 && naxis < 2) cout << " FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
1000 if (n3 > 0 && naxis < 3) cout << " FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
1001 fits_create_img(fptr, bitpix, naxis, naxes, &status);
1002 if( status ) printerror( status );
1003
1004
1005
1006 long nelements= naxes[0];
1007 if (naxis >=2) nelements*=naxes[1];
1008 if (naxis == 3) nelements*=naxes[2];
1009 switch (FITS_tab_typ_)
1010 {
1011 case TDOUBLE :
1012 fits_write_img(fptr, TDOUBLE, 1, nelements, r_8tab_, &status);
1013 if( status ) printerror( status );
1014 break;
1015 case TFLOAT :
1016 fits_write_img(fptr, TFLOAT, 1, nelements, r_4tab_, &status);
1017 if( status ) printerror( status );
1018 break;
1019 case TINT :
1020 fits_write_img(fptr, TINT, 1, nelements, i_4tab_, &status);
1021 if( status ) printerror( status );
1022 break;
1023 default :
1024 cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
1025 break;
1026 }
1027
1028 // write the current date
1029 fits_write_date(fptr, &status);
1030 if( status ) printerror( status );
1031 // on convient d ecrire apres la date, les mots cles utilisateur.
1032 // la date servira de repere pour relire ces mots cles.
1033
1034 //dvl.Print();
1035 DVList::ValList::const_iterator it;
1036 for(it = dvl.Begin(); it != dvl.End(); it++)
1037 {
1038 int datatype= key_type_PL2FITS( (*it).second.typ);
1039 char keyname[]="";
1040 int flen_keyword = 9;
1041 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
1042 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
1043 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
1044 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
1045 char comment[FLEN_COMMENT]="";
1046 int ival=0;
1047 double dval=0.;
1048 char strval[]="";
1049 switch (datatype)
1050 {
1051 case TINT :
1052 ival=(*it).second.mtv.iv;
1053 strcpy(comment,"I entier");
1054 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
1055 break;
1056 case TDOUBLE :
1057 dval=(*it).second.mtv.dv;
1058 strcpy(comment,"D double");
1059 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
1060 break;
1061 case TSTRING :
1062 strcpy(strval, (*it).second.mtv.strv);
1063 strcpy(comment,"S character string");
1064 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
1065 break;
1066 default :
1067 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
1068 break;
1069 }
1070 if( status ) printerror( status );
1071 }
1072 // close the file
1073 fits_close_file(fptr, &status);
1074 if( status ) printerror( status );
1075}
1076
1077void FitsIoServer::planck_write_bntbl(char flnm[], int npixels, char typeOfContent[], char extname[], char comment1[], DVList& dvl)
1078{
1079
1080
1081 // pointer to the FITS file, defined in fitsio.h
1082 fitsfile *fptr;
1083
1084 // initialize status before calling fitsio routines
1085 int status = 0;
1086
1087 // delete old file if it already exists
1088 remove(flnm);
1089
1090 // create new FITS file
1091 fits_create_file(&fptr, flnm, &status);
1092 if( status ) printerror( status );
1093
1094 // primary header
1095 int bitpix=LONG_IMG;
1096 int naxis=0;
1097 fits_create_img(fptr, bitpix, naxis, NULL, &status);
1098 // write the current date
1099 fits_write_date(fptr, &status);
1100 if( status ) printerror( status );
1101
1102
1103 long nrows=npixels/1024;
1104 if (nrows==0) nrows=1;
1105 if (npixels%1024 !=0) nrows++;
1106 int tfields=1;
1107 char **ttype, **tform;
1108 ttype= new char*[tfields];
1109 tform= new char*[tfields];
1110 char* format;
1111 if ( FITS_tab_typ_ == TDOUBLE) format="1024D";
1112 if ( FITS_tab_typ_ == TFLOAT) format="1024E";
1113 if ( FITS_tab_typ_ == TINT) format="1024I";
1114 for(int i = 0; i < tfields; i++)
1115 {
1116 ttype[i]= new char[FLEN_VALUE];
1117 strcpy(ttype[i], typeOfContent);
1118 tform[i]= new char[FLEN_VALUE];
1119 strcpy(tform[i], format);
1120 }
1121 // create a new empty binary table onto the FITS file
1122 // physical units if they exist, are defined in the DVList object
1123 // so the null pointer is given for the tunit parameters.
1124 fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
1125 NULL,extname,&status);
1126 if( status ) printerror( status );
1127 for( int ii = 0; ii < tfields; ii++)
1128 {
1129 delete [] ttype[ii];
1130 delete [] tform[ii];
1131 }
1132 delete [] ttype;
1133 delete [] tform;
1134 long nelements= npixels;
1135 switch (FITS_tab_typ_)
1136 {
1137 case TDOUBLE :
1138 fits_write_col(fptr, TDOUBLE, 1, 1, 1, nelements, r_8tab_, &status);
1139 if( status ) printerror( status, "probleme dans ecriture du tableau de doubles" );
1140 break;
1141
1142
1143 case TFLOAT :
1144 for (int kk=0; kk<10; kk++) cout << r_4tab_[kk] << endl;
1145 fits_write_col(fptr, TFLOAT, 1, 1, 1, nelements, r_4tab_, &status);
1146 if( status ) printerror( status );
1147 break;
1148 case TINT :
1149 fits_write_col(fptr, TINT, 1, 1, 1, nelements, i_4tab_, &status);
1150 if( status ) printerror( status );
1151 break;
1152 default :
1153 cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
1154 break;
1155 }
1156 // write supplementary keywords
1157 fits_write_comment(fptr, " ", &status);
1158 fits_write_comment(fptr, " ", &status);
1159 if( status ) printerror( status );
1160 fits_write_comment(fptr,"--------------------------------------------", &status);
1161 fits_write_comment(fptr, comment1, &status);
1162 fits_write_comment(fptr,"--------------------------------------------", &status);
1163 if( status ) printerror( status );
1164
1165
1166 int flen_keyword = 9;
1167 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
1168 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
1169 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
1170 DVList::ValList::const_iterator it;
1171 for(it = dvl.Begin(); it != dvl.End(); it++)
1172 {
1173 int datatype= key_type_PL2FITS( (*it).second.typ);
1174 char keyname[]="";
1175 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
1176 char comment[FLEN_COMMENT]="";
1177 int ival=0;
1178 double dval=0.;
1179 char strval[]="";
1180 switch (datatype)
1181 {
1182 case TINT :
1183 ival=(*it).second.mtv.iv;
1184 strcpy(comment," ");
1185 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
1186 break;
1187 case TDOUBLE :
1188 dval=(*it).second.mtv.dv;
1189 strcpy(comment," ");
1190 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
1191 break;
1192 case TSTRING :
1193 strcpy(strval, (*it).second.mtv.strv);
1194 strcpy(comment," ");
1195 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
1196 break;
1197 default :
1198 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
1199 break;
1200 }
1201 if( status ) printerror( status );
1202 }
1203 //char keyname[]="";
1204 //char strval[]="";
1205 //char comment[FLEN_COMMENT]="";
1206 //strncpy(keyname, "CREATOR",flen_keyword-1);
1207 //strcpy(strval, "SOPHYA");
1208 //strcpy(comment," Orsay");
1209 //fits_write_key(fptr, TSTRING, keyname, &strval, comment, &status);
1210 //if( status ) printerror( status );
1211 // close the file
1212 fits_close_file(fptr, &status);
1213 if( status ) printerror( status );
1214}
1215
1216void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
1217{
1218 int status=0;
1219 long bitpix;
1220 long naxes[3]={0,0,0};
1221 char* comment=NULL;
1222
1223 // pointer to the FITS file, defined in fitsio.h
1224 fitsfile *fptr;
1225 // initialize status before calling fitsio routines
1226 fits_open_file(&fptr, flnm, READONLY, &status);
1227 if( status ) printerror( status );
1228
1229
1230 fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
1231 if( status ) printerror( status );
1232 fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
1233 if( status ) printerror( status );
1234 int nfound;
1235 int nkeys=(int)naxis;
1236 fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
1237 if( status ) printerror( status );
1238
1239 n1 = naxes[0] ;
1240 n2 = naxes[1] ;
1241 n3 = naxes[2] ;
1242
1243
1244 long nelements= naxes[0];
1245 if (naxis >=2) nelements*=naxes[1];
1246 if (naxis == 3) nelements*=naxes[2];
1247 int anynull;
1248 r_8 dnullval=0.;
1249 r_4 fnullval=0.;
1250 int_4 inullval=0;
1251 // on laisse a fits le soin de convertir le type du tableau lu vers
1252 // le type suppose par l'utilisateur de fitsioserver
1253 //
1254 switch ( FITS_tab_typ_)
1255 {
1256 case TDOUBLE :
1257 if (bitpix != DOUBLE_IMG)
1258 {
1259 cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
1260 << " conversion to double will be achieved by cfitsio lib " << endl;
1261 }
1262 if (r_8tab_ != NULL) delete [] r_8tab_;
1263 r_8tab_=new r_8[nelements];
1264 fits_read_img(fptr, TDOUBLE, 1, nelements, &dnullval, r_8tab_,
1265 &anynull, &status);
1266 if( status ) printerror( status );
1267 break;
1268 case TFLOAT :
1269 if (bitpix != FLOAT_IMG)
1270 {
1271 cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
1272 << " conversion to float will be achieved by cfitsio lib " << endl;
1273 }
1274 if (r_4tab_ != NULL) delete [] r_4tab_;
1275 r_4tab_=new r_4[nelements];
1276 fits_read_img(fptr, TFLOAT, 1, nelements, &fnullval, r_4tab_,
1277 &anynull, &status);
1278 if( status ) printerror( status );
1279 break;
1280
1281
1282 case TINT :
1283 if (bitpix != LONG_IMG)
1284 {
1285 cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
1286 << " conversion to long will be achieved by cfitsio lib " << endl;
1287 }
1288 if (i_4tab_ != NULL) delete [] i_4tab_;
1289 i_4tab_=new int_4[nelements];
1290 fits_read_img(fptr, TINT, 1, nelements, &inullval, i_4tab_,
1291 &anynull, &status);
1292 if( status ) printerror( status );
1293 break;
1294
1295
1296 default :
1297 cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_ << endl;
1298 break;
1299 }
1300 status = 0;
1301 char card[FLEN_CARD];
1302 int num = 0;
1303 char comment2[FLEN_COMMENT] = "x";
1304 char keyname[]= "";
1305 char datekey[]= "DATE";
1306 char endkey[] = "END";
1307 char typ='x';
1308 int ival;
1309 double dval;
1310 char strval[]="";
1311 // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
1312 // on va jusqu'au mot cle DATE
1313 int flen_keyword = 9;
1314 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
1315 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
1316 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
1317 while (status == 0 && strncmp(keyname, datekey,4) != 0 )
1318 {
1319 num++;
1320 fits_read_record(fptr, num, card, &status);
1321 strncpy(keyname,card,flen_keyword-1);
1322 }
1323 if (status != 0 )
1324 {
1325 cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
1326 }
1327 // on recupere la liste des mots-cles utilisateurs
1328 while (status == 0)
1329 {
1330 num++;
1331 // on lit un record pour recuperer le nom du mot-cle
1332 fits_read_record(fptr, num, card, &status);
1333 strncpy(keyname,card,flen_keyword-1);
1334 char value[FLEN_VALUE];
1335 // on recupere le premier caractere du commentaire, qui contient
1336 // le code du type de la valeur
1337 // (tant que ce n est pas le mot cle END)
1338 fits_read_keyword(fptr, keyname, value, comment2, &status);
1339 if ( strncmp(keyname, endkey,flen_keyword-1) != 0)
1340 {
1341 typ = comment2[0];
1342 // quand le type est connu, on lit la valeur
1343 switch (typ)
1344 {
1345 case 'I' :
1346 fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
1347 if( status ) printerror( status );
1348 strip (keyname, 'B',' ');
1349 dvl[keyname] = (int_4)ival;
1350 break;
1351 case 'D' :
1352 fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
1353 if( status ) printerror( status );
1354 strip (keyname, 'B',' ');
1355 dvl[keyname] = dval;
1356 break;
1357 case 'S' :
1358 fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
1359 if( status ) printerror( status );
1360 strip (keyname, 'B',' ');
1361 strip(strval, 'B',' ');
1362 dvl[keyname]=strval;
1363 break;
1364 default :
1365 cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
1366 break;
1367 }
1368 }
1369 }
1370
1371
1372 // close the file
1373 status=0;
1374 fits_close_file(fptr, &status);
1375 if( status ) printerror( status );
1376}
1377
1378
1379void FitsIoServer::planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl)
1380{
1381 int status=0;
1382 int nkeys,keypos;
1383 int hdutype;
1384 int tfields;
1385 int datype;
1386 long lastpix;
1387 long repeat, width;
1388 long nrows;
1389 long extend;
1390 char* comment=NULL;
1391
1392 // pointer to the FITS file, defined in fitsio.h
1393 fitsfile *fptr;
1394 // initialize status before calling fitsio routines
1395 fits_open_file(&fptr, flnm, READONLY, &status);
1396 if( status ) printerror( status );
1397 fits_read_key_lng(fptr, "EXTEND", &extend, comment, &status);
1398 if( status ) printerror( status );
1399 if (extend!=1)
1400 {
1401 cout << "FitsIoServer:: le fichier fits ne contient pas d'extension binary table" << endl;
1402 return;
1403 }
1404 fits_movabs_hdu(fptr, hdunum,&hdutype,&status);
1405 if( status ) printerror( status );
1406 if (hdutype!=BINARY_TBL)
1407 {
1408 cout << "FitsIoServer:: this HDU is not a binary table " << endl;
1409 exit(status);
1410 }
1411 char xtension[FLEN_VALUE];
1412 fits_read_key_str(fptr,"XTENSION",xtension,NULL,&status);
1413 if( status ) printerror( status );
1414
1415 char binta[] = "BINTABLE";
1416 if ( strncmp(xtension, binta,8) != 0)
1417 // if (xtension !="BINTABLE")
1418 {
1419 cout << "FitsIoServer:: not a binary table " << endl;
1420 exit(status);
1421 }
1422 fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
1423 if( status ) printerror( status );
1424 //cout << " nombre de mots-cles : " << nkeys << endl;
1425 fits_get_num_cols(fptr, &tfields, &status);
1426 if (tfields != 1)
1427 {
1428 cout << "FitsIoServer:: il y a plus d'une colonne" << endl;
1429 return;
1430 }
1431 fits_get_num_rows(fptr, &nrows, &status);
1432 //cout << "nblignes= " << nrows << endl;
1433 fits_get_coltype(fptr, 1, &datype, &repeat, &width, &status);
1434 if( status ) printerror( status );
1435 //cout << " type de donnees : " << datype << endl;
1436 //cout << " repeat : " << repeat << endl;
1437 //cout << " width : " << width << endl;
1438 fits_read_key_lng(fptr, "LASTPIX", &lastpix, comment, &status);
1439 if( status ) printerror( status," mot cle LASTPIX" );
1440
1441 long nelements= nrows*repeat;
1442 if (nelements!=lastpix+1)
1443 {
1444 cout << " erreur sur longueur du vecteur " << endl;
1445 cout << " nelements= " << nelements << " lastpix+1=" << lastpix+1 << endl;
1446 }
1447 npixels=nelements;
1448 int anynull;
1449 r_8 dnullval=0.;
1450 r_4 fnullval=0.;
1451 int_4 inullval=0;
1452 // on laisse a fits le soin de convertir le type du tableau lu vers
1453 // le type suppose par l'utilisateur de fitsioserver
1454 //
1455 switch ( FITS_tab_typ_)
1456 {
1457 case TDOUBLE :
1458 if (datype != TDOUBLE)
1459 {
1460 cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
1461 << " conversion to double will be achieved by cfitsio lib " << endl;
1462 }
1463 if (r_8tab_ != NULL) delete [] r_8tab_;
1464 r_8tab_=new r_8[ npixels];
1465 fits_read_col(fptr, TDOUBLE, 1, 1, 1, nelements, &dnullval,
1466 r_8tab_,
1467 &anynull, &status);
1468 if( status ) printerror( status, "probleme dans lecture du tableau de doubles" );
1469 break;
1470 case TFLOAT :
1471 if (datype != TFLOAT)
1472 {
1473
1474 cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
1475 << " conversion to float will be achieved by cfitsio lib " << endl;
1476 }
1477 if (r_4tab_ != NULL) delete [] r_4tab_;
1478 r_4tab_=new r_4[nelements];
1479 fits_read_col(fptr, TFLOAT, 1, 1, 1, nelements, &fnullval,
1480 r_4tab_, &anynull, &status);
1481 if( status ) printerror( status,"probleme dans lecture du tableau de floats" );
1482 break;
1483
1484
1485 case TINT :
1486 if (datype != TLONG)
1487 {
1488 cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
1489 << " conversion to long will be achieved by cfitsio lib " << endl;
1490 }
1491 if (i_4tab_ != NULL) delete [] i_4tab_;
1492 i_4tab_=new int_4[nelements];
1493 fits_read_col(fptr, TLONG, 1, 1, 1, nelements, &inullval,
1494 i_4tab_, &anynull, &status);
1495 //fits_read_img(fptr, TINT, 1, nelements, &inullval, i_4tab_,
1496 // &anynull, &status);
1497 if( status ) printerror( status,"probleme dans lecture du tableau de ints" );
1498 break;
1499
1500
1501 default :
1502 cout << " FitsIOServer::read_bntbl : type non traite: " << FITS_tab_typ_ << endl;
1503 break;
1504 }
1505 char card[FLEN_CARD];
1506 char keyname[LEN_KEYWORD]= "";
1507 char strval[FLEN_VALUE];
1508 char comment1[FLEN_COMMENT];
1509 char dtype;
1510 //char bidon[LEN_KEYWORD];
1511 char comkey[] = "COMMENT";
1512
1513 for(int j = 1; j <= nkeys; j++)
1514 {
1515 // fits_read_record(fptr, j, card, &status);
1516 // strncpy(keyname,card,LEN_KEYWORD-1);
1517 // cout << " bidon= " << keyname << endl;
1518 // if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
1519 fits_read_keyn(fptr,j,card,strval,comment1,&status);
1520 strncpy(keyname,card,LEN_KEYWORD-1);
1521
1522 if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
1523 {
1524 fits_get_keytype(strval,&dtype,&status);
1525 // cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl;
1526 strip (keyname, 'B',' ');
1527 strip(strval, 'B',' ');
1528 switch( dtype )
1529 {
1530 case 'C':
1531 dvl[keyname]= strval;
1532 break;
1533 case 'I':
1534 int ival;
1535 ctoi(strval,&ival);
1536 dvl[keyname]= (int_4)ival;
1537 break;
1538 case 'F':
1539 double dval;
1540 ctof(strval,&dval);
1541 dvl[keyname]= dval;
1542 break;
1543 default :
1544 cout << " FitsIOServer : mot-cle bizarre " << endl;
1545 break;
1546 }
1547 }
1548 }
1549
1550 // close the file
1551 status=0;
1552 fits_close_file(fptr, &status);
1553 if( status ) printerror( status );
1554}
1555
1556
1557// projects a SphericalMap<double>, according sinus-method, and saves onto
1558// a FITS-file
1559void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
1560{
1561
1562 long naxes[2]={600, 300};
1563 float* map =new float[ 600*300 ];
1564 int npixels= naxes[0]*naxes[1];
1565
1566 cout << " image FITS en projection SINUS" << endl;
1567 // table will have npixels rows
1568 for(int j=0; j < npixels; j++) map[j]=0.;
1569 for(int j=0; j<naxes[1]; j++)
1570 {
1571 double yd = (j+0.5)/naxes[1]-0.5;
1572 double theta = (0.5-yd)*Pi;
1573 double facteur=1./sin(theta);
1574 for(int i=0; i<naxes[0]; i++)
1575 {
1576 double xa = (i+0.5)/naxes[0]-0.5;
1577 double phi = 2.*Pi*xa*facteur+Pi;
1578 float th=float(theta);
1579 float ff=float(phi);
1580 if (phi<2*Pi && phi>= 0)
1581 {
1582 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1583 }
1584 }
1585 }
1586
1587 write_picture(naxes, map, filename);
1588 delete [] map;
1589}
1590
1591// projects a SphericalMap<double>, according sinus-method, and saves onto
1592// a FITS-file
1593void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
1594{
1595 // le code de cete methode duplique celui de la precedente, la seule
1596 //difference etant le type de sphere en entree. Ces methodes de projection
1597 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1598 // ne me suis pas casse la tete, pour l instant
1599
1600 long naxes[2]={600, 300};
1601 float* map = new float[ 600*300 ];
1602 int npixels= naxes[0]*naxes[1];
1603
1604 cout << " image FITS en projection SINUS" << endl;
1605 // table will have npixels rows
1606 for(int j=0; j < npixels; j++) map[j]=0.;
1607 for(int j=0; j<naxes[1]; j++)
1608 {
1609 double yd = (j+0.5)/naxes[1]-0.5;
1610 double theta = (0.5-yd)*Pi;
1611 double facteur=1./sin(theta);
1612 for(int i=0; i<naxes[0]; i++)
1613 {
1614 double xa = (i+0.5)/naxes[0]-0.5;
1615 double phi = 2.*Pi*xa*facteur+Pi;
1616 float th=float(theta);
1617 float ff=float(phi);
1618 if (phi<2*Pi && phi>= 0)
1619 {
1620 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1621 }
1622 }
1623 }
1624
1625 write_picture(naxes, map, filename);
1626 delete [] map;
1627
1628}
1629
1630// projects a SphericalMap<float>, according Mollweide-method, and saves onto
1631// a FITS-file
1632void FitsIoServer::Mollweide_picture_projection(SphericalMap<float>& sph, char filename[])
1633{
1634 // le code de cete methode duplique celui de la precedente, la seule
1635 //difference etant le type de sphere en entree. Ces methodes de projection
1636 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1637 // ne me suis pas casse la tete, pour l instant
1638
1639 long naxes[2]={600, 300};
1640 float* map = new float[ 600*300 ];
1641 int npixels= naxes[0]*naxes[1];
1642
1643 cout << " image FITS en projection MOLLWEIDE" << endl;
1644 // table will have npixels rows
1645 for(int j=0; j < npixels; j++) map[j]=0.;
1646 for(int j=0; j<naxes[1]; j++)
1647 {
1648 double yd = (j+0.5)/naxes[1]-0.5;
1649 double facteur=2.*Pi/sin(acos(yd*2));
1650 double theta = (0.5-yd)*Pi;
1651 for(int i=0; i<naxes[0]; i++)
1652 {
1653 double xa = (i+0.5)/naxes[0]-0.5;
1654 double phi = xa*facteur+Pi;
1655 float th=float(theta);
1656 float ff=float(phi);
1657 if (phi<2*Pi && phi>= 0)
1658 {
1659 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1660 }
1661 }
1662 }
1663
1664 write_picture(naxes, map, filename);
1665 delete [] map;
1666
1667}
1668// projects a SphericalMap<double>, according Mollweide-method, and saves onto
1669// a FITS-file
1670void FitsIoServer::Mollweide_picture_projection(SphericalMap<double>& sph, char filename[])
1671{
1672 // le code de cete methode duplique celui de la precedente, la seule
1673 //difference etant le type de sphere en entree. Ces methodes de projection
1674 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1675 // ne me suis pas casse la tete, pour l instant
1676
1677 long naxes[2]={600, 300};
1678 float* map = new float[ 600*300 ];
1679 int npixels= naxes[0]*naxes[1];
1680
1681 cout << " image FITS en projection MOLLWEIDE" << endl;
1682 // table will have npixels rows
1683 for(int j=0; j < npixels; j++) map[j]=0.;
1684 for(int j=0; j<naxes[1]; j++)
1685 {
1686 double yd = (j+0.5)/naxes[1]-0.5;
1687 double facteur=2.*Pi/sin(acos(yd*2));
1688 double theta = (0.5-yd)*Pi;
1689 for(int i=0; i<naxes[0]; i++)
1690 {
1691 double xa = (i+0.5)/naxes[0]-0.5;
1692 double phi = xa*facteur+Pi;
1693 float th=float(theta);
1694 float ff=float(phi);
1695 if (phi<2*Pi && phi>= 0)
1696 {
1697 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1698 }
1699 }
1700 }
1701
1702 write_picture(naxes, map, filename);
1703 delete [] map;
1704
1705}
1706
1707
1708
1709// saves a (LocalMap<double> on a FITS-file in order to be visualized
1710// (for tests)
1711void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
1712{
1713
1714 long naxes[2];
1715 naxes[0] = lcm.Size_x();
1716 naxes[1] = lcm.Size_y();
1717 int npixels= naxes[0]*naxes[1];
1718 float* map = new float[npixels];
1719
1720 // table will have npixels rows
1721 for(int j=0; j < npixels; j++) map[j]=0.;
1722 for(int j=0; j<naxes[1]; j++)
1723 {
1724 for(int i=0; i<naxes[0]; i++)
1725 {
1726 map[j*naxes[0]+i] = lcm(i, j);
1727 }
1728 }
1729
1730 write_picture(naxes, map, filename);
1731 delete [] map;
1732}
1733
1734
1735
1736void FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
1737{
1738
1739 int bitpix = FLOAT_IMG;
1740 long naxis = 2;
1741
1742 //pointer to the FITS file, defined in fitsio.h
1743 fitsfile *fptr;
1744 // delete old file if it already exists
1745 remove(filename);
1746 // initialize status before calling fitsio routines
1747 int status = 0;
1748
1749 // create new FITS file
1750 fits_create_file(&fptr, filename, &status);
1751 if( status ) printerror( status );
1752
1753 // write the required header keywords
1754 fits_create_img(fptr, bitpix, naxis, naxes, &status);
1755 if( status ) printerror( status );
1756
1757 // write the current date
1758 fits_write_date(fptr, &status);
1759 if( status ) printerror( status );
1760
1761
1762 // first row in table to write
1763 long firstrow = 1;
1764 // first element in row
1765 long firstelem = 1;
1766 int colnum = 1;
1767 int nelements=naxes[0]*naxes[1];
1768 fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
1769 if( status ) printerror( status );
1770
1771 // close the file
1772 fits_close_file(fptr, &status);
1773 if( status ) printerror( status );
1774}
1775
1776
1777bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
1778
1779 //*****************************************************/
1780 //* check if the specified keyword exits in the CHU */
1781 //*****************************************************/
1782{
1783
1784 bool KEY_EXIST = false;
1785 int status = 0;
1786 char strbide[FLEN_VALUE];
1787 char keybide[LEN_KEYWORD]= "";
1788 for(int jj = 1; jj <= nkeys; jj++)
1789 {
1790 if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
1791 printerror( status );
1792 if( !strcmp(keybide,keyword) )
1793 {
1794 KEY_EXIST= true;
1795 break;
1796 }
1797 }
1798 return(KEY_EXIST);
1799}
1800
1801void FitsIoServer::readheader ( char filename[] )
1802
1803 //**********************************************************************/
1804 //* Print out all the header keywords in all extensions of a FITS file */
1805 //**********************************************************************/
1806{
1807
1808 // standard string lengths defined in fitsioc.h
1809 char card[FLEN_CARD];
1810
1811 // pointer to the FITS file, defined in fitsio.h
1812 fitsfile *fptr;
1813
1814 int status = 0;
1815 if ( fits_open_file(&fptr, filename, READONLY, &status) )
1816 printerror( status );
1817
1818 // attempt to move to next HDU, until we get an EOF error
1819 int hdutype;
1820 for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
1821 {
1822 if (hdutype == ASCII_TBL)
1823 printf("\nReading ASCII table in HDU %d:\n", ii);
1824 else if (hdutype == BINARY_TBL)
1825 printf("\nReading binary table in HDU %d:\n", ii);
1826 else if (hdutype == IMAGE_HDU)
1827 printf("\nReading FITS image in HDU %d:\n", ii);
1828 else
1829 {
1830 printf("Error: unknown type of this HDU \n");
1831 printerror( status );
1832 }
1833
1834 // get the number of keywords
1835 int nkeys, keypos;
1836 if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
1837 printerror( status );
1838
1839 printf("Header listing for HDU #%d:\n", ii);
1840 for (int jj = 1; jj <= nkeys; jj++)
1841 {
1842 if ( fits_read_record(fptr, jj, card, &status) )
1843 printerror( status );
1844
1845 // print the keyword card
1846 printf("%s\n", card);
1847 }
1848 printf("END\n\n");
1849 }
1850
1851 // got the expected EOF error; reset = 0
1852 if (status == END_OF_FILE)
1853 status = 0;
1854 else
1855 printerror( status );
1856
1857 if ( fits_close_file(fptr, &status) )
1858 printerror( status );
1859
1860 return;
1861}
1862
1863void FitsIoServer::printerror(int& status) const
1864
1865 //*****************************************************/
1866 //* Print out cfitsio error messages and exit program */
1867 //*****************************************************/
1868{
1869
1870 // print out cfitsio error messages and exit program
1871 // print error report
1872 fits_report_error(stderr, status);
1873 // terminate the program, returning error status
1874 // exit( status );
1875 status=0;
1876}
1877void FitsIoServer::printerror(int& status, char* texte) const
1878
1879 //*****************************************************/
1880 //* Print out cfitsio error messages and exit program */
1881 //*****************************************************/
1882{
1883
1884 // print out cfitsio error messages and exit program
1885 // print error report
1886 fits_report_error(stderr, status);
1887 cout << " erreur : " << texte << endl;
1888 status=0;
1889}
1890
1891
1892
Note: See TracBrowser for help on using the repository browser.