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

Last change on this file since 1004 was 923, checked in by ansari, 25 years ago

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