Changeset 603 in Sophya for trunk/SophyaExt

Timestamp:

Nov 19, 1999, 6:33:34 PM (26 years ago)

Author:

ansari

Message:

ajout save sphereGorski

Location:

trunk/SophyaExt/FitsIOServer

Files:

• fitsioserver.cc (modified) (13 diffs)
• fitsioserver.h (modified) (3 diffs)

trunk/SophyaExt/FitsIOServer/fitsioserver.cc

@@ -333 +333 @@
         }
     }
+  //  cout << " fin relecture debut transfert" << endl;
   for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
+}
+void FitsIoServer::load(SphereGorski<double>& sph, char flnm[], int hdunum) 
+{
+  int npixels=0;
+  int nside=0;
+
+  FITS_tab_typ_ = TDOUBLE;
+
+  DVList dvl;
+  planck_read_bntbl(flnm, hdunum, npixels, dvl);
+  //dvl.Print();
+  nside= dvl.GetI("NSIDE");
+  const char* ordering= dvl.GetS("ORDERING").c_str();
+  char* ring = "'RING";
+  if (  strncmp(ordering, ring,3) != 0) 
+
+    //  if (ordering!="RING    ")
+    {
+      cout << " numerotation non RING" << endl;
+    }
+  // number of pixels in the sphere
+  if (sph.NbPixels() != npixels)
+    {
+      cout << " found " << npixels << " pixels" << endl;
+      cout << " expected " << sph.NbPixels() << endl;
+      if (nside <= 0 ) 
+        {
+          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
+          exit(0);
+        }
+      if (nside != sph.SizeIndex())
+        { 
+          sph.Resize(nside);
+          cout << " resizing the sphere to nside=  " << nside  << endl;
+        }
+      else
+        {
+          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
+          exit(0);
+        }
+    }
+  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
 }
 
@@ -731 +774 @@
 
 }
+void FitsIoServer::save(SphereGorski<float>& sph, char filename[])   
+{
+  int npixels = sph.NbPixels();
+  FITS_tab_typ_ = TFLOAT;
+  if (r_4tab_ != NULL ) delete[] r_4tab_;
+  r_4tab_=new r_4[npixels];
+
+
+  for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
+  DVList dvl;
+  dvl.SetS("PIXTYPE","HEALPIX ");
+  dvl["ORDERING"]=sph.TypeOfMap();
+  dvl["NSIDE"] = (int_4)sph.SizeIndex();
+  dvl["FIRSTPIX"]=0;
+  dvl["LASTPIX"]=npixels-1;
+  char* typeOfContent="TEMPERATURE";
+  char* extname="SIMULATION";
+  char* comment1="     Sky Map Pixelisation Specific Keywords";
+  planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
+  delete[] r_4tab_;
+
+}
+void FitsIoServer::save(SphereGorski<double>& sph, char filename[])   
+{
+  int npixels = sph.NbPixels();
+  FITS_tab_typ_ = TDOUBLE;
+  if (r_8tab_ != NULL ) delete[] r_8tab_;
+  r_8tab_=new r_8[npixels];
+
+
+  for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
+  DVList dvl;
+  dvl.SetS("PIXTYPE","HEALPIX ");
+  dvl["ORDERING"]=sph.TypeOfMap();
+  dvl["NSIDE"] = (int_4)sph.SizeIndex();
+  dvl["FIRSTPIX"]=0;
+  dvl["LASTPIX"]=npixels-1;
+  char* typeOfContent="TEMPERATURE";
+  char* extname="SIMULATION";
+  char* comment1="     Sky Map Pixelisation Specific Keywords";
+  planck_write_bntbl(filename, npixels, typeOfContent, extname, comment1, dvl);
+  delete[] r_8tab_;
+
+}
 
 void FitsIoServer::save(LocalMap<double>& locm, char filename[])   
@@ -983 +1070 @@
       if( status ) printerror( status );
     }
