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

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

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