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

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

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