+  // close the file
+  fits_close_file(fptr, &status);
+  if( status )  printerror( status );
+}
+
+void FitsIoServer::planck_write_bntbl(char flnm[], int npixels, char typeOfContent[], char extname[], char comment1[], DVList& dvl)
+{
+
+
+  // pointer to the FITS file, defined in fitsio.h
+  fitsfile *fptr;
+
+  // initialize status before calling fitsio routines
+  int status = 0;         
+
+  // delete old file if it already exists
+  remove(flnm);
+
+  // create new FITS file  
+  fits_create_file(&fptr, flnm, &status); 
+  if( status )  printerror( status );
+
+  // primary header
+  int bitpix=LONG_IMG;
+  int naxis=0;
+  fits_create_img(fptr, bitpix, naxis, NULL, &status);
+  // write the current date
+  fits_write_date(fptr, &status);
+  if( status )  printerror( status );
+
+
+  long nrows=npixels/1024;
+  if (nrows==0) nrows=1;
+  if (npixels%1024 !=0) nrows++;
+  int tfields=1;
+  char **ttype, **tform;
+  ttype= new char*[tfields];
+  tform= new char*[tfields];
+  char* format;
+  if ( FITS_tab_typ_ == TDOUBLE) format="1024D";
+  if ( FITS_tab_typ_ == TFLOAT)  format="1024E";
+  if ( FITS_tab_typ_ == TINT)    format="1024I";
+  for(int i = 0; i < tfields; i++)
+    {
+      ttype[i]= new char[FLEN_VALUE];
+      strcpy(ttype[i],  typeOfContent);
+      tform[i]= new char[FLEN_VALUE];
+      strcpy(tform[i], format);
+    }
+  // create a new empty binary table onto the FITS file
+  // physical units if they exist, are defined in the DVList object
+  // so the null pointer is given for the tunit parameters.
+ fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
+                 NULL,extname,&status);
+  if( status )  printerror( status );
+  for( int ii = 0; ii < tfields; ii++)
+    {
+      delete [] ttype[ii];
+      delete [] tform[ii];
+    }
+  delete [] ttype;
+  delete [] tform;
+  long nelements= npixels;
+  switch (FITS_tab_typ_) 
+    {
+    case TDOUBLE :
+      fits_write_col(fptr, TDOUBLE, 1, 1, 1, nelements, r_8tab_, &status);
+      if( status )  printerror( status, "probleme dans ecriture du tableau de doubles" );
+      break;
+
+
+    case TFLOAT :
+      for (int kk=0; kk<10; kk++) cout << r_4tab_[kk] << endl;
+      fits_write_col(fptr, TFLOAT, 1, 1, 1, nelements, r_4tab_, &status);
+      if( status ) printerror( status );
+      break;
+    case TINT :
+      fits_write_col(fptr, TINT, 1, 1, 1, nelements, i_4tab_, &status);
+      if( status ) printerror( status );
+      break;
+    default :
+      cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
+      break;
+    }
+  // write supplementary keywords
+  fits_write_comment(fptr, " ", &status);
+  fits_write_comment(fptr, " ", &status);
+  if( status ) printerror( status );
+  fits_write_comment(fptr,"--------------------------------------------", &status);
+  fits_write_comment(fptr, comment1, &status);
+  fits_write_comment(fptr,"--------------------------------------------", &status);
+  if( status ) printerror( status );
+
+
+  int flen_keyword = 9;
+  // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une 
+  // erreur dans la librairie fits qui donne FLEN_KEYWORD=72 
+  // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
+  DVList::ValList::const_iterator it;
+  for(it = dvl.Begin(); it != dvl.End(); it++)  
+    {
+      int datatype=  key_type_PL2FITS( (*it).second.typ);
+      char keyname[]="";
+      strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
+      char comment[FLEN_COMMENT]="";
+      int ival=0;
+      double dval=0.;
+      char strval[]="";
+      switch (datatype)
+        {
+        case TINT :
+          ival=(*it).second.mtv.iv;
+          strcpy(comment," ");
+          fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
+          break;
+        case TDOUBLE :
+          dval=(*it).second.mtv.dv;
+          strcpy(comment," ");
+          fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
+          break;
+        case TSTRING :
+          strcpy(strval, (*it).second.mtv.strv);
+          strcpy(comment," ");
+          fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
+          break;
+        default :
+          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
+          break;
+        }
+      if( status ) printerror( status );
+    }
+  //char keyname[]="";
+  //char strval[]="";
+  //char comment[FLEN_COMMENT]="";
+  //strncpy(keyname, "CREATOR",flen_keyword-1);
+  //strcpy(strval, "SOPHYA");
+  //strcpy(comment," Orsay");
+  //fits_write_key(fptr, TSTRING, keyname, &strval, comment, &status);
+  //if( status ) printerror( status );
   // close the file
   fits_close_file(fptr, &status);
@@ -1211 +1437 @@
   //cout << " width : " << width << endl;
   fits_read_key_lng(fptr, "LASTPIX", &lastpix, comment, &status);
-  if( status )  printerror( status );
+  if( status )  printerror( status," mot cle LASTPIX" );
 
   long nelements= nrows*repeat;
@@ -1240 +1466 @@
                     r_8tab_, 
                     &anynull, &status);
-      if( status )  printerror( status );
+      if( status )  printerror( status, "probleme dans lecture du tableau de doubles" );
       break;
     case TFLOAT :
@@ -1253 +1479 @@
       fits_read_col(fptr, TFLOAT, 1, 1, 1, nelements,  &fnullval, 
                     r_4tab_, &anynull, &status);
-      if( status )  printerror( status );
+      if( status )  printerror( status,"probleme dans lecture du tableau de floats" );
       break;
 
@@ -1269 +1495 @@
       //fits_read_img(fptr, TINT, 1, nelements,  &inullval, i_4tab_, 
       //                    &anynull, &status);
-  //if( status )  printerror( status );
+      if( status )  printerror( status,"probleme dans lecture du tableau de ints" );
       break;
 
