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

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

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