Index: trunk/SophyaProg/PMixer/skymixer.cc
===================================================================
--- trunk/SophyaProg/PMixer/skymixer.cc	(revision 903)
+++ trunk/SophyaProg/PMixer/skymixer.cc	(revision 904)
@@ -59,13 +59,14 @@
 void integratedMap(SpectralResponse&  sr,   
 		   SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
-  //template <class T> 
-  //void integratedMap(SpectralResponse&  sr,   
-  //			      PixelMap<T>& betaMap, 
-  //			      double normFreq, 
-  //			      PixelMap<T>& outMap);
+
+//
 template <class T>
 void addComponentBeta(SphereHEALPix<T>& finalMap, 
 		      SphereHEALPix<T>& mapToAdd,
 		      SphereHEALPix<T>& intBetaMap, double K);
+//
+template <class T> 
+void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap, 
+	       double theta,double phi,double amp,double temp);
 //
 // -----------------------------------------------------------------
@@ -89,56 +90,56 @@
     exit(0);
   }
-
-InitTim();
-
-string msg;
-int rc = 0;
-RadSpectra * es = NULL;
-SpectralResponse * sr = NULL;
-double moy, sig;
-
-DataCards dc;
-so = NULL;
-
- try {
-   string dcard = arg[1];
-   if(printlev > 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);
- }   
-
-
-cout << " skymix/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
-cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev 
-     << "  PrintLev= " << printlev << endl;
-
-// We create an output persist file for writing debug objects
-if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
-
- 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;
-     }
-     if(printlev>0)
-       cout << " Output HealPIx Map read - NbPixels= " << 
-	 outgs.NbPixels() << endl;
+  
+  InitTim();
+  
+  string msg;
+  int rc = 0;
+  RadSpectra * es = NULL;
+  SpectralResponse * sr = NULL;
+  double moy, sig;
+  
+  DataCards dc;
+  so = NULL;
+  
+  try {
+    string dcard = arg[1];
+    if(printlev > 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);
+  }   
+  
+  
+  cout << " skymix/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
+  cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev 
+       << "  PrintLev= " << printlev << endl;
+  
+  // We create an output persist file for writing debug objects
+  if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
+  
+  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;
+      }
+      if(printlev>0)
+	cout << " Output HealPIx Map read - NbPixels= " << 
+	  outgs.NbPixels() << endl;
       if (printlev > 0) {
 	MeanSig(outgs.DataBlock(), moy, sig );
@@ -161,5 +162,4 @@
     string key;
     
-    //    SphereHEALPix<float> ings(hp_nside);  // The input map
     double K = 1.;
     double freqOfMap = 1.;
@@ -170,98 +170,120 @@
       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.);
-      
-      
-      cout << " Reading Input FITS map " << (string)flnm << endl;
-      SphereHEALPix<float> ings(hp_nside);
-      {
-	FITS_SphereHEALPix<float> fiosIn(flnm);
-	ings = (SphereHEALPix<float>)fiosIn;
-      }
-      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;
-      }
-      bool mapDependentOfFreq = false;
-      key = (string)"BETAFITSFILE"+ buff;
+      key = (string)"DIPOLE" + buff;
+      // check for an eventual dipole
       if(dc.HasKey(key)) 
 	{
-	  mapDependentOfFreq = true;
-	}
-
-      // getting Emission spectra  
-      if(!mapDependentOfFreq)
-	{
 	  if (es) { delete es; es = NULL; }
-	  es = getEmissionSpectra(dc, sk);
-	  addComponent(*sr, outgs, ings, *es, K);
+	  double temp = -10.;
+	  double theta=  dc.DParam(key,1,1.);
+	  double phi  =  dc.DParam(key,2,1.);
+	  double amp  =  dc.DParam(key,3,1.);
+	  if(dc.NbParam(key)>3) 
+	    {
+	      temp =  dc.DParam(key,4,1.);
+	    }
+	  cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
+	  addDipole(*sr, outgs,theta,phi,amp,temp);
 	}
       else
 	{
-	  key = (string)"BETAFITSFILE"+ buff;
-	  //SphereHEALPix<float> betaMap;
+	  sprintf(buff, "%d", sk+1);
+	  key = (string)"MAPFITSFILE" + buff;
 	  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;
-	  SphereHEALPix<float> betaMap(hp_nside);
+	  
+	  K = dc.DParam(key, 1, 1.);
+	  
+	  cout << " Reading Input FITS map " << (string)flnm << endl;
+	  SphereHEALPix<float> ings(hp_nside);
 	  {
-	    FITS_SphereHEALPix<float> fiosBM(flnm);
-	    betaMap = (SphereHEALPix<float>)fiosBM;
+	    FITS_SphereHEALPix<float> fiosIn(flnm);
+	    ings = (SphereHEALPix<float>)fiosIn;
 	  }
 	  if (debuglev > 4) {  // Writing the input map to the outppf
-	    FIO_SphereHEALPix<float> fiogs(betaMap);
-	    fiogs.Write(*so, key);
+	    FIO_SphereHEALPix<float> fiog(ings);
+	    fiog.Write(*so, key);
 	  }
-
-	  if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
-	  bool bydefault = true;
-	  if(!bydefault)
-	    addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
+	  if (printlev > 2) {
+	    MeanSig(ings.DataBlock(), moy, sig );
+	    cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
+	  }
+	  bool mapDependentOfFreq = false;
+	  key = (string)"BETAFITSFILE"+ buff;
+	  if(dc.HasKey(key)) 
+	    {
+	      mapDependentOfFreq = true;
+	    }
+	  
+	  // getting Emission spectra  
+	  if(!mapDependentOfFreq)
+	    {
+	      if (es) { delete es; es = NULL; }
+	      es = getEmissionSpectra(dc, sk);
+	      addComponent(*sr, outgs, ings, *es, K);
+	    }
 	  else
 	    {
-	      // integrate the betamap over the SpectralResponse
-	      SphereHEALPix<float> intBetaMap(nSideForInt);
-	      integratedMap(*sr, betaMap, normFreq, intBetaMap);
+	      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;
+	      SphereHEALPix<float> betaMap(hp_nside);
+	      {
+		FITS_SphereHEALPix<float> fiosBM(flnm);
+		betaMap = (SphereHEALPix<float>)fiosBM;
+	      }
 	      if (debuglev > 4) {  // Writing the input map to the outppf
-		FIO_SphereHEALPix<float> fiogs2(intBetaMap);
-		fiogs2.Write(*so, "INTBETAMAP");
+		FIO_SphereHEALPix<float> fiogs(betaMap);
+		fiogs.Write(*so, key);
 	      }
-
-	      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);
+	      
+	      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);
+		  if (debuglev > 4) {  // Writing the input map to the outppf
+		    FIO_SphereHEALPix<float> fiogs2(intBetaMap);
+		    fiogs2.Write(*so, "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;
-      cout << "-------------------------------------------------" << endl;
-      
-      sprintf(buff, "End of Processing Component %d ", sk+1);
-      PrtTim(buff);
+	  
+	  MeanSig(outgs.DataBlock(), moy, sig );
+	  cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
+	  cout << "-------------------------------------------------" << endl;
+	  
+	  sprintf(buff, "End of Processing Component %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);
     } 
- }
- 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 " 
@@ -297,5 +319,5 @@
 
 int rc = 0;
-string key, key2;
+string key, key2,key3;
 
   //  Cheking datacards
@@ -318,6 +340,7 @@
   key = "GAUSSFILTER";
   key2 = "FILTERFITSFILE";
-  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) { 
-    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
+  key3 = "DIPOLE";
+  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) { 
+    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
     rc = 73;  return(rc); 
     }
