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

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

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