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

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

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