Ignore:
Timestamp:
Oct 20, 1999, 3:30:37 PM (26 years ago)
Author:
ansari
Message:

dans meth. sinus_picture_projection, declaration map=new... au lieu de map[...] 20-OCT-99 GLM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/FitsIOServer/fitsioserver.cc

    r479 r482  
     1//************************************************************************
     2// Class for loadind and saving from FITS-formatted file to DPC objects
     3//
     4// methods 'load(X& x, char f[])' get from FITS file "f" a DPC object x
     5// from DPC  class X.
     6// methods 'save(X& x, char f[])' save  a DPC object x from DPC // class X
     7// onto a FITS file "f" .
     8
     9//************************************************************************
    110
    211#include <iostream.h>
    … …  
    1120  int nbcols=0;
    1221  FITS_tab_typ_ = TDOUBLE;
    13   int status=0;
    1422  long naxis;
    1523  int n1, n2, n3;
    1624  DVList dvl;
    17   // pointer to the FITS file, defined in fitsio.h
    18   fitsfile *fptr;
    19   fits_open_file(&fptr, flnm, READONLY, &status);
    20   if( status )  printerror( status );
    21   planck_read_img(fptr, naxis, n1, n2, n3, dvl);
    22   // close the file
    23   fits_close_file(fptr, &status);
    24   if( status )  printerror( status );
     25  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
    2526
    2627  nbrows=n1;
    … …  
    7374    for (int i = 0; i < nbrows; i++)  mat(i,j) = dtab_[ij++];
    7475}
    75 
    76 void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
    77 {
    78   int npixels=0;
    79   int nside=0;
    80   FITS_tab_typ_ = TDOUBLE;
    81   read_sphere(flnm, npixels, nside);
    82 
    83   // number of pixels in the sphere
    84   if (sph.NbPixels() != npixels)
    85     {
    86       cout << " found " << npixels << " pixels" << endl;
    87       cout << " expected " << sph.NbPixels() << endl;
    88       if (nside <= 0 )
    89         {
    90           cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
    91           exit(0);
    92         }
    93       if (nside != sph.SizeIndex())
    94         {
    95           sph.Resize(nside);
    96           cout << " resizing the sphere to nside=  " << nside  << endl;
    97         }
    98       else
    99         {
    100           cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
    101           exit(0);
    102         }
    103     }
    104   for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
    105 }
    106 void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
    107 {
    108   int npixels=0;
    109   int nside=0;
    110   FITS_tab_typ_ = TFLOAT;
    111   read_sphere(flnm, npixels, nside);
    112   // number of pixels in the sphere
    113   if (sph.NbPixels() != npixels)
    114     {
    115       cout << " found " << npixels << " pixels" << endl;
    116       cout << " expected " << sph.NbPixels() << endl;
    117       if (nside <= 0 )
    118         {
    119           cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
    120           exit(0);
    121         }
    122       if (nside != sph.SizeIndex())
    123         {
    124           sph.Resize(nside);
    125           cout << " resizing the sphere to nside=  " << nside  << endl;
    126         }
    127       else
    128         {
    129           cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
    130           exit(0);
    131         }
    132     }
    133   for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
    134 }
    135 
    136 void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
    137 {
    138   int nside;
    139   int nbrows=0;
    140   int nbcols=0;
    141   FITS_tab_typ_ = TDOUBLE;
    142   int status=0;
    143   long naxis;
    144   int n1, n2, n3;
    145   DVList dvl;
    146   // pointer to the FITS file, defined in fitsio.h
    147   fitsfile *fptr;
    148   fits_open_file(&fptr, flnm, READONLY, &status);
    149   if( status )  printerror( status );
    150   planck_read_img(fptr, naxis, n1, n2, n3, dvl);
    151   // close the file
    152   fits_close_file(fptr, &status);
    153   if( status )  printerror( status );
    154 
    155   nbrows=n1;
    156   nbcols=n2;
    157   if (naxis != 2)
    158     {
    159       cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
    160     }
    161   //dvl.Print();
    162   DVList::ValList::const_iterator it;
    163   for(it = dvl.Begin(); it != dvl.End(); it++) 
    164     {
    165       char datatype= (*it).second.typ;
    166       char keyname[]="";
    167       strcpy(keyname, (*it).first.substr(0,64).c_str());   
    168       int ival=0;
    169       double dval=0.;
    170       char strval[]="";
    171       switch (datatype)
    172         {
    173         case 'I' :
    174           ival=(*it).second.mtv.iv;
    175           break;
    176         case 'D' :
    177           dval=(*it).second.mtv.dv;
    178           break;
    179         case 'S' :
    180           strcpy(strval, (*it).second.mtv.strv);
    181           break;
    182         default :
    183           cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    184           break;
    185         }
    186       if (!strcmp(keyname, "NSIDE") )
    187         {
    188           nside = ival;
    189         }
    190       //if (!strcmp(keyname, "ORDERING") )
    191       //        {
    192       //  cout << " j'ai trouve ORDERING " << endl;
    193       //        }
    194     }
    195   float theta0 = dvl.GetD("THETA0");
    196   float phi0 = dvl.GetD("PHI0");
    197   int x0 = dvl.GetI("X0");
    198   int y0 = dvl.GetI("Y0");
    199   int xo = dvl.GetI("XO");
    200   int yo = dvl.GetI("YO");
    201   float anglex=dvl.GetD("ANGLEX");
    202   float angley=dvl.GetD("ANGLEY");
    203   float angle = dvl.GetD("ANGLE");
    204 
    205   // number of components
    206 
    207   if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
    208     {
    209       cout << " found " << nbrows << " x-pixels ";
    210       cout << "  expected " << lcm.Size_x()  << endl;
    211       cout << " found " << nbcols << " y-pixels " ;
    212       cout << "  expected " <<   lcm.Size_y()  << endl;
    213       lcm.ReSize(nbrows,nbcols);
    214       cout << " resizing the map to x-size=  " << nbrows  << " y-size= " << nbcols << endl;
    215     }
    216   lcm.SetSize(anglex, angley);
    217   lcm.SetOrigin(theta0, phi0, x0, y0, angle);
    218   int ij=0;
    219   for (int j=0; j< nbcols; j++)
    220     for (int i = 0; i < nbrows; i++)  lcm(i,j) = dtab_[ij++];
    221 }
    222 
    223 
    224 void FitsIoServer::save( TMatrix<double>& mat, char filename[])   
    225 {
    226   int nbrows = mat.NRows();
    227   int nbcols = mat.NCols();
    228   long naxis = nbcols > 1 ? 2 : 1;
    229   cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
    230   FITS_tab_typ_ = TDOUBLE;
    231   if (dtab_ != NULL ) delete[] dtab_;
    232   dtab_=new double[nbrows*nbcols];
    233 
    234 
    235   //pointer to the FITS file, defined in fitsio.h
    236   fitsfile *fptr;     
    237  
    238  
    239   // initialize status before calling fitsio routines
    240   int status = 0;         
    241 
    242   // delete old file if it already exists
    243   remove(filename);
    244 
    245   // create new FITS file 
    246   fits_create_file(&fptr, filename, &status);
    247   if( status )  printerror( status );
    248   int ij=0;
    249   for (int j=0; j< nbcols; j++)
    250     for (int i = 0; i < nbrows; i++) dtab_[ij++]= mat(i,j);
    251 
    252   DVList dvl;
    253 
    254   planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
    255 
    256 
    257   // close the file
    258   fits_close_file(fptr, &status);
    259   if( status )  printerror( status );
    260 
    261 }
    262 
    263 //void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
    264 void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
    265 {
    266   int npixels = sph.NbPixels();
    267   FITS_tab_typ_ = TDOUBLE;
    268   if (dtab_ != NULL ) delete[] dtab_;
    269   dtab_=new double[npixels];
    270 
    271   //pointer to the FITS file, defined in fitsio.h
    272   fitsfile *fptr;     
    273  
    274  
    275   // initialize status before calling fitsio routines
    276   int status = 0;         
    277 
    278   // delete old file if it already exists
    279   remove(filename);
    280 
    281   // create new FITS file 
    282   fits_create_file(&fptr, filename, &status);
    283   if( status )  printerror( status );
    284 
    285   for (int j = 0; j < npixels; j++) dtab_[j]= sph(j);
    286   DVList dvl;
    287   dvl["NSIDE"] = sph.SizeIndex();
    288   dvl["ORDERING"]=sph.TypeOfMap();
    289 
    290   planck_write_img(fptr, 1, npixels, 0, 0, dvl);
    291 
    292   // decider ulterieurement ce qu'on fait de ce qui suit, specifique
    293   // pour l'instant, aux spheres gorski. Actuellement, on ne sauve pas
    294   // les informations d'analyse harmonique
    295   /*
    296   // write supplementary keywords
    297   fits_write_comment(fptr, " ", &status);
    298   if( status ) printerror( status );
    299 
    300 
    301   strcpy(comment,"HEALPIX Pixelisation");
    302   strcpy(svalue, "HEALPIX");
    303   fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
    304   if( status )  printerror( status );
    305 
    306 
    307   strcpy(comment,"pixel ordering scheme, either RING or NESTED");
    308   strcpy(svalue, "RING");
    309   fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
    310   if( status )  printerror( status );
    311 
    312 
    313   strcpy(comment,"Random generator seed");
    314   int iseed= sph.iseed();
    315   fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
    316   if( status )  printerror( status );
    317 
    318 
    319   strcpy(comment,"--------------------------------------------");
    320   fits_write_comment(fptr, comment, &status);
    321   if( status ) printerror( status );
    322 
    323   strcpy(comment," Above keywords are still likely to change !");
    324   fits_write_comment(fptr, comment, &status);
    325   if( status ) printerror( status );
    326 
    327   strcpy(comment,"--------------------------------------------");
    328   fits_write_comment(fptr, comment, &status);
    329   if( status ) printerror( status );
    330  
    331   */
    332 
    333   // close the file
    334   fits_close_file(fptr, &status);
    335   if( status )  printerror( status );
    336 
    337 }
    338 
    339 void FitsIoServer::save(SphericalMap<float>& sph, char filename[])   
    340 {
    341   int npixels = sph.NbPixels();
    342   FITS_tab_typ_ = TFLOAT;
    343   if (ftab_ != NULL ) delete[] ftab_;
    344   ftab_=new float[npixels];
    345 
    346 
    347   //pointer to the FITS file, defined in fitsio.h
    348   fitsfile *fptr;     
    349  
    350  
    351   // initialize status before calling fitsio routines
    352   int status = 0;         
    353 
    354   // delete old file if it already exists
    355   remove(filename);
    356 
    357   // create new FITS file 
    358   fits_create_file(&fptr, filename, &status);
    359   if( status )  printerror( status );
    360   for (int j = 0; j < npixels; j++) ftab_[j]= sph(j);
    361   DVList dvl;
    362   dvl["NSIDE"] = sph.SizeIndex();
    363   dvl["ORDERING"]=sph.TypeOfMap();
    364 
    365   planck_write_img(fptr, 1, npixels, 0, 0, dvl);
    366 
    367 
    368   // close the file
    369   fits_close_file(fptr, &status);
    370   if( status )  printerror( status );
    371 
    372 }
    373 
    374 void FitsIoServer::save(LocalMap<double>& locm, char filename[])   
    375 {
    376   int nbrows = locm.Size_x();
    377   int nbcols = locm.Size_y();
    378   long naxis = 2;
    379   cout << " nombre de pts en x : " << nbrows << " en y  " << nbcols << endl;
    380   FITS_tab_typ_ = TDOUBLE;
    381   if (dtab_ != NULL ) delete[] dtab_;
    382   dtab_=new double[nbrows*nbcols];
    383 
    384   //pointer to the FITS file, defined in fitsio.h
    385   fitsfile *fptr;     
    386  
    387  
    388   // initialize status before calling fitsio routines
    389   int status = 0;         
    390 
    391   // delete old file if it already exists
    392   remove(filename);
    393 
    394   // create new FITS file 
    395   fits_create_file(&fptr, filename, &status);
    396   if( status )  printerror( status );
    397   int ij=0;
    398   for (int j=0; j< nbcols; j++)
    399     for (int i = 0; i < nbrows; i++) dtab_[ij++]= locm(i,j);
    400 
    401   DVList dvl;
    402   dvl["NSIDE"] = locm.SizeIndex();
    403   dvl["ORDERING"]=locm.TypeOfMap();
    404   double theta0;
    405   double phi0;
    406   int x0;
    407   int y0;
    408   double angle;
    409   locm.Origin(theta0,phi0,x0,y0,angle);
    410   double anglex;
    411   double angley;
    412   locm.Aperture(anglex,angley);
    413   dvl["THETA0"] = theta0;
    414   dvl["PHI0"] = phi0;
    415   dvl["X0"] = (int_4)x0;
    416   dvl["Y0"] = (int_4)y0;
    417   dvl["ANGLE"] = angle;
    418   dvl["ANGLEX"] = anglex;
    419   dvl["ANGLEY"] = angley;
    420   planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
    421 
    422 
    423   // close the file
    424   fits_close_file(fptr, &status);
    425   if( status )  printerror( status );
    426 
    427 }
    428 void FitsIoServer::planck_write_img( fitsfile *fptr, int naxis,int n1, int n2, int n3, DVList& dvl)
    429 {
    430 
    431   int status=0;
    432   int bitpix=0;
    433   if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
    434   if ( FITS_tab_typ_ == TFLOAT)  bitpix = FLOAT_IMG;
    435   if ( FITS_tab_typ_ == TINT)   bitpix = LONG_IMG;
    436   long naxes[3];
    437   naxes[0] = n1;
    438   naxes[1] = n2;
    439   naxes[2] = n3;
    440   if (n2 > 0 && naxis < 2) cout << "  FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
    441   if (n3 > 0 && naxis < 3) cout << "  FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
    442   fits_create_img(fptr, bitpix, naxis, naxes, &status);
    443   if( status )  printerror( status );
    444 
    445 
    446 
    447   long nelements= naxes[0];
    448   if (naxis >=2) nelements*=naxes[1];
    449   if (naxis == 3)  nelements*=naxes[2];
    450   switch (FITS_tab_typ_)
    451     {
    452     case TDOUBLE :
    453       fits_write_img(fptr, TDOUBLE, 1, nelements, dtab_,  &status);
    454       if( status ) printerror( status );
    455       break;
    456     case TFLOAT :
    457       fits_write_img(fptr, TFLOAT, 1, nelements, ftab_,  &status);
    458       if( status ) printerror( status );
    459       break;
    460     case TINT :
    461       fits_write_img(fptr, TINT, 1, nelements, i_4tab_,  &status);
    462       if( status ) printerror( status );
    463       break;
    464     default :
    465       cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
    466       break;
    467     }
    468 
    469   // write the current date
    470   fits_write_date(fptr, &status);
    471   if( status )  printerror( status );
    472   // on convient d ecrire apres la date, les mots cles utilisateur.
    473   // la date servira de repere pour relire ces mots cles.
    474 
    475   //dvl.Print();
    476   DVList::ValList::const_iterator it;
    477   for(it = dvl.Begin(); it != dvl.End(); it++) 
    478     {
    479       int datatype=  key_type_PL2FITS( (*it).second.typ);
    480       char keyname[]="";
    481       int flen_keyword = 9;
    482       // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
    483       // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
    484       // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
    485       strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
    486       char comment[FLEN_COMMENT]="";
    487       int ival=0;
    488       double dval=0.;
    489       char strval[]="";
    490       switch (datatype)
    491         {
    492         case TINT :
    493           ival=(*it).second.mtv.iv;
    494           strcpy(comment,"I entier");
    495           fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
    496           break;
    497         case TDOUBLE :
    498           dval=(*it).second.mtv.dv;
    499           strcpy(comment,"D double");
    500           fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
    501           break;
    502         case TSTRING :
    503           strcpy(strval, (*it).second.mtv.strv);
    504           strcpy(comment,"S character string");
    505           fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
    506           break;
    507         default :
    508           cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    509           break;
    510         }
    511       if( status ) printerror( status );
    512     }
    513 }
    514 
    515 void FitsIoServer::planck_read_img( fitsfile *fptr, long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
    516 {
    517   int status=0;
    518   long bitpix;
    519   long naxes[3]={0,0,0};
    520   char* comment=NULL;
    521 
    522   fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
    523   if( status )  printerror( status );
    524   fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
    525   if( status )  printerror( status );
    526   int nfound;
    527   int nkeys=(int)naxis;
    528   fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
    529   if( status )  printerror( status );
    530 
    531   n1 = naxes[0] ;
    532   n2 = naxes[1] ;
    533   n3 = naxes[2] ;
    534 
    535 
    536   long nelements= naxes[0];
    537   if (naxis >=2) nelements*=naxes[1];
    538   if (naxis == 3)  nelements*=naxes[2];
    539   int anynull;
    540   double dnullval=0.;
    541   float fnullval=0.;
    542   // on laisse a fits le soin de convertir le type du tableau lu vers
    543   // le type suppose par l'utilisateur de fitsioserver
    544   //
    545   switch ( FITS_tab_typ_)
    546     {
    547     case TDOUBLE :
    548       if (bitpix !=  DOUBLE_IMG)
    549         {
    550           cout << " FitsIOServer : the data type on fits file is not double, "
    551                << " conversion to double will be achieved by cfitsio lib " << endl;
    552         }
    553       if (dtab_ != NULL) delete [] dtab_;
    554       dtab_=new double[nelements];
    555       fits_read_img(fptr, TDOUBLE, 1, nelements,  &dnullval, dtab_,
    556                     &anynull, &status);
    557       if( status )  printerror( status );
    558       break;
    559     case TFLOAT :
    560       if (bitpix !=  FLOAT_IMG)
    561         {
    562           cout << " FitsIOServer : the data type on fits file is not float, "
    563                << " conversion to float will be achieved by cfitsio lib " << endl;
    564         }
    565       if (ftab_ != NULL) delete [] ftab_;
    566       ftab_=new float[nelements];
    567       fits_read_img(fptr, TFLOAT, 1, nelements,  &fnullval, ftab_,
    568                     &anynull, &status);
    569       if( status )  printerror( status );
    570       break;
    571 
    572 
    573     case TINT :
    574       if (bitpix !=  LONG_IMG)
    575         {
    576           cout << " FitsIOServer : the data type on fits file is not long, "
    577                << " conversion to long will be achieved by cfitsio lib " << endl;
    578         }
    579       if (i_4tab_ != NULL) delete [] i_4tab_;
    580       i_4tab_=new int_4[nelements];
    581       fits_read_img(fptr, TINT, 1, nelements,  &fnullval, i_4tab_,
    582                     &anynull, &status);
    583       if( status )  printerror( status );
    584       break;
    585 
    586 
    587     default :
    588       cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
    589       break;
    590     }
    591   status = 0;
    592   char card[FLEN_CARD];
    593   int num = 0;
    594   char comment2[FLEN_COMMENT] = "x";
    595   char keyname[]= "";
    596   char datekey[]= "DATE";
    597   char endkey[] = "END";
    598   char typ='x';
    599   int ival;
    600   double dval;
    601   char strval[]="";
    602   // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
    603   // on va jusqu'au mot cle DATE
    604   int flen_keyword = 9;
    605   // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
    606   // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
    607   // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
    608   while (status == 0 &&  strncmp(keyname, datekey,4) != 0 )
    609     {
    610       num++;
    611       fits_read_record(fptr, num, card, &status);
    612       strncpy(keyname,card,flen_keyword-1);
    613     }
    614   if (status != 0 )
    615     {
    616       cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
    617     }
    618   // on recupere la liste des mots-cles utilisateurs
    619   while (status == 0)
    620     {
    621       num++;
    622       // on lit un record pour recuperer le nom du mot-cle
    623       fits_read_record(fptr, num, card, &status);
    624       strncpy(keyname,card,flen_keyword-1);
    625       char value[FLEN_VALUE];
    626       // on recupere le premier caractere du commentaire, qui contient
    627       // le code du type de la valeur
    628       // (tant que ce n est pas le mot cle END)
    629       fits_read_keyword(fptr, keyname, value, comment2, &status);
    630       if (  strncmp(keyname, endkey,flen_keyword-1) != 0)
    631         {
    632           typ = comment2[0];
    633           // quand le type est connu, on lit la valeur
    634           switch (typ)
    635             {
    636             case 'I' :
    637               fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
    638               if( status )  printerror( status );
    639               strip (keyname, 'B',' ');
    640               dvl[keyname] = ival;
    641               break;
    642             case 'D' :
    643               fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
    644               if( status )  printerror( status );
    645               strip (keyname, 'B',' ');
    646               dvl[keyname] = dval;
    647               break;
    648             case 'S' :
    649               fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
    650               if( status )  printerror( status );
    651               strip (keyname, 'B',' ');
    652               strip(strval, 'B',' ');
    653               dvl[keyname]=strval;
    654               break;
    655             default :
    656               cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
    657               break;
    658             }
    659         }
    660     }
    661 
    662 }
    663 
    664 
    665 void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
    666 {
    667    
    668   long naxes[2]={600, 300};
    669   float map [ 600*300 ];
    670   int npixels=  naxes[0]*naxes[1];
    671 
    672    cout << " image FITS en projection SINUS" << endl;
    673   // table will have npixels rows
    674   for(int j=0; j < npixels; j++) map[j]=0.;
    675     for(int j=0; j<naxes[1]; j++)
    676       {
    677         double yd = (j+0.5)/naxes[1]-0.5;
    678         double theta = (0.5-yd)*Pi;
    679         double facteur=1./sin(theta);
    680         for(int i=0; i<naxes[0]; i++) 
    681           {
    682             double xa = (i+0.5)/naxes[0]-0.5;
    683             double phi = 2.*Pi*xa*facteur+Pi;
    684             float th=float(theta);
    685             float ff=float(phi);
    686             if (phi<2*Pi && phi>= 0)
    687               {
    688                 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
    689               }
    690           }
    691       }
    692  
    693     write_picture(naxes, map,  filename);
    694 
    695 }
    696 void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
    697 {
    698   // le code de cete methode duplique celui de la precedente, la seule
    699   //difference etant le type de sphere en entree. Ces methodes de projection
    700   // sont provisoires, et ne servent que pour les tests. C est pourquoi je
    701   // ne me suis pas casse la tete, pour l instant
    702    
    703   long naxes[2]={600, 300};
    704   float map [ 600*300 ];
    705   int npixels=  naxes[0]*naxes[1];
    706 
    707    cout << " image FITS en projection SINUS" << endl;
    708   // table will have npixels rows
    709   for(int j=0; j < npixels; j++) map[j]=0.;
    710     for(int j=0; j<naxes[1]; j++)
    711       {
    712         double yd = (j+0.5)/naxes[1]-0.5;
    713         double theta = (0.5-yd)*Pi;
    714         double facteur=1./sin(theta);
    715         for(int i=0; i<naxes[0]; i++) 
    716           {
    717             double xa = (i+0.5)/naxes[0]-0.5;
    718             double phi = 2.*Pi*xa*facteur+Pi;
    719             float th=float(theta);
    720             float ff=float(phi);
    721             if (phi<2*Pi && phi>= 0)
    722               {
    723                 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
    724               }
    725           }
    726       }
    727  
    728     write_picture(naxes, map,  filename);
    729 
    730 
    731 }
    732 
    733 void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
    734 {
    735    
    736   long naxes[2];
    737   naxes[0] = lcm.Size_x();
    738   naxes[1] = lcm.Size_y();
    739   int npixels=  naxes[0]*naxes[1];
    740   float* map = new float[npixels];
    741 
    742   // table will have npixels rows
    743   for(int j=0; j < npixels; j++) map[j]=0.;
    744     for(int j=0; j<naxes[1]; j++)
    745       {
    746         for(int i=0; i<naxes[0]; i++) 
    747           {
    748                 map[j*naxes[0]+i] = lcm(i, j);
    749           }
    750       }
    751  
    752     write_picture(naxes, map,  filename);
    753     delete [] map;
    754 }
    755 
    756 void  FitsIoServer::read_sphere(char flnm[], int& npixels, int& nside)
    757 {
    758   int status=0;
    759   long naxis;
    760   nside=0;
    761   int n1, n2, n3;
    762   DVList dvl;
    763   // pointer to the FITS file, defined in fitsio.h
    764   fitsfile *fptr;
    765   fits_open_file(&fptr, flnm, READONLY, &status);
    766   if( status )  printerror( status );
    767   planck_read_img(fptr, naxis, n1, n2, n3, dvl);
    768   // close the file
    769   fits_close_file(fptr, &status);
    770   if( status )  printerror( status );
    771 
    772   npixels=n1;
    773   //dvl.Print();
    774   DVList::ValList::const_iterator it;
    775   for(it = dvl.Begin(); it != dvl.End(); it++) 
    776     {
    777       char datatype= (*it).second.typ;
    778       char keyname[]="";
    779       strcpy(keyname, (*it).first.substr(0,64).c_str());   
    780       int ival=0;
    781       double dval=0.;
    782       char strval[]="";
    783       switch (datatype)
    784         {
    785         case 'I' :
    786           ival=(*it).second.mtv.iv;
    787           break;
    788         case 'D' :
    789           dval=(*it).second.mtv.dv;
    790           break;
    791         case 'S' :
    792           strcpy(strval, (*it).second.mtv.strv);
    793           break;
    794         default :
    795           cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    796           break;
    797         }
    798       if (!strcmp(keyname, "NSIDE") )
    799         {
    800           nside = ival;
    801         }
    802       //if (!strcmp(keyname, "ORDERING") )
    803       //        {
    804       //          cout << " j'ai trouve ORDERING " << endl;
    805       //        }
    806     }
    807 
    808   if (naxis != 1)
    809     {
    810       cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
    811     }
    812 }
    813 
    814 
    815 void  FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
    816 {
    817   //pointer to the FITS file, defined in fitsio.h
    818   fitsfile *fptr; 
    819   int bitpix = FLOAT_IMG;
    820   long naxis = 2;
    821 
    822   // delete old file if it already exists
    823   remove(filename);
    824  
    825   // initialize status before calling fitsio routines
    826   int status = 0;         
    827 
    828   // create new FITS file 
    829   fits_create_file(&fptr, filename, &status);
    830   if( status )  printerror( status );
    831 
    832   // write the required header keywords
    833   fits_create_img(fptr, bitpix, naxis, naxes, &status);
    834   if( status )  printerror( status );
    835 
    836   // write the current date
    837   fits_write_date(fptr, &status);
    838   if( status )  printerror( status );
    839 
    840 
    841   // first row in table to write
    842   long firstrow  = 1;
    843   // first element in row   
    844   long firstelem = 1; 
    845   int colnum = 1;
    846   int nelements=naxes[0]*naxes[1];
    847   fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
    848   if( status )  printerror( status );
    849 
    850   // close the file
    851   fits_close_file(fptr, &status);
    852   if( status )  printerror( status );
    853 }
    854 
    855 void FitsIoServer::save(NTuple& ntpl,char flnm[])
    856 
    857   //****************************************************/
    858   //* read the elements of the NTuple ntpl, and create */
    859   //* a FITS file with a binary table extension        */
    860   //****************************************************/
    861 {
    862 
    863   // delete old file if it already exists
    864   remove(flnm);
    865 
    866   // pointer to the FITS file, defined in fitsio.h
    867   fitsfile *fptr;
    868   int status = 0;
    869   if( fits_create_file(&fptr, flnm, &status) )
    870     printerror( status );
    871 
    872   // table will have tfields columns
    873   int tfields= ntpl.NVar();
    874 
    875   // table will have nrows rows
    876   int nrows= ntpl.NEntry();
    877 
    878   // extension name
    879   char extname[] = "NTuple_Binary_tbl";
    880 
    881   // define the name, and the datatype for the tfields columns
    882   char **ttype, **tform;
    883   ttype= new char*[tfields];
    884   tform= new char*[tfields];
    885   for(int i = 0; i < tfields; i++)
    886     {
    887       ttype[i]= new char[FLEN_VALUE];
    888       strcpy(ttype[i], ntpl.NomIndex(i));
    889       tform[i]= new char[FLEN_VALUE];
    890       strcpy(tform[i], "1E");
    891     }
    892 
    893   // create a new empty binary table onto the FITS file
    894   // physical units if they exist, are defined in the DVList object
    895   // so the null pointer is given for the tunit parameters.
    896   if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
    897                        NULL,extname,&status) )
    898     printerror( status );
    899 
    900   // first row in table to write
    901   int firstrow  = 1;
    902 
    903   // first element in row  (ignored in ASCII tables)
    904   int firstelem = 1;
    905 
    906   for(int i = 0; i < tfields; i++)
    907     {
    908       float *dens= new float[nrows];
    909       for(int j = 0; j < nrows; j++)
    910         {
    911           dens[j]= ntpl.GetVal(j,i);
    912         }
    913       fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
    914       delete []dens;
    915     }
    916 
    917   // number of blocks per event
    918   int blk= ntpl.BLock();
    919   fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
    920 
    921   // get names and values from the join DVList object
    922   DVList dvl= ntpl.Info();
    923   dvl.Print();
    924   DVList::ValList::const_iterator it;
    925   for(it = dvl.Begin(); it != dvl.End(); it++)
    926     {
    927       char keytype= (*it).second.typ;
    928       char keyname[9]="";
    929       strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
    930       char comment[FLEN_COMMENT]="";
    931 
    932       switch (keytype)
    933         {
    934         case 'I' :
    935           {
    936             int ival=(*it).second.mtv.iv;
    937             strcpy(comment,"I entier");
    938             fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
    939             break;
    940           }
    941         case 'D' :
    942           {
    943             double dval=(*it).second.mtv.dv;
    944             strcpy(comment,"D double");
    945             fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
    946             break;
    947           }
    948         case 'S' :
    949           {
    950             char strval[]="";
    951             strcpy(strval,(*it).second.mtv.strv);
    952             strcpy(comment,"S character string");
    953             fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
    954             break;
    955           }
    956         }
    957     }
    958 
    959   //close the FITS file
    960   if ( fits_close_file(fptr, &status) )
    961     printerror( status );
    962   return;
    963 } 
    96476
    96577void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum)
    … …  
    1161273}
    1162274
    1163 bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
    1164 
    1165   //*****************************************************/
    1166   //* check if the specified keyword exits in the CHU   */
    1167   //*****************************************************/
    1168 {
    1169 
    1170   bool KEY_EXIST = false;
     275
     276void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
     277{
     278  int npixels=0;
     279  int nside=0;
     280  FITS_tab_typ_ = TDOUBLE;
     281  read_sphere(flnm, npixels, nside);
     282
     283  // number of pixels in the sphere
     284  if (sph.NbPixels() != npixels)
     285    {
     286      cout << " found " << npixels << " pixels" << endl;
     287      cout << " expected " << sph.NbPixels() << endl;
     288      if (nside <= 0 )
     289        {
     290          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
     291          exit(0);
     292        }
     293      if (nside != sph.SizeIndex())
     294        {
     295          sph.Resize(nside);
     296          cout << " resizing the sphere to nside=  " << nside  << endl;
     297        }
     298      else
     299        {
     300          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
     301          exit(0);
     302        }
     303    }
     304  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
     305}
     306void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
     307{
     308  int npixels=0;
     309  int nside=0;
     310  FITS_tab_typ_ = TFLOAT;
     311  read_sphere(flnm, npixels, nside);
     312  // number of pixels in the sphere
     313  if (sph.NbPixels() != npixels)
     314    {
     315      cout << " found " << npixels << " pixels" << endl;
     316      cout << " expected " << sph.NbPixels() << endl;
     317      if (nside <= 0 )
     318        {
     319          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
     320          exit(0);
     321        }
     322      if (nside != sph.SizeIndex())
     323        {
     324          sph.Resize(nside);
     325          cout << " resizing the sphere to nside=  " << nside  << endl;
     326        }
     327      else
     328        {
     329          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
     330          exit(0);
     331        }
     332    }
     333  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
     334}
     335
     336void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
     337{
     338  int nside;
     339  int nbrows=0;
     340  int nbcols=0;
     341  FITS_tab_typ_ = TDOUBLE;
     342  long naxis;
     343  int n1, n2, n3;
     344  DVList dvl;
     345  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
     346
     347  nbrows=n1;
     348  nbcols=n2;
     349  if (naxis != 2)
     350    {
     351      cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
     352    }
     353  //dvl.Print();
     354  DVList::ValList::const_iterator it;
     355  for(it = dvl.Begin(); it != dvl.End(); it++) 
     356    {
     357      char datatype= (*it).second.typ;
     358      char keyname[]="";
     359      strcpy(keyname, (*it).first.substr(0,64).c_str());   
     360      int ival=0;
     361      double dval=0.;
     362      char strval[]="";
     363      switch (datatype)
     364        {
     365        case 'I' :
     366          ival=(*it).second.mtv.iv;
     367          break;
     368        case 'D' :
     369          dval=(*it).second.mtv.dv;
     370          break;
     371        case 'S' :
     372          strcpy(strval, (*it).second.mtv.strv);
     373          break;
     374        default :
     375          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
     376          break;
     377        }
     378      if (!strcmp(keyname, "NSIDE") )
     379        {
     380          nside = ival;
     381        }
     382      //if (!strcmp(keyname, "ORDERING") )
     383      //        {
     384      //  cout << " j'ai trouve ORDERING " << endl;
     385      //        }
     386    }
     387  float theta0 = dvl.GetD("THETA0");
     388  float phi0 = dvl.GetD("PHI0");
     389  int x0 = dvl.GetI("X0");
     390  int y0 = dvl.GetI("Y0");
     391  int xo = dvl.GetI("XO");
     392  int yo = dvl.GetI("YO");
     393  float anglex=dvl.GetD("ANGLEX");
     394  float angley=dvl.GetD("ANGLEY");
     395  float angle = dvl.GetD("ANGLE");
     396
     397  // number of components
     398
     399  if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
     400    {
     401      cout << " found " << nbrows << " x-pixels ";
     402      cout << "  expected " << lcm.Size_x()  << endl;
     403      cout << " found " << nbcols << " y-pixels " ;
     404      cout << "  expected " <<   lcm.Size_y()  << endl;
     405      lcm.ReSize(nbrows,nbcols);
     406      cout << " resizing the map to x-size=  " << nbrows  << " y-size= " << nbcols << endl;
     407    }
     408  lcm.SetSize(anglex, angley);
     409  lcm.SetOrigin(theta0, phi0, x0, y0, angle);
     410  int ij=0;
     411  for (int j=0; j< nbcols; j++)
     412    for (int i = 0; i < nbrows; i++)  lcm(i,j) = dtab_[ij++];
     413}
     414
     415void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
     416{
     417  FITS_tab_typ_ = TFLOAT;
     418  long naxis;
     419  int siz_x;
     420  int siz_y;
     421  int n3;
     422  DVList dvl;
     423
     424  planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
     425
     426 
     427  if (naxis != 2)
     428    {
     429      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
     430    }
     431
     432  DpcImg.Allocate(siz_x, siz_y, 0., 0);
     433  float* pixelPtr=DpcImg.ImagePtr();
     434  //for (int i=0; i < siz_x*siz_y; i++) pixelPtr[i]=ftab_[i];
     435  // verifications de type
     436  PBaseDataTypes  dataT=DataType((r_4)0);
     437  int TypeDonnees=dvl.GetI("DATATYPE");
     438  if (int(dataT) !=  TypeDonnees)
     439    {
     440      cout << " FitsIOServer : parameter DATATYPE on file is not float " << endl;
     441      cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
     442    }
     443
     444  memcpy(pixelPtr, ftab_, siz_x*siz_y*DataSize(dataT));
     445  const char* nom=dvl.GetS("NAME").c_str();
     446  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
     447  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
     448  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
     449  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
     450       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
     451       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
     452
     453}
     454
     455
     456void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
     457{
     458  FITS_tab_typ_ = TINT;
     459  long naxis;
     460  int siz_x;
     461  int siz_y;
     462  int n3;
     463  DVList dvl;
     464
     465  // pointer to the FITS file, defined in fitsio.h
     466  //fitsfile *fptr;
     467  // initialize status before calling fitsio routines
     468  //int status = 0;
     469       
     470  //fits_open_file(&fptr, flnm, READONLY, &status);
     471  //if( status )  printerror( status );
     472  planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
     473  // close the file
     474  //fits_close_file(fptr, &status);
     475  //if( status )  printerror( status );
     476 
     477  if (naxis != 2)
     478    {
     479      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
     480    }
     481
     482  DpcImg.Allocate(siz_x, siz_y, 0., 0);
     483  int_4* pixelPtr=DpcImg.ImagePtr();
     484
     485
     486  // verifications de type
     487  PBaseDataTypes  dataT=DataType((int_4)0);
     488  int TypeDonnees=dvl.GetI("DATATYPE");
     489  if (int(dataT) !=  TypeDonnees)
     490    {
     491      cout << " FitsIOServer : parameter DATATYPE on file is not int_4 " << endl;
     492      cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
     493    }
     494  memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
     495  const char* nom=dvl.GetS("NAME").c_str();
     496  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
     497  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
     498  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
     499  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
     500       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
     501       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
     502}
     503
     504
     505void FitsIoServer::save( TMatrix<double>& mat, char filename[])   
     506{
     507  int nbrows = mat.NRows();
     508  int nbcols = mat.NCols();
     509  long naxis = nbcols > 1 ? 2 : 1;
     510  cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
     511  FITS_tab_typ_ = TDOUBLE;
     512  if (dtab_ != NULL ) delete[] dtab_;
     513  dtab_=new double[nbrows*nbcols];
     514
     515
     516  int ij=0;
     517  for (int j=0; j< nbcols; j++)
     518    for (int i = 0; i < nbrows; i++) dtab_[ij++]= mat(i,j);
     519
     520  DVList dvl;
     521
     522  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
     523
     524}
     525
     526void FitsIoServer::save(NTuple& ntpl,char flnm[])
     527
     528  //****************************************************/
     529  //* read the elements of the NTuple ntpl, and create */
     530  //* a FITS file with a binary table extension        */
     531  //****************************************************/
     532{
     533
     534  // delete old file if it already exists
     535  remove(flnm);
     536
     537  // pointer to the FITS file, defined in fitsio.h
     538  fitsfile *fptr;
    1171539  int status = 0;
    1172   char strbide[FLEN_VALUE];
    1173   char keybide[LEN_KEYWORD]= "";
    1174   for(int jj = 1; jj <= nkeys; jj++)
    1175     {
    1176       if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
    1177         printerror( status );
    1178       if( !strcmp(keybide,keyword) )
    1179         {
    1180           KEY_EXIST= true;
    1181           break;
    1182         }
    1183     }
    1184   return(KEY_EXIST);
    1185 }
    1186 
    1187 void FitsIoServer::readheader ( char filename[] )
    1188 
    1189   //**********************************************************************/
    1190   //* Print out all the header keywords in all extensions of a FITS file */
    1191   //**********************************************************************/
    1192 {
    1193 
    1194   // standard string lengths defined in fitsioc.h
    1195   char card[FLEN_CARD];
    1196 
    1197   // pointer to the FITS file, defined in fitsio.h
    1198   fitsfile *fptr;
    1199 
    1200   int status = 0;
    1201   if ( fits_open_file(&fptr, filename, READONLY, &status) )
     540  if( fits_create_file(&fptr, flnm, &status) )
    1202541    printerror( status );
    1203  
    1204   // attempt to move to next HDU, until we get an EOF error
    1205   int hdutype;
    1206   for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
    1207     {
    1208       if (hdutype == ASCII_TBL)
    1209         printf("\nReading ASCII table in HDU %d:\n", ii);
    1210       else if (hdutype == BINARY_TBL)
    1211         printf("\nReading binary table in HDU %d:\n", ii);
    1212       else if (hdutype == IMAGE_HDU)
    1213         printf("\nReading FITS image in HDU %d:\n", ii);
    1214       else
    1215         {
    1216           printf("Error: unknown type of this HDU \n");
    1217           printerror( status );
    1218         }
    1219 
    1220       // get the number of keywords
    1221       int nkeys, keypos;
    1222       if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
    1223         printerror( status );
    1224      
    1225       printf("Header listing for HDU #%d:\n", ii);
    1226       for (int jj = 1; jj <= nkeys; jj++) 
    1227         {
    1228           if ( fits_read_record(fptr, jj, card, &status) )
    1229             printerror( status );
    1230 
    1231           // print the keyword card
    1232           printf("%s\n", card);
    1233         }
    1234       printf("END\n\n");
    1235     }
    1236 
    1237   // got the expected EOF error; reset = 0
    1238   if (status == END_OF_FILE)   
    1239     status = 0;             
    1240   else
     542
     543  // table will have tfields columns
     544  int tfields= ntpl.NVar();
     545
     546  // table will have nrows rows
     547  int nrows= ntpl.NEntry();
     548
     549  // extension name
     550  char extname[] = "NTuple_Binary_tbl";
     551
     552  // define the name, and the datatype for the tfields columns
     553  char **ttype, **tform;
     554  ttype= new char*[tfields];
     555  tform= new char*[tfields];
     556  for(int i = 0; i < tfields; i++)
     557    {
     558      ttype[i]= new char[FLEN_VALUE];
     559      strcpy(ttype[i], ntpl.NomIndex(i));
     560      tform[i]= new char[FLEN_VALUE];
     561      strcpy(tform[i], "1E");
     562    }
     563
     564  // create a new empty binary table onto the FITS file
     565  // physical units if they exist, are defined in the DVList object
     566  // so the null pointer is given for the tunit parameters.
     567  if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
     568                       NULL,extname,&status) )
    1241569    printerror( status );
    1242570
     571  // first row in table to write
     572  int firstrow  = 1;
     573
     574  // first element in row  (ignored in ASCII tables)
     575  int firstelem = 1;
     576
     577  for(int i = 0; i < tfields; i++)
     578    {
     579      float *dens= new float[nrows];
     580      for(int j = 0; j < nrows; j++)
     581        {
     582          dens[j]= ntpl.GetVal(j,i);
     583        }
     584      fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
     585      delete []dens;
     586    }
     587
     588  // number of blocks per event
     589  int blk= ntpl.BLock();
     590  fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
     591
     592  // get names and values from the join DVList object
     593  DVList dvl= ntpl.Info();
     594  dvl.Print();
     595  DVList::ValList::const_iterator it;
     596  for(it = dvl.Begin(); it != dvl.End(); it++)
     597    {
     598      char keytype= (*it).second.typ;
     599      char keyname[9]="";
     600      strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
     601      char comment[FLEN_COMMENT]="";
     602
     603      switch (keytype)
     604        {
     605        case 'I' :
     606          {
     607            int ival=(*it).second.mtv.iv;
     608            strcpy(comment,"I entier");
     609            fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
     610            break;
     611          }
     612        case 'D' :
     613          {
     614            double dval=(*it).second.mtv.dv;
     615            strcpy(comment,"D double");
     616            fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
     617            break;
     618          }
     619        case 'S' :
     620          {
     621            char strval[]="";
     622            strcpy(strval,(*it).second.mtv.strv);
     623            strcpy(comment,"S character string");
     624            fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
     625            break;
     626          }
     627        }
     628    }
     629
     630  //close the FITS file
    1243631  if ( fits_close_file(fptr, &status) )
    1244632    printerror( status );
     633  return;
     634} 
     635
     636
     637
     638
     639//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
     640void FitsIoServer::save(SphericalMap<double>& sph, char filename[])   
     641{
     642  int npixels = sph.NbPixels();
     643  FITS_tab_typ_ = TDOUBLE;
     644  if (dtab_ != NULL ) delete[] dtab_;
     645  dtab_=new double[npixels];
     646
     647
     648  for (int j = 0; j < npixels; j++) dtab_[j]= sph(j);
     649  DVList dvl;
     650  dvl["NSIDE"] = sph.SizeIndex();
     651  dvl["ORDERING"]=sph.TypeOfMap();
     652
     653  planck_write_img(filename, 1, npixels, 0, 0, dvl);
     654
     655  // decider ulterieurement ce qu'on fait de ce qui suit, specifique
     656  // pour l'instant, aux spheres gorski.
     657
     658  /*
     659  // write supplementary keywords
     660  fits_write_comment(fptr, " ", &status);
     661  if( status ) printerror( status );
     662
     663
     664  strcpy(comment,"HEALPIX Pixelisation");
     665  strcpy(svalue, "HEALPIX");
     666  fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
     667  if( status )  printerror( status );
     668
     669
     670  strcpy(comment,"pixel ordering scheme, either RING or NESTED");
     671  strcpy(svalue, "RING");
     672  fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
     673  if( status )  printerror( status );
     674
     675
     676  strcpy(comment,"Random generator seed");
     677  int iseed= sph.iseed();
     678  fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
     679  if( status )  printerror( status );
     680
     681
     682  strcpy(comment,"--------------------------------------------");
     683  fits_write_comment(fptr, comment, &status);
     684  if( status ) printerror( status );
     685
     686  strcpy(comment," Above keywords are still likely to change !");
     687  fits_write_comment(fptr, comment, &status);
     688  if( status ) printerror( status );
     689
     690  strcpy(comment,"--------------------------------------------");
     691  fits_write_comment(fptr, comment, &status);
     692  if( status ) printerror( status );
    1245693 
    1246   return;
    1247 }
    1248 
    1249 void FitsIoServer::printerror(int status) const
    1250 
    1251   //*****************************************************/
    1252   //* Print out cfitsio error messages and exit program */
    1253   //*****************************************************/
    1254 {
    1255 
    1256   // print out cfitsio error messages and exit program
    1257   if( status )
    1258     {
    1259       // print error report 
    1260       fits_report_error(stderr, status);
    1261       // terminate the program, returning error status
    1262       exit( status );
    1263     }
    1264   return;
    1265 }
     694  */
     695
     696
     697}
     698
     699void FitsIoServer::save(SphericalMap<float>& sph, char filename[])   
     700{
     701  int npixels = sph.NbPixels();
     702  FITS_tab_typ_ = TFLOAT;
     703  if (ftab_ != NULL ) delete[] ftab_;
     704  ftab_=new float[npixels];
     705
     706
     707  for (int j = 0; j < npixels; j++) ftab_[j]= sph(j);
     708  DVList dvl;
     709  dvl["NSIDE"] = sph.SizeIndex();
     710  dvl["ORDERING"]=sph.TypeOfMap();
     711
     712  planck_write_img(filename, 1, npixels, 0, 0, dvl);
     713
     714
     715}
     716
     717void FitsIoServer::save(LocalMap<double>& locm, char filename[])   
     718{
     719  int nbrows = locm.Size_x();
     720  int nbcols = locm.Size_y();
     721  long naxis = 2;
     722  cout << " nombre de pts en x : " << nbrows << " en y  " << nbcols << endl;
     723  FITS_tab_typ_ = TDOUBLE;
     724  if (dtab_ != NULL ) delete[] dtab_;
     725  dtab_=new double[nbrows*nbcols];
     726
     727  int ij=0;
     728  for (int j=0; j< nbcols; j++)
     729    for (int i = 0; i < nbrows; i++) dtab_[ij++]= locm(i,j);
     730
     731  DVList dvl;
     732  dvl["NSIDE"] = locm.SizeIndex();
     733  dvl["ORDERING"]=locm.TypeOfMap();
     734  double theta0;
     735  double phi0;
     736  int x0;
     737  int y0;
     738  double angle;
     739  locm.Origin(theta0,phi0,x0,y0,angle);
     740  double anglex;
     741  double angley;
     742  locm.Aperture(anglex,angley);
     743  dvl["THETA0"] = theta0;
     744  dvl["PHI0"] = phi0;
     745  dvl["X0"] = (int_4)x0;
     746  dvl["Y0"] = (int_4)y0;
     747  dvl["ANGLE"] = angle;
     748  dvl["ANGLEX"] = anglex;
     749  dvl["ANGLEY"] = angley;
     750  planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
     751}
     752
    1266753
    1267754void FitsIoServer::save(ImageR4& DpcImg,char flnm[])
    … …  
    1273760  FITS_tab_typ_ = TFLOAT;
    1274761
    1275   // pointer to the FITS file, defined in fitsio.h
    1276   fitsfile *fptr;
    1277 
    1278   // initialize status before calling fitsio routines
    1279   int status = 0;         
    1280 
    1281   // delete old file if it already exists
    1282   remove(flnm);
    1283 
    1284   // create new FITS file 
    1285   fits_create_file(&fptr, flnm, &status);
    1286   if( status )  printerror( status );
    1287762
    1288763  // get the values of the DpcImage
    … …  
    1310785
    1311786
    1312   planck_write_img(fptr, naxis, siz_x, siz_y, 0, dvl);
    1313 
    1314   // close the file
    1315   fits_close_file(fptr, &status);
    1316   if( status )  printerror( status );
     787  planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
     788
    1317789
    1318790  ftab_=NULL;
    1319791}
    1320792
    1321 void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
    1322 {
    1323   FITS_tab_typ_ = TFLOAT;
    1324   long naxis;
    1325   int siz_x;
    1326   int siz_y;
    1327   int n3;
    1328   DVList dvl;
    1329 
    1330   // pointer to the FITS file, defined in fitsio.h
    1331   fitsfile *fptr;
    1332   // initialize status before calling fitsio routines
    1333   int status = 0;
    1334        
    1335   fits_open_file(&fptr, flnm, READONLY, &status);
    1336   if( status )  printerror( status );
    1337 
    1338   planck_read_img(fptr, naxis, siz_x, siz_y, n3, dvl);
    1339 
    1340   // close the file
    1341   fits_close_file(fptr, &status);
    1342   if( status )  printerror( status );
    1343  
    1344   if (naxis != 2)
    1345     {
    1346       cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
    1347     }
    1348 
    1349   DpcImg.Allocate(siz_x, siz_y, 0., 0);
    1350   float* pixelPtr=DpcImg.ImagePtr();
    1351   //for (int i=0; i < siz_x*siz_y; i++) pixelPtr[i]=ftab_[i];
    1352   // verifications de type
    1353   PBaseDataTypes  dataT=DataType((r_4)0);
    1354   int TypeDonnees=dvl.GetI("DATATYPE");
    1355   if (int(dataT) !=  TypeDonnees)
    1356     {
    1357       cout << " FitsIOServer : parameter DATATYPE on file is not float " << endl;
    1358       cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
    1359     }
    1360 
    1361   memcpy(pixelPtr, ftab_, siz_x*siz_y*DataSize(dataT));
    1362   const char* nom=dvl.GetS("NAME").c_str();
    1363   DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
    1364   DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
    1365   DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
    1366   DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
    1367        dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
    1368        dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
    1369 
    1370 }
    1371793
    1372794void FitsIoServer::save(ImageI4& DpcImg,char flnm[])
    … …  
    1379801
    1380802  // pointer to the FITS file, defined in fitsio.h
    1381   fitsfile *fptr;
     803  //fitsfile *fptr;
    1382804
    1383805  // initialize status before calling fitsio routines
    1384   int status = 0;         
     806  //int status = 0;         
    1385807
    1386808  // delete old file if it already exists
    1387   remove(flnm);
     809  //remove(flnm);
    1388810
    1389811  // create new FITS file 
    1390   fits_create_file(&fptr, flnm, &status);
    1391   if( status )  printerror( status );
     812  //fits_create_file(&fptr, flnm, &status);
     813  //if( status )  printerror( status );
    1392814
    1393815  // get the values of the DpcImage
    … …  
    1415837
    1416838
    1417   planck_write_img(fptr, naxis, siz_x, siz_y, 0, dvl);
    1418 
     839  planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
     840
     841  // close the file
     842  //fits_close_file(fptr, &status);
     843  //if( status )  printerror( status );
     844
     845  i_4tab_=NULL;
     846}
     847
     848
     849
     850
     851void FitsIoServer::planck_write_img(char flnm[], int naxis,int n1, int n2, int n3, DVList& dvl)
     852{
     853
     854  cout << " nouveau write_img" << endl;
     855
     856  // pointer to the FITS file, defined in fitsio.h
     857  fitsfile *fptr;
     858
     859  // initialize status before calling fitsio routines
     860  int status = 0;         
     861
     862  // delete old file if it already exists
     863  remove(flnm);
     864
     865  // create new FITS file 
     866  fits_create_file(&fptr, flnm, &status);
     867  if( status )  printerror( status );
     868
     869  int bitpix=0;
     870  if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
     871  if ( FITS_tab_typ_ == TFLOAT)  bitpix = FLOAT_IMG;
     872  if ( FITS_tab_typ_ == TINT)   bitpix = LONG_IMG;
     873  long naxes[3];
     874  naxes[0] = n1;
     875  naxes[1] = n2;
     876  naxes[2] = n3;
     877  if (n2 > 0 && naxis < 2) cout << "  FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
     878  if (n3 > 0 && naxis < 3) cout << "  FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
     879  fits_create_img(fptr, bitpix, naxis, naxes, &status);
     880  if( status )  printerror( status );
     881
     882
     883
     884  long nelements= naxes[0];
     885  if (naxis >=2) nelements*=naxes[1];
     886  if (naxis == 3)  nelements*=naxes[2];
     887  switch (FITS_tab_typ_)
     888    {
     889    case TDOUBLE :
     890      fits_write_img(fptr, TDOUBLE, 1, nelements, dtab_,  &status);
     891      if( status ) printerror( status );
     892      break;
     893    case TFLOAT :
     894      fits_write_img(fptr, TFLOAT, 1, nelements, ftab_,  &status);
     895      if( status ) printerror( status );
     896      break;
     897    case TINT :
     898      fits_write_img(fptr, TINT, 1, nelements, i_4tab_,  &status);
     899      if( status ) printerror( status );
     900      break;
     901    default :
     902      cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
     903      break;
     904    }
     905
     906  // write the current date
     907  fits_write_date(fptr, &status);
     908  if( status )  printerror( status );
     909  // on convient d ecrire apres la date, les mots cles utilisateur.
     910  // la date servira de repere pour relire ces mots cles.
     911
     912  //dvl.Print();
     913  DVList::ValList::const_iterator it;
     914  for(it = dvl.Begin(); it != dvl.End(); it++) 
     915    {
     916      int datatype=  key_type_PL2FITS( (*it).second.typ);
     917      char keyname[]="";
     918      int flen_keyword = 9;
     919      // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
     920      // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
     921      // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
     922      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
     923      char comment[FLEN_COMMENT]="";
     924      int ival=0;
     925      double dval=0.;
     926      char strval[]="";
     927      switch (datatype)
     928        {
     929        case TINT :
     930          ival=(*it).second.mtv.iv;
     931          strcpy(comment,"I entier");
     932          fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
     933          break;
     934        case TDOUBLE :
     935          dval=(*it).second.mtv.dv;
     936          strcpy(comment,"D double");
     937          fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
     938          break;
     939        case TSTRING :
     940          strcpy(strval, (*it).second.mtv.strv);
     941          strcpy(comment,"S character string");
     942          fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
     943          break;
     944        default :
     945          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
     946          break;
     947        }
     948      if( status ) printerror( status );
     949    }
    1419950  // close the file
    1420951  fits_close_file(fptr, &status);
    1421952  if( status )  printerror( status );
    1422 
    1423   i_4tab_=NULL;
    1424 }
    1425 
    1426 void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
    1427 {
    1428   FITS_tab_typ_ = TINT;
    1429   long naxis;
    1430   int siz_x;
    1431   int siz_y;
    1432   int n3;
    1433   DVList dvl;
     953}
     954
     955void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
     956{
     957  int status=0;
     958  long bitpix;
     959  long naxes[3]={0,0,0};
     960  char* comment=NULL;
    1434961
    1435962  // pointer to the FITS file, defined in fitsio.h
    1436963  fitsfile *fptr;
    1437964  // initialize status before calling fitsio routines
    1438   int status = 0;
    1439        
    1440965  fits_open_file(&fptr, flnm, READONLY, &status);
    1441966  if( status )  printerror( status );
    1442   planck_read_img(fptr, naxis, siz_x, siz_y, n3, dvl);
     967
     968
     969  fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
     970  if( status )  printerror( status );
     971  fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
     972  if( status )  printerror( status );
     973  int nfound;
     974  int nkeys=(int)naxis;
     975  fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
     976  if( status )  printerror( status );
     977
     978  n1 = naxes[0] ;
     979  n2 = naxes[1] ;
     980  n3 = naxes[2] ;
     981
     982
     983  long nelements= naxes[0];
     984  if (naxis >=2) nelements*=naxes[1];
     985  if (naxis == 3)  nelements*=naxes[2];
     986  int anynull;
     987  double dnullval=0.;
     988  float fnullval=0.;
     989  // on laisse a fits le soin de convertir le type du tableau lu vers
     990  // le type suppose par l'utilisateur de fitsioserver
     991  //
     992  switch ( FITS_tab_typ_)
     993    {
     994    case TDOUBLE :
     995      if (bitpix !=  DOUBLE_IMG)
     996        {
     997          cout << " FitsIOServer : the data type on fits file is not double, "
     998               << " conversion to double will be achieved by cfitsio lib " << endl;
     999        }
     1000      if (dtab_ != NULL) delete [] dtab_;
     1001      dtab_=new double[nelements];
     1002      fits_read_img(fptr, TDOUBLE, 1, nelements,  &dnullval, dtab_,
     1003                    &anynull, &status);
     1004      if( status )  printerror( status );
     1005      break;
     1006    case TFLOAT :
     1007      if (bitpix !=  FLOAT_IMG)
     1008        {
     1009          cout << " FitsIOServer : the data type on fits file is not float, "
     1010               << " conversion to float will be achieved by cfitsio lib " << endl;
     1011        }
     1012      if (ftab_ != NULL) delete [] ftab_;
     1013      ftab_=new float[nelements];
     1014      fits_read_img(fptr, TFLOAT, 1, nelements,  &fnullval, ftab_,
     1015                    &anynull, &status);
     1016      if( status )  printerror( status );
     1017      break;
     1018
     1019
     1020    case TINT :
     1021      if (bitpix !=  LONG_IMG)
     1022        {
     1023          cout << " FitsIOServer : the data type on fits file is not long, "
     1024               << " conversion to long will be achieved by cfitsio lib " << endl;
     1025        }
     1026      if (i_4tab_ != NULL) delete [] i_4tab_;
     1027      i_4tab_=new int_4[nelements];
     1028      fits_read_img(fptr, TINT, 1, nelements,  &fnullval, i_4tab_,
     1029                    &anynull, &status);
     1030      if( status )  printerror( status );
     1031      break;
     1032
     1033
     1034    default :
     1035      cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
     1036      break;
     1037    }
     1038  status = 0;
     1039  char card[FLEN_CARD];
     1040  int num = 0;
     1041  char comment2[FLEN_COMMENT] = "x";
     1042  char keyname[]= "";
     1043  char datekey[]= "DATE";
     1044  char endkey[] = "END";
     1045  char typ='x';
     1046  int ival;
     1047  double dval;
     1048  char strval[]="";
     1049  // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
     1050  // on va jusqu'au mot cle DATE
     1051  int flen_keyword = 9;
     1052  // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
     1053  // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
     1054  // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
     1055  while (status == 0 &&  strncmp(keyname, datekey,4) != 0 )
     1056    {
     1057      num++;
     1058      fits_read_record(fptr, num, card, &status);
     1059      strncpy(keyname,card,flen_keyword-1);
     1060    }
     1061  if (status != 0 )
     1062    {
     1063      cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
     1064    }
     1065  // on recupere la liste des mots-cles utilisateurs
     1066  while (status == 0)
     1067    {
     1068      num++;
     1069      // on lit un record pour recuperer le nom du mot-cle
     1070      fits_read_record(fptr, num, card, &status);
     1071      strncpy(keyname,card,flen_keyword-1);
     1072      char value[FLEN_VALUE];
     1073      // on recupere le premier caractere du commentaire, qui contient
     1074      // le code du type de la valeur
     1075      // (tant que ce n est pas le mot cle END)
     1076      fits_read_keyword(fptr, keyname, value, comment2, &status);
     1077      if (  strncmp(keyname, endkey,flen_keyword-1) != 0)
     1078        {
     1079          typ = comment2[0];
     1080          // quand le type est connu, on lit la valeur
     1081          switch (typ)
     1082            {
     1083            case 'I' :
     1084              fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
     1085              if( status )  printerror( status );
     1086              strip (keyname, 'B',' ');
     1087              dvl[keyname] = (int_4)ival;
     1088              break;
     1089            case 'D' :
     1090              fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
     1091              if( status )  printerror( status );
     1092              strip (keyname, 'B',' ');
     1093              dvl[keyname] = dval;
     1094              break;
     1095            case 'S' :
     1096              fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
     1097              if( status )  printerror( status );
     1098              strip (keyname, 'B',' ');
     1099              strip(strval, 'B',' ');
     1100              dvl[keyname]=strval;
     1101              break;
     1102            default :
     1103              cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
     1104              break;
     1105            }
     1106        }
     1107    }
     1108
     1109  // close the file
     1110  status=0;
     1111  fits_close_file(fptr, &status);
     1112  if( status )  printerror( status );
     1113}
     1114
     1115
     1116
     1117// projects a SphericalMap<double>, according sinus-method, and saves onto
     1118// a FITS-file
     1119void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
     1120{
     1121   
     1122  long naxes[2]={600, 300};
     1123  float* map =new float[ 600*300 ];
     1124  int npixels=  naxes[0]*naxes[1];
     1125
     1126   cout << " image FITS en projection SINUS" << endl;
     1127  // table will have npixels rows
     1128  for(int j=0; j < npixels; j++) map[j]=0.;
     1129    for(int j=0; j<naxes[1]; j++)
     1130      {
     1131        double yd = (j+0.5)/naxes[1]-0.5;
     1132        double theta = (0.5-yd)*Pi;
     1133        double facteur=1./sin(theta);
     1134        for(int i=0; i<naxes[0]; i++) 
     1135          {
     1136            double xa = (i+0.5)/naxes[0]-0.5;
     1137            double phi = 2.*Pi*xa*facteur+Pi;
     1138            float th=float(theta);
     1139            float ff=float(phi);
     1140            if (phi<2*Pi && phi>= 0)
     1141              {
     1142                map[j*naxes[0]+i] = sph.PixValSph(th, ff);
     1143              }
     1144          }
     1145      }
     1146 
     1147    write_picture(naxes, map,  filename);
     1148    delete [] map;
     1149}
     1150
     1151// projects a SphericalMap<double>, according sinus-method, and saves onto
     1152// a FITS-file
     1153void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
     1154{
     1155  // le code de cete methode duplique celui de la precedente, la seule
     1156  //difference etant le type de sphere en entree. Ces methodes de projection
     1157  // sont provisoires, et ne servent que pour les tests. C est pourquoi je
     1158  // ne me suis pas casse la tete, pour l instant
     1159   
     1160  long naxes[2]={600, 300};
     1161  float* map = new float[ 600*300 ];
     1162  int npixels=  naxes[0]*naxes[1];
     1163
     1164   cout << " image FITS en projection SINUS" << endl;
     1165  // table will have npixels rows
     1166  for(int j=0; j < npixels; j++) map[j]=0.;
     1167    for(int j=0; j<naxes[1]; j++)
     1168      {
     1169        double yd = (j+0.5)/naxes[1]-0.5;
     1170        double theta = (0.5-yd)*Pi;
     1171        double facteur=1./sin(theta);
     1172        for(int i=0; i<naxes[0]; i++) 
     1173          {
     1174            double xa = (i+0.5)/naxes[0]-0.5;
     1175            double phi = 2.*Pi*xa*facteur+Pi;
     1176            float th=float(theta);
     1177            float ff=float(phi);
     1178            if (phi<2*Pi && phi>= 0)
     1179              {
     1180                map[j*naxes[0]+i] = sph.PixValSph(th, ff);
     1181              }
     1182          }
     1183      }
     1184 
     1185    write_picture(naxes, map,  filename);
     1186    delete [] map;
     1187
     1188}
     1189
     1190// saves a (LocalMap<double> on a FITS-file in order to be visualized
     1191// (for tests)
     1192void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
     1193{
     1194   
     1195  long naxes[2];
     1196  naxes[0] = lcm.Size_x();
     1197  naxes[1] = lcm.Size_y();
     1198  int npixels=  naxes[0]*naxes[1];
     1199  float* map = new float[npixels];
     1200
     1201  // table will have npixels rows
     1202  for(int j=0; j < npixels; j++) map[j]=0.;
     1203    for(int j=0; j<naxes[1]; j++)
     1204      {
     1205        for(int i=0; i<naxes[0]; i++) 
     1206          {
     1207                map[j*naxes[0]+i] = lcm(i, j);
     1208          }
     1209      }
     1210 
     1211    write_picture(naxes, map,  filename);
     1212    delete [] map;
     1213}
     1214
     1215void  FitsIoServer::read_sphere(char flnm[], int& npixels, int& nside)
     1216{
     1217  long naxis;
     1218  nside=0;
     1219  int n1, n2, n3;
     1220  DVList dvl;
     1221  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
     1222
     1223  npixels=n1;
     1224  //dvl.Print();
     1225  DVList::ValList::const_iterator it;
     1226  for(it = dvl.Begin(); it != dvl.End(); it++) 
     1227    {
     1228      char datatype= (*it).second.typ;
     1229      char keyname[]="";
     1230      strcpy(keyname, (*it).first.substr(0,64).c_str());   
     1231      int ival=0;
     1232      double dval=0.;
     1233      char strval[]="";
     1234      switch (datatype)
     1235        {
     1236        case 'I' :
     1237          ival=(*it).second.mtv.iv;
     1238          break;
     1239        case 'D' :
     1240          dval=(*it).second.mtv.dv;
     1241          break;
     1242        case 'S' :
     1243          strcpy(strval, (*it).second.mtv.strv);
     1244          break;
     1245        default :
     1246          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
     1247          break;
     1248        }
     1249      if (!strcmp(keyname, "NSIDE") )
     1250        {
     1251          nside = ival;
     1252        }
     1253      //if (!strcmp(keyname, "ORDERING") )
     1254      //        {
     1255      //          cout << " j'ai trouve ORDERING " << endl;
     1256      //        }
     1257    }
     1258
     1259  if (naxis != 1)
     1260    {
     1261      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
     1262    }
     1263}
     1264
     1265
     1266void  FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
     1267{
     1268 
     1269  int bitpix = FLOAT_IMG;
     1270  long naxis = 2;
     1271
     1272  //pointer to the FITS file, defined in fitsio.h
     1273  fitsfile *fptr;
     1274  // delete old file if it already exists
     1275  remove(filename);
     1276    // initialize status before calling fitsio routines
     1277  int status = 0;         
     1278
     1279  // create new FITS file 
     1280  fits_create_file(&fptr, filename, &status);
     1281  if( status )  printerror( status );
     1282
     1283  // write the required header keywords
     1284  fits_create_img(fptr, bitpix, naxis, naxes, &status);
     1285  if( status )  printerror( status );
     1286
     1287  // write the current date
     1288  fits_write_date(fptr, &status);
     1289  if( status )  printerror( status );
     1290
     1291
     1292  // first row in table to write
     1293  long firstrow  = 1;
     1294  // first element in row   
     1295  long firstelem = 1; 
     1296  int colnum = 1;
     1297  int nelements=naxes[0]*naxes[1];
     1298  fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
     1299  if( status )  printerror( status );
    14431300
    14441301  // close the file
    14451302  fits_close_file(fptr, &status);
    14461303  if( status )  printerror( status );
     1304}
     1305
     1306
     1307bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
     1308
     1309  //*****************************************************/
     1310  //* check if the specified keyword exits in the CHU   */
     1311  //*****************************************************/
     1312{
     1313
     1314  bool KEY_EXIST = false;
     1315  int status = 0;
     1316  char strbide[FLEN_VALUE];
     1317  char keybide[LEN_KEYWORD]= "";
     1318  for(int jj = 1; jj <= nkeys; jj++)
     1319    {
     1320      if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
     1321        printerror( status );
     1322      if( !strcmp(keybide,keyword) )
     1323        {
     1324          KEY_EXIST= true;
     1325          break;
     1326        }
     1327    }
     1328  return(KEY_EXIST);
     1329}
     1330
     1331void FitsIoServer::readheader ( char filename[] )
     1332
     1333  //**********************************************************************/
     1334  //* Print out all the header keywords in all extensions of a FITS file */
     1335  //**********************************************************************/
     1336{
     1337
     1338  // standard string lengths defined in fitsioc.h
     1339  char card[FLEN_CARD];
     1340
     1341  // pointer to the FITS file, defined in fitsio.h
     1342  fitsfile *fptr;
     1343
     1344  int status = 0;
     1345  if ( fits_open_file(&fptr, filename, READONLY, &status) )
     1346    printerror( status );
    14471347 
    1448   if (naxis != 2)
    1449     {
    1450       cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
    1451     }
    1452 
    1453   DpcImg.Allocate(siz_x, siz_y, 0., 0);
    1454   int_4* pixelPtr=DpcImg.ImagePtr();
    1455 
    1456 
    1457   // verifications de type
    1458   PBaseDataTypes  dataT=DataType((int_4)0);
    1459   int TypeDonnees=dvl.GetI("DATATYPE");
    1460   if (int(dataT) !=  TypeDonnees)
    1461     {
    1462       cout << " FitsIOServer : parameter DATATYPE on file is not int_4 " << endl;
    1463       cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
    1464     }
    1465   memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
    1466   const char* nom=dvl.GetS("NAME").c_str();
    1467   DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
    1468   DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
    1469   DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
    1470   DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
    1471        dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
    1472        dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
    1473 }
    1474 
    1475 
     1348  // attempt to move to next HDU, until we get an EOF error
     1349  int hdutype;
     1350  for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
     1351    {
     1352      if (hdutype == ASCII_TBL)
     1353        printf("\nReading ASCII table in HDU %d:\n", ii);
     1354      else if (hdutype == BINARY_TBL)
     1355        printf("\nReading binary table in HDU %d:\n", ii);
     1356      else if (hdutype == IMAGE_HDU)
     1357        printf("\nReading FITS image in HDU %d:\n", ii);
     1358      else
     1359        {
     1360          printf("Error: unknown type of this HDU \n");
     1361          printerror( status );
     1362        }
     1363
     1364      // get the number of keywords
     1365      int nkeys, keypos;
     1366      if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
     1367        printerror( status );
     1368     
     1369      printf("Header listing for HDU #%d:\n", ii);
     1370      for (int jj = 1; jj <= nkeys; jj++) 
     1371        {
     1372          if ( fits_read_record(fptr, jj, card, &status) )
     1373            printerror( status );
     1374
     1375          // print the keyword card
     1376          printf("%s\n", card);
     1377        }
     1378      printf("END\n\n");
     1379    }
     1380
     1381  // got the expected EOF error; reset = 0
     1382  if (status == END_OF_FILE)   
     1383    status = 0;             
     1384  else
     1385    printerror( status );
     1386
     1387  if ( fits_close_file(fptr, &status) )
     1388    printerror( status );
     1389 
     1390  return;
     1391}
     1392
     1393void FitsIoServer::printerror(int status) const
     1394
     1395  //*****************************************************/
     1396  //* Print out cfitsio error messages and exit program */
     1397  //*****************************************************/
     1398{
     1399
     1400  // print out cfitsio error messages and exit program
     1401  if( status )
     1402    {
     1403      // print error report 
     1404      fits_report_error(stderr, status);
     1405      // terminate the program, returning error status
     1406      exit( status );
     1407    }
     1408  return;
     1409}
     1410
     1411
     1412
Note: See TracChangeset for help on using the changeset viewer.