Changeset 886 in Sophya for trunk/SophyaProg

Timestamp:

Apr 11, 2000, 6:11:05 PM (25 years ago)

Author:

ansari

Message:

Sophie: change to FITS_SphereHealpix

File:

• trunk/SophyaProg/PMixer/skymixer.cc (modified) (3 diffs)

trunk/SophyaProg/PMixer/skymixer.cc

@@ -17 +17 @@
 #include "spherehealpix.h"
 #include "fiospherehealpix.h"
-
+#include "fitsspherehealpix.h"
 #include "radspecvector.h"
 #include "blackbody.h"
@@ -101 +101 @@
 so = NULL;
 
-try {
-  string dcard = arg[1];
-  cout << " Decoding parameters from file " << dcard << endl;
-  dc.ReadFile(dcard);
-
-  rc = CheckCards(dc, msg);
-  if (rc) { 
-    cerr << " Error condition -> Rc= " << rc << endl;
-    cerr << " Msg= " << msg << endl;
-    return(rc);
-  }
-}
-catch (PException exc) {
-  msg = exc.Msg();
-  cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
-  return(90);
-}   
+ try {
+   string dcard = arg[1];
+   cout << " Decoding parameters from file " << dcard << endl;
+   dc.ReadFile(dcard);
+  
+   rc = CheckCards(dc, msg);
+   if (rc) { 
+     cerr << " Error condition -> Rc= " << rc << endl;
+     cerr << " Msg= " << msg << endl;
+     return(rc);
+   }
+ }
+ catch (PException exc) {
+   msg = exc.Msg();
+   cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
+   return(90);
+ }   
 
 
@@ -127 +127 @@
 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
 
