Changeset 876 in Sophya

Timestamp:

Apr 11, 2000, 3:23:03 PM (25 years ago)

Author:

ansari

Message:

sophie: new functionnalities added according to Garching meeting:
:

File:

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

trunk/SophyaProg/PMixer/skymixer.cc

@@ -15 +15 @@
 #include "fitsioserver.h"
 
-#include "spherehealpix.h"
-#include "fiospherehealpix.h"
+#include "spheregorski.h"
 
 #include "radspecvector.h"
 #include "blackbody.h"
+#include "derivblackbody.h"
 #include "nupower.h"
 
@@ -36 +36 @@
 void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
 
+// 
+// treating maps
+//---------------------------------------------------------
+template <class T> void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
+template <class T> void Sph2Sph(PixelMap<T>& in, PixelMap<T>& out);
+//
+// to add different sky components and corresponding tools
+//----------------------------------------------------------
+template <class T> 
+void addComponent(SpectralResponse&  sr, 
+                                     PixelMap<T>& finalMap, 
+                                     PixelMap<T>& mapToAdd, 
+                                     RadSpectra& rs, double K=1.);
+//
 template <class T>
-void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
-                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K=1.);
+void addComponentBeta(SphereGorski<T>& finalMap, 
+                      SphereGorski<T>& mapToAdd,SpectralResponse& sr, 
+                      SphereGorski<T>& betaMap, double normFreq, double K);
+//
 template <class T>
-void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
+void integratedMap(SpectralResponse&  sr,   
+                   SphereGorski<T>& betaMap, double normFreq, SphereGorski<T>& intBetaMap);
+  //template <class T> 
+  //void integratedMap(SpectralResponse&  sr,   
+  //                          PixelMap<T>& betaMap, 
+  //                          double normFreq, 
+  //                          PixelMap<T>& outMap);
+template <class T>
+void addComponentBeta(SphereGorski<T>& finalMap, 
+                      SphereGorski<T>& mapToAdd,
+                      SphereGorski<T>& intBetaMap, double K);
+//
 // -----------------------------------------------------------------
 
@@ -92 +119 @@
 
 
-cout << " skymix/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
+cout << " skymix/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
 cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev 
      << "  PrintLev= " << printlev << endl;
@@ -99 +126 @@
 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
 
-SphereHEALPix<float> outgs(hp_nside);
+SphereGorski<float> outgs(hp_nside);
 bool okout = false;
 
@@ -129 +156 @@
   string key;
 
-  SphereHEALPix<float> ings(hp_nside);  // The input map
+  SphereGorski<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;
     fios.load(ings, flnm, 2);
-    K = dc.DParam(key, 1, 1.);
-    if (debuglev > 4) {  // Writing tne input map to the outppf
-      FIO_SphereHEALPix<float> fiog(ings);
+    if (debuglev > 4) {  // Writing the input map to the outppf
+      FIO_SphereGorski<float> fiog(ings);
       fiog.Write(*so, key);
-      }
-    // getting Emission spectra  
-    if (es) { delete es; es = NULL; }
-    es = getEmissionSpectra(dc, sk);
-    
+    }
     if (printlev > 2) {
       MeanSig(ings.DataBlock(), moy, sig );
       cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
     }
-    addComponent(*sr, outgs, ings, *es, 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;
+        SphereGorski<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(betaMap.NbPixels()/12);
+
+        bool bydefault = true;
+        if(!bydefault)
+          addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
+        else
+          {
+            // integrate the betamap over the SpectralResponse
+            SphereGorski<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;
+      }
+
     okout = true;
     if (printlev > 1) {
@@ -175 +248 @@
   FitsIoServer fios;
   fios.save(outgs, arg[2]);
-  cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
+  cout << "Output Map (SphereGorski<float>) written to FITS file " 
        << (string)(arg[2]) << endl;
   PrtTim("End of WriteFITS ");
-// Saving the output map in PPF format 
+  // Saving the output map in PPF format 
   if (narg > 3) {
    POutPersist s(arg[3]); 
-   FIO_SphereHEALPix<float> fiog(&outgs) ;
+   FIO_SphereGorski<float> fiog(&outgs) ;
    fiog.Write(s);
-   cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "  
+   cout << "Output Map (SphereGorski<float>) written to POutPersist file "  
         << (string)(arg[3]) << endl;
    PrtTim("End of WritePPF ");
@@ -242 +315 @@
 
   // Checking detection filter specification
-//  Checking input FITS file specifications 
+  //  Checking input FITS file specifications 
   int kc;
   char buff[256];
   bool pb = false;
   string key3;
+  string key4;
+  string key5;
   for(kc=0; kc<ncomp; kc++) {
     sprintf(buff, "MAPFITSFILE%d", kc+1);
     key = buff;
-    if (dc.NbParam(key) < 1) {   
-      msg = "Missing or invalid card : " + key;
+    if (dc.NbParam(key) < 1 ) {   
+      msg = "Missing or invalid card : " + key + " " + key2 ;
       pb = true;  break;
     }
@@ -260 +335 @@
     sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
     key3 = buff;
-    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6))  {   
-      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3;
+    sprintf(buff, "BETAFITSFILE%d", kc+1);
+    key4 = buff;
+    sprintf(buff, "DIPOLE%d", kc+1);
+    key5 = buff;
+    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
+         && (dc.NbParam(key5)<1))  {   
+      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
       pb = true;  break;
     }
@@ -352 +432 @@
   string key = (string)"SPECTRAFITSFILE" + numb;
   string key2 = (string)"BLACKBODY" + numb;
+  string key5 = (string)"DIPOLE" + numb;
   string key3 = (string)"POWERLAWSPECTRA" + numb;
   string ppfname = "espectra";
@@ -377 +458 @@
     ppfname = key2;
   }
+  else if (dc.HasKey(key5) ) { // Creating Dipole
+    rs = new DerivBlackBody(dc.DParam(key5, 0, 2.726));
+    ppfname = key5;
+  }
   else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
     double a = dc.DParam(key3, 0, 1.);
@@ -416 +501 @@
     }
 }
