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

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

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