-SphereHEALPix<float> outgs(hp_nside);
-bool okout = false;
-
-try {
-  if (rdmap) {  // Reading map from FITS file
-    FitsIoServer fios;
-    char ifnm[256];
-    strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
-    ifnm[255] = '\0';
-    cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
-    fios.load(outgs, ifnm, 2);
-    cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
-    if (printlev > 0) {
-      MeanSig(outgs.DataBlock(), moy, sig );
-      cout << " MeanSig for outpout map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-  }
-  else {
-    cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
-    outgs.SetPixels(0.);
-  }
-
-  // Decoding detection pass-band filter 
-  sr = getSpectralResponse(dc);
-  PrtTim(" After FilterCreation ");
-
-  FitsIoServer fios; // Our FITS IO Server
-  char * flnm, buff[90];
-  string key;
-
-  SphereHEALPix<float> ings(hp_nside);  // The input map
-  double K = 1.;
-  double freqOfMap = 1.;
-  // Loop over sky component 
-  int sk;
-  for(sk = 0; sk<nskycomp; sk++) {
-    cout << " Processing sky component No " << sk+1 << endl;
+ SphereHEALPix<float> outgs(hp_nside);
+ try{
+   if (rdmap) {  // Reading map from FITS file
+     char ifnm[256];
+     strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
+     ifnm[255] = '\0';
+     cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
+     FITS_SphereHEALPix<float> fios(ifnm);
+     outgs = (SphereHEALPix<float>)fios;
+      cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
+      if (printlev > 0) {
+        MeanSig(outgs.DataBlock(), moy, sig );
+        cout << " MeanSig for outpout map - Mean= " << moy << " Sigma= " << sig << endl;
+      }
+    }
+    else {
+      cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
+      outgs.SetPixels(0.);
+    }
     
-    sprintf(buff, "%d", sk+1);
-    key = (string)"MAPFITSFILE" + buff;
-    flnm = BuildFITSFileName(dc.SParam(key, 0));
+    // Decoding detection pass-band filter 
+    sr = getSpectralResponse(dc);
+    PrtTim(" After FilterCreation ");
     
-    K = dc.DParam(key, 1, 1.);
+    char * flnm, buff[90];
+    string key;
     
-    bool mapDependentOfFreq = false;
-    key = (string)"BETAFITSFILE"+ buff;
-    if(dc.HasKey(key)) 
-      {
-        mapDependentOfFreq = true;
+    //    SphereHEALPix<float> ings(hp_nside);  // The input map
+    double K = 1.;
+    double freqOfMap = 1.;
+    // Loop over sky component 
+    int sk;
+    for(sk = 0; sk<nskycomp; sk++) {
+      cout << " Processing sky component No " << sk+1 << endl;
+      
+      sprintf(buff, "%d", sk+1);
+      key = (string)"MAPFITSFILE" + buff;
+      flnm = BuildFITSFileName(dc.SParam(key, 0));
+      
+      K = dc.DParam(key, 1, 1.);
+      
+      bool mapDependentOfFreq = false;
+      key = (string)"BETAFITSFILE"+ buff;
+      if(dc.HasKey(key)) 
+        {
+          mapDependentOfFreq = true;
+        }
+      cout << "the map has a mapDependentOfFreq = " << 
+        mapDependentOfFreq << endl;
+      
+      cout << " Reading Input FITS map " << (string)flnm << endl;
+      FITS_SphereHEALPix<float> fiosIn(flnm);
+      SphereHEALPix<float> ings = (SphereHEALPix<float>)fiosIn;
+      
+      if (debuglev > 4) {  // Writing the input map to the outppf
+        FIO_SphereHEALPix<float> fiog(ings);
+        fiog.Write(*so, key);
       }
-    cout << "the map has a mapDependentOfFreq = " << 
-      mapDependentOfFreq << endl;
-
-    cout << " Reading Input FITS map " << (string)flnm << endl;
-    fios.load(ings, flnm, 2);
-    if (debuglev > 4) {  // Writing the input map to the outppf
-      FIO_SphereHEALPix<float> fiog(ings);
-      fiog.Write(*so, key);
-    }
-    if (printlev > 2) {
-      MeanSig(ings.DataBlock(), moy, sig );
-      cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-    // getting Emission spectra  
-    if(!mapDependentOfFreq)
-      {
-        if (es) { delete es; es = NULL; }
-        es = getEmissionSpectra(dc, sk);
-        addComponent(*sr, outgs, ings, *es, K);
+      if (printlev > 2) {
+        MeanSig(ings.DataBlock(), moy, sig );
+        cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
       }
-    else
-      {
-        key = (string)"BETAFITSFILE"+ buff;
-        SphereHEALPix<float> betaMap;
-        flnm = BuildFITSFileName(dc.SParam(key, 0));
-        double normFreq = dc.DParam(key, 1, 1.);
-        if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
-        int nSideForInt = dc.DParam(key, 2, 1.);
-        if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
-        cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
-        fios.load(betaMap, flnm, 2);
-        if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
-
-        bool bydefault = true;
-        if(!bydefault)
-          addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
-        else
-          {
-            // integrate the betamap over the SpectralResponse
-            SphereHEALPix<float> intBetaMap(nSideForInt);
-            integratedMap(*sr, betaMap, normFreq, intBetaMap);
-            betaMap.Resize(8);
-            MeanSig(intBetaMap.DataBlock(), moy, sig );
-            if (printlev > 4) cout << "....BetaFits... MeanSig for intBetaMap - Mean= " << moy << " Sigma= " << sig << endl;
-            // add the integrated beta map 
-            addComponentBeta(outgs,ings,intBetaMap, K);
-          }
+      // getting Emission spectra  
+      if(!mapDependentOfFreq)
+        {
+          if (es) { delete es; es = NULL; }
+          es = getEmissionSpectra(dc, sk);
+          addComponent(*sr, outgs, ings, *es, K);
+        }
+      else
+        {
+          key = (string)"BETAFITSFILE"+ buff;
+          //SphereHEALPix<float> betaMap;
+          flnm = BuildFITSFileName(dc.SParam(key, 0));
+          double normFreq = dc.DParam(key, 1, 1.);
+          if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
+          int nSideForInt = dc.DParam(key, 2, 1.);
+          if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
+          cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
+          
+          FITS_SphereHEALPix<float> fiosBM(flnm);
+          SphereHEALPix<float> betaMap = (SphereHEALPix<float>)fiosBM;
+          
+          if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
+          
+          bool bydefault = true;
+          if(!bydefault)
+            addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
+          else
+            {
+              // integrate the betamap over the SpectralResponse
+              SphereHEALPix<float> intBetaMap(nSideForInt);
+              integratedMap(*sr, betaMap, normFreq, intBetaMap);
+              betaMap.Resize(8);
+              MeanSig(intBetaMap.DataBlock(), moy, sig );
+              if (printlev > 4) cout << "....BetaFits... MeanSig for intBetaMap - Mean= " << moy << " Sigma= " << sig << endl;
+              // add the integrated beta map 
+              addComponentBeta(outgs,ings,intBetaMap, K);
+            }
+          MeanSig(outgs.DataBlock(), moy, sig );
+          cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
+        }
+      
+      if (printlev > 1) {
         MeanSig(outgs.DataBlock(), moy, sig );
         cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
       }
-
-    okout = true;
-    if (printlev > 1) {
-      MeanSig(outgs.DataBlock(), moy, sig );
-      cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-    sprintf(buff, "End of Proc. Comp. %d ", sk+1);
-    cout << " --------------------------------------- \n" << endl;    
-    PrtTim(buff);
-  }   // End of sky component loop 
-}   // End of try block
-
-catch (PException exc) {
-  msg = exc.Msg();
-  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
-  rc = 50;
-  }   
-
-// Saving the output map in FITS format 
-if (okout) {
-  FitsIoServer fios;
-  fios.save(outgs, arg[2]);
-  cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
-       << (string)(arg[2]) << endl;
-  PrtTim("End of WriteFITS ");
-  // Saving the output map in PPF format 
-  if (narg > 3) {
+      sprintf(buff, "End of Proc. Comp. %d ", sk+1);
+      PrtTim(buff);
+    } 
+ }
+ catch(PException exc) 
+   {
+     cout << "catched PException" << endl;
+     msg = exc.Msg();
+     cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
+     rc = 50;
+     return(50);
+   } 
+ 
+ // Saving the output map in FITS format 
+ cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
+      << (string)(arg[2]) << endl;
+ FITS_SphereHEALPix<float> fios2(outgs);
+ fios2.Write(arg[2]);
+ PrtTim("End of WriteFITS ");
+ // Saving the output map in PPF format 
+ if (narg > 3) {
    POutPersist s(arg[3]); 
    FIO_SphereHEALPix<float> fiog(&outgs) ;
    fiog.Write(s);
    cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "  
-        << (string)(arg[3]) << endl;
+        << (string)(arg[3]) << endl;
    PrtTim("End of WritePPF ");
-   }
-  }
-  if (so) delete so;  //  Closing the debug ppf file 
-  return(rc);
+ }
+ if (so) delete so;  //  Closing the debug ppf file 
+ return(rc);
 }