+/* Nouvelle-Fonction */
+template <class T> 
+void addComponentBeta(SphereGorski<T>& finalMap, 
+                      SphereGorski<T>& mapToAdd,SpectralResponse& sr, 
+                      SphereGorski<T>& betaMap, double normFreq, double K)
+{
+  // finalMap = finalMap + coeff* mapToAdd
+  // coeff    =  convolution of sr and rs 
+  // compute the coefficient corresponding to mapToAdd
+  // betaMap is the map of (beta(theta,phi)) 
+
+  int nbpix = finalMap.NbPixels();
+  if (nbpix != mapToAdd.NbPixels())   
+    throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
+  if (printlev > 1) 
+    {
+      cout << "addComponentBeta - Coeff= "  << K  << endl; 
+      cout << "nb pixels: " << finalMap.NbPixels() << endl;
+    }
+  SphereGorski<T> bigBetaMap(sqrt(nbpix/12));
+  if(nbpix != betaMap.NbPixels())
+    {
+       Sph2Sph(betaMap,bigBetaMap);
+     }
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
+      RadSpectra* rs =  new PowerLawSpectra
+        (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
+      double coeff = rs->filteredIntegratedFlux(sr);
+      finalMap(i) += coeff*K*mapToAdd(i);
+    }
+}
+
+template <class T>
+void integratedMap(SpectralResponse&  sr,   
+                   SphereGorski<T>& betaMap, 
+                   double normFreq, 
+                   SphereGorski<T>& intBetaMap)
+{
+  //  int nbpix = intBetaMap.NbPixels();
+  //  SphereGorski<T> newMap(sqrt(nbpix/12));
+  Sph2Sph(betaMap,intBetaMap);
+  for(int i=0; i<intBetaMap.NbPixels(); i++)
+    {
+      RadSpectra* rs =  new PowerLawSpectra
+        (1.,-intBetaMap(i), 0., normFreq);
+      double coeff = rs->filteredIntegratedFlux(sr);
+      intBetaMap(i) = coeff;
+    }
+}
+
+template <class T>
+void addComponentBeta(SphereGorski<T>& finalMap, 
+                      SphereGorski<T>& mapToAdd,SphereGorski<T>& intBetaMap, double K)
+{
+  // finalMap = finalMap + coeff* mapToAdd
+  // coeff    =  convolution of sr and rs 
+  // compute the coefficient corresponding to mapToAdd
+  // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
+  // the spectralResponse
+  // different from addComponentBeta(PixelMap<T>& finalMap, 
+  //    PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
+  //   since it permits to use a intBetaMap with a different number of pixels than
+  //   the other maps
+
+  int nbpix = finalMap.NbPixels();
+  if (nbpix != mapToAdd.NbPixels())   
+    throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
+  double coeff = K;
+
+  if(nbpix != intBetaMap.NbPixels())
+    {
+      cout << "nbpix != intBetaMap.NbPixels()" << endl;
+      SphereGorski<T> bigBetaMap(sqrt(nbpix/12));
+      cout << "new map with size corresponding to mapToAdd" <<  sqrt(nbpix/12) << endl;
+      Sph2Sph(intBetaMap,bigBetaMap);
+      for(int i=0; i<finalMap.NbPixels();i++)
+        {
+          finalMap(i) += coeff*mapToAdd(i)*bigBetaMap(i);
+        }
+    }
+  else
+    {
+      for(int i=0; i<finalMap.NbPixels();i++)
+        {
+          finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
+        }
+    }
+  if (printlev > 1) 
+    {
+      cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= "  << K  << endl;     
+    }
+}
+
+/* Nouvelle-Fonction */
+
+template <class T>
+void Sph2Sph(PixelMap<T>& in, PixelMap<T>& out)
+  // Cette fonction remplit la sphere out a partir de la sphere in
+  // Les spheres peuvent etre de type et de pixelisations differentes
+{
+  int k,kin,kout;
+  double teta,phi;
+  
+  int* cnt = new int[out.NbPixels()+1];
+  for(kout=0; kout<out.NbPixels(); kout++)
+    { cnt[kout] = 0; out(kout) = 0.; }
+  
+  for(kin=0; kin<in.NbPixels(); kin++) {
+    in.PixThetaPhi(kin, teta, phi);
+    kout = out.PixIndexSph(teta, phi);
+    out(kout) += in(kin);
+    
+    cnt[kout] ++;
+    
+  }
+  
+  double moy, sig, dcn;
+  moy = 0.; sig = 0.;
+  for(kout=0; kout<out.NbPixels(); kout++) {
+    dcn = cnt[kout];  moy += dcn;        sig += (dcn*dcn);
+    if (cnt[kout] > 0)  out(kout) /= dcn;
+    else {
+      out.PixThetaPhi(kout, teta, phi);
+      int pixel = in.PixIndexSph(teta,phi);
+      out(kout) = in.PixVal(pixel);
+    }
+  }
+  
+  moy /= out.NbPixels();
+  sig = sig/out.NbPixels() - moy*moy;
+  if (sig >= 0.)        sig = sqrt(sig);
+  
+  delete[] cnt;
+  
+  printf("Sph2Sph_Info NbPix In= %d  Out= %d CntMoy,Sig= %g %g\n",
+         in.NbPixels(), out.NbPixels(), moy, sig);
+  PrtTim("End of  Sph2Sph() ");
+}
+
+
 
 /* Nouvelle-Fonction */