@@ -1297 +1523 @@
           {
             fits_get_keytype(strval,&dtype,&status);
-            //       cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
+            //  cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
             strip (keyname, 'B',' '); 
             strip(strval, 'B',' ');
@@ -1316 +1542 @@
                 break;
               default :
-                cout << " FitsIOServer : mot-cle bizerre " << endl;
+                cout << " FitsIOServer : mot-cle bizarre " << endl;
                 break;
               }
@@ -1402 +1628 @@
 }
 
-// projects a SphericalMap<double>, according Mollweide-method, and saves onto 
+// projects a SphericalMap<float>, according Mollweide-method, and saves onto 
 // a FITS-file
 void FitsIoServer::Mollweide_picture_projection(SphericalMap<float>& sph, char filename[])
@@ -1440 +1666 @@
    
 }
+// projects a SphericalMap<double>, according Mollweide-method, and saves onto 
+// a FITS-file
+void FitsIoServer::Mollweide_picture_projection(SphericalMap<double>& sph, char filename[])
+{
+  // le code de cete methode duplique celui de la precedente, la seule 
+  //difference etant le type de sphere en entree. Ces methodes de projection
+  // sont provisoires, et ne servent que pour les tests. C est pourquoi je 
+  // ne me suis pas casse la tete, pour l instant
+    
+  long naxes[2]={600, 300};
+  float* map = new float[ 600*300 ];
+  int npixels=  naxes[0]*naxes[1];
+
+   cout << " image FITS en projection MOLLWEIDE" << endl;
+  // table will have npixels rows
+   for(int j=0; j < npixels; j++) map[j]=0.;
+   for(int j=0; j<naxes[1]; j++) 
+     {
+       double yd = (j+0.5)/naxes[1]-0.5;
+       double facteur=2.*Pi/sin(acos(yd*2));
+       double theta = (0.5-yd)*Pi;
+       for(int i=0; i<naxes[0]; i++)  
+         {
+           double xa = (i+0.5)/naxes[0]-0.5;
+           double phi =  xa*facteur+Pi;
+           float th=float(theta);
+           float ff=float(phi);
+           if (phi<2*Pi && phi>= 0) 
+             {
+               map[j*naxes[0]+i] = sph.PixValSph(th, ff);
+             }
+         }
+     }
+   
+   write_picture(naxes, map,  filename);
+   delete [] map;
+   
+}
 
 
@@ -1597 +1861 @@
 }
 
-void FitsIoServer::printerror(int status) const
+void FitsIoServer::printerror(int& status) const
 
   //*****************************************************/
@@ -1605 +1869 @@
 
   // print out cfitsio error messages and exit program 
-  if( status )
-    {
       // print error report  
       fits_report_error(stderr, status);
       // terminate the program, returning error status 
-      exit( status ); 
-    }
-  return;
-}
-
-
-
+      //      exit( status ); 
+      status=0;
+}
+void FitsIoServer::printerror(int& status, char* texte) const
+
+  //*****************************************************/
+  //* Print out cfitsio error messages and exit program */
+  //*****************************************************/
+{
+
+  // print out cfitsio error messages and exit program 
+      // print error report  
+      fits_report_error(stderr, status);
+      cout << " erreur : " << texte << endl;
+      status=0;
+}
+
+
+

trunk/SophyaExt/FitsIOServer/fitsioserver.h

@@ -35 +35 @@
  void load(ImageI4& DpcImg,char flnm[]);
  void load(SphereGorski<float>& sph, char flnm[], int hdunum); 
+ void load(SphereGorski<double>& sph, char flnm[], int hdunum); 
  void save(TMatrix<double>& mat, char flnm[]);
  void save(NTuple& ntpl,char flnm[]);
  void save(SphericalMap<double>& sph, char flnm[]);
  void save(SphericalMap<float>& sph, char flnm[]);
+ void save(SphereGorski<float>& sph, char filename[]);   
+ void save(SphereGorski<double>& sph, char filename[]);   
  void save(LocalMap<double>& locm, char flnm[]);
  void save(const ImageR4& DpcImg,char flnm[]);
@@ -45 +48 @@
  void sinus_picture_projection(SphericalMap<float>& sph, char flnm[]);
  void Mollweide_picture_projection(SphericalMap<float>& sph, char flnm[]);
+ void Mollweide_picture_projection(SphericalMap<double>& sph, char flnm[]);
  void picture(LocalMap<double>& lcm, char flnm[]);
 
@@ -60 +64 @@
  void planck_read_img(char flnm[],long& naxis,int& n1,int& n2,int& n3,DVList& dvl);
  void planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl);
+ void planck_write_bntbl(char flnm[], int npixels, char typeOfContent[], char extname[], char comment1[], DVList& dvl);
  void write_picture(long* naxes,float* map,char* filename) const;
 
- void printerror(int status) const;
-
+ void printerror(int& status) const;
+ void printerror(int& status, char* texte) const;
  bool check_keyword(fitsfile*,int,char keyword[]);