@@ -339,12 +362,14 @@
   char buff[256];
   bool pb = false;
-  string key3;
   string key4;
   string key5;
+  string key6;
   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 + " " + key2 ;
+    sprintf(buff, "DIPOLE%d", kc+1);
+    key3 = buff;
+    if (dc.NbParam(key) < 1 && dc.NbParam(key3)<1) {   
+      msg = "Missing or invalid card : " + key + " " + key2 + " " + key3;
       pb = true;  break;
     }
@@ -359,6 +384,8 @@
     sprintf(buff, "DIPOLE%d", kc+1);
     key5 = buff;
+    sprintf(buff, "DERIVBB%d", kc+1);
+    key6 = buff;
     if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
-	 && (dc.NbParam(key5)<1))  {   
+	 &&  (dc.NbParam(key6)<1) &&  (dc.NbParam(key5)<3))  {   
       msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
       pb = true;  break;
@@ -455,5 +482,5 @@
   string key = (string)"SPECTRAFITSFILE" + numb;
   string key2 = (string)"BLACKBODY" + numb;
-  string key5 = (string)"DIPOLE" + numb;
+  string key5 = (string)"DERIVBB" + numb;
   string key3 = (string)"POWERLAWSPECTRA" + numb;
   string ppfname = "espectra";
@@ -482,5 +509,5 @@
   }
   else if (dc.HasKey(key5) ) { // Creating Dipole
-    rs = new DerivBlackBody(dc.DParam(key5, 0, 2.726));
+    rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
     ppfname = key5;
   }
@@ -506,4 +533,35 @@
 
 
+template <class T> 
+void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap, 
+	       double theta,double phi,double amp,double temp)
+{
+  DerivBlackBody dbb;
+  if(temp>0) dbb.setTemperature(temp);
+  double coeff = dbb.filteredIntegratedFlux(sr) * amp;
+  UnitVector vd(theta,phi);
+  UnitVector vc(theta,phi);
+  
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      double thetar,phir;
+      finalMap.PixThetaPhi(i,thetar,phir);      
+      vc.SetThetaPhi(thetar,  phir);
+      finalMap(i) += vd.Psc(vc)*coeff;
+    }
+  if (debuglev > 4) {  // Writing the input map to the outppf
+    SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
+    for(int i=0; i<finalMap.NbPixels(); i++)
+      {
+	double thetar,phir;
+	finalMap.PixThetaPhi(i,thetar,phir);      
+	vc.SetThetaPhi(thetar,  phir);
+	ings(i) = vd.Psc(vc);
+      }
+    FIO_SphereHEALPix<float> fiog(ings);
+    fiog.Write(*so, "dipole");
+    cout << "Debug the dipole map....saved in debug file !" << endl;
+  }
+}
 /* Nouvelle-Fonction */
 template <class T> 
