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

Last change on this file since 903 was 854, checked in by ansari, 25 years ago

Adaptation s SphereHEALPix et deplacement SkyMap/fitsspherehealpix ici

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