Index: trunk/SophyaLib/SkyT/Makefile
===================================================================
--- trunk/SophyaLib/SkyT/Makefile	(revision 608)
+++ trunk/SophyaLib/SkyT/Makefile	(revision 610)
@@ -24,5 +24,5 @@
 myobjs = $(OBJ)radspec.o $(OBJ)radspecvector.o $(OBJ)nupower.o $(OBJ)blackbody.o $(OBJ)specresp.o $(OBJ)specrespvector.o $(OBJ)squarefilt.o  $(OBJ)trianglefilt.o $(OBJ)gaussfilt.o $(OBJ)convtools.o 
 
-myexe = $(EXE)easyTest $(EXE)skymixer
+myexe = $(EXE)easyTest $(EXE)skymixer $(EXE)tgsky $(EXE)tgrsr 
 
 all : $(myexe)
@@ -36,5 +36,11 @@
 	$(LINK.cc)  $^ $(LIBS) -o $@
 
+$(EXE)tgsky :  $(OBJ)tgsky.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)tgrsr :  $(OBJ)tgrsr.o $(myobjs)
+	$(LINK.cc)  $^ $(LIBS) -o $@
 
 
 
+
Index: trunk/SophyaLib/SkyT/gaussfilt.cc
===================================================================
--- trunk/SophyaLib/SkyT/gaussfilt.cc	(revision 608)
+++ trunk/SophyaLib/SkyT/gaussfilt.cc	(revision 610)
@@ -1,5 +1,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: gaussfilt.cc,v 1.2 1999-11-20 21:00:48 ansari Exp $
+//      $Id: gaussfilt.cc,v 1.3 1999-11-21 23:25:45 ansari Exp $
 //
 // Description:
@@ -50,5 +50,5 @@
   else { 
     double tmp = (nu-_nu0)/_s;
-    return(_a * exp(tmp*tmp));
+    return(_a * exp(-tmp*tmp));
   }
 }
@@ -70,5 +70,5 @@
 { 
   os << "GaussianFilter::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
-  os << " T = A * Exp(((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= " 
+  os << " T = A * Exp(-((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= " 
      << _s << "  A= " << _a << endl;  
   os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;
Index: trunk/SophyaLib/SkyT/gaussfilt.h
===================================================================
--- trunk/SophyaLib/SkyT/gaussfilt.h	(revision 608)
+++ trunk/SophyaLib/SkyT/gaussfilt.h	(revision 610)
@@ -2,5 +2,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: gaussfilt.h,v 1.2 1999-11-20 21:00:49 ansari Exp $
+//      $Id: gaussfilt.h,v 1.3 1999-11-21 23:25:45 ansari Exp $
 //
 // Description:
@@ -27,5 +27,5 @@
 //              ---------------------
 
-// Spectral response in the form  A Exp(((nu-nu0)/s)^2) 
+// Spectral response in the form  A Exp(-((nu-nu0)/s)^2) 
 
 class GaussianFilter:public SpectralResponse
Index: trunk/SophyaLib/SkyT/nupower.cc
===================================================================
--- trunk/SophyaLib/SkyT/nupower.cc	(revision 608)
+++ trunk/SophyaLib/SkyT/nupower.cc	(revision 610)
@@ -1,5 +1,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: nupower.cc,v 1.2 1999-11-20 21:00:49 ansari Exp $
+//      $Id: nupower.cc,v 1.3 1999-11-21 23:25:45 ansari Exp $
 //
 // Description:
@@ -24,5 +24,5 @@
 // Constructor --
 //----------------
-PowerLawSpectra::PowerLawSpectra(double a, double b, double nu0, double numin, double numax)
+PowerLawSpectra::PowerLawSpectra(double a, double b, double nu0, double dnu, double numin, double numax)
   : RadSpectra(numin, numax)
 {
@@ -30,4 +30,5 @@
   _b = b;
   _nu0 = nu0;
+  _dnu = (dnu > 1.e-19) ? dnu : 1.;
 }
 
@@ -40,6 +41,8 @@
 PowerLawSpectra::flux(double nu) const
 {
-  if ((nu < _numin) || (nu > _numax)) return(0.);
-  else return(  _a * pow((nu-_nu0), _b) );
+  if ((nu< _numin) || (nu > _numax)) return(0.);
+  double tmp = (nu-_nu0)/_dnu;
+  if (tmp <= 0.)  return(0.); 
+  else return(  _a * pow(tmp, _b) );
 }
 
@@ -47,6 +50,6 @@
 PowerLawSpectra::Print(ostream& os) const
 {
-  os << "PowerLawSpectra::Print f = a (nu-nu0)^b " << _a << "(nu-" << _nu0 
-     << ") ^ " << _b << endl;
+  os << "PowerLawSpectra::Print f = a ((nu-nu0)/dnu)^b " << _a << "((nu-" << _nu0 
+     << ") / " << _dnu << ") ^ " << _b << endl;
   os << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
   os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
Index: trunk/SophyaLib/SkyT/nupower.h
===================================================================
--- trunk/SophyaLib/SkyT/nupower.h	(revision 608)
+++ trunk/SophyaLib/SkyT/nupower.h	(revision 610)
@@ -2,5 +2,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: nupower.h,v 1.2 1999-11-20 21:00:49 ansari Exp $
+//      $Id: nupower.h,v 1.3 1999-11-21 23:25:46 ansari Exp $
 //
 // Description:
@@ -17,10 +17,10 @@
 #include "convtools.h"
 
-// Power law spectra  f = a (nu-nu0)^b   
+// Power law spectra  f = a ((nu-nu0)/dnu)^b   
 class PowerLawSpectra : public RadSpectra
 {
 public:   //Constructor
   
-  PowerLawSpectra(double a, double b, double nu0, double numin=0., double numax=9.e49);
+  PowerLawSpectra(double a, double b, double nu0, double dnu, double numin=0., double numax=9.e49);
 
   virtual ~PowerLawSpectra();
@@ -32,5 +32,5 @@
 
 protected:
-  double _a, _b, _nu0;
+  double _a, _b, _nu0, _dnu;
 };
 
Index: trunk/SophyaLib/SkyT/radspecvector.cc
===================================================================
--- trunk/SophyaLib/SkyT/radspecvector.cc	(revision 608)
+++ trunk/SophyaLib/SkyT/radspecvector.cc	(revision 610)
@@ -1,5 +1,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: radspecvector.cc,v 1.2 1999-11-20 21:00:50 ansari Exp $
+//      $Id: radspecvector.cc,v 1.3 1999-11-21 23:25:46 ansari Exp $
 //
 // Description:
@@ -57,13 +57,12 @@
 {
   if ( (nu < _numin) || (nu > _numax) ) return(0.);
-  double value = -99;
+  double value = 0.;
   int sizeVecOfNu = _vecOfNu.NElts();
-  double minVal = _vecOfNu(0);
-  if(nu <= minVal) return _vecOfFDeNu(0);
+  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
   if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
   
   for (int i=1; i<sizeVecOfNu; i++)
     {
-      if(nu>= minVal && nu < _vecOfNu(i))
+      if(nu < _vecOfNu(i))
 	{
 	  double up = _vecOfFDeNu(i) ;
@@ -75,10 +74,6 @@
 	  return value;
 	}
-      else 
-	{
-	  minVal = _vecOfNu(i);
-	}	
     }
-  return value;
+  return _vecOfFDeNu(sizeVecOfNu-1);
 }
 
Index: trunk/SophyaLib/SkyT/skm.d
===================================================================
--- trunk/SophyaLib/SkyT/skm.d	(revision 608)
+++ trunk/SophyaLib/SkyT/skm.d	(revision 610)
@@ -1,18 +1,43 @@
-#  MAPPATH  PathForFITS files
-@MAPPATH /exp/eros/Reza/CartesRT
-#  SKYMIX  NbComponents  M_Gorski
+#  SKYMIX  NbComponents  NSide_HealPix
 @SKYMIX  2  64
-#  FILTER  Nu0 Sigma_Nu Tmax NuMin NuMax
-@FILTER 143 12 0.84 50. 300.
 
-#  Sky components 
-#  MAPFITSFILEi FITSfilename [Normalisation]
-@MAPFITSFILE1  map_essai.fits 
+#  READMAP  FitsName
+#  To Read the already prepared (mixed) Output MAP from FITS file
+#  FitsName should be the complete path to the FITS file
+READMAP xxx.fits
+
+#  Defining the Detection filter pass - band
+#  GAUSSFILTER or FILTERFITSFILE  card
+#  GAUSSFILTER -> Gaussian filter
+#  FILTERFITSFILE -> Filter (nu, T(nu)) from FITS file
+#  GAUSSFILTER  Nu0 Sigma_Nu Tmax NuMin NuMax
+#  FILTERFITSFILE FileName NuMin NuMax
+#       FileName should be a complete path
+@GAUSSFILTER 143 12 0.84 50. 300.
+
+#  MAPPATH  PathForFITS   
+#  path for input Sky Component FITS files , and EmissionSpectra files 
+MAPPATH /exp/eros/Reza/CartesRT
+
+
+#  ----        Sky components  ---------
+#      ----     HealPix maps   ----
+#  MAPFITSFILEi FITSfilename Normalisation
+@MAPFITSFILE1  map_essai.fits  1.0 
 @MAPFITSFILE2    comp2.fits    0.25 
 
+#      ---   Emission Spectra definition -----
 #   SPECTRAFITSFILEi  FITSfilename  Fmin Fmax 
 @SPECTRAFITSFILE1  spec1.fits  
-@SPECTRAFUNC2  3. 4. 5.
+
+#   Other possible definition of emission spectra
+#   BLACKBODYi  temperature
+BLACKBODY1 2.726
+#   For power-law spectra  f = a ((nu-nu0)/dnu)^b)
+#   POWERLAWSPECTRAi  a nu0 dnu b Fmin Fmax
+POWERLAWSPECTRA1 1. 150. 50. -0.5 100. 500.
 
 #  Define the Debug level
-@DEBUGLEVEL  5
+@DEBUGLEVEL  0
+#  Define the Print level
+@PRINTLEVEL  0
Index: trunk/SophyaLib/SkyT/skymixer.cc
===================================================================
--- trunk/SophyaLib/SkyT/skymixer.cc	(revision 608)
+++ trunk/SophyaLib/SkyT/skymixer.cc	(revision 610)
@@ -26,20 +26,28 @@
 #include "gaussfilt.h"
 
-// ------ Function declaration 
-void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
-                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K=1.);
-void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
-                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K=1.);
-
+// -----------------------------------------------------------------
+// ------------- Function declaration ------------------------------
 int CheckCards(DataCards & dc, string & msg);
 char * BuildFITSFileName(string const & fname);
+SpectralResponse * getSpectralResponse(DataCards & dc);
+RadSpectra * getEmissionSpectra(DataCards & dc, int nc);
 void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
 void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
 
-//  ----- Variable globale ------------
+template <class T>
+void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
+                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K=1.);
+template <class T>
+void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
+// -----------------------------------------------------------------
+
+//  ----- Global (static) variables ------------
+static bool rdmap = false;    // true -> Read map first 
 static char mapPath[256];     // Path for input maps
 static int  hp_nside = 32;    // HealPix NSide
 static int  nskycomp = 0;     // Number of sky components
-static int  debuglev = 0;     // debuglevel
+static int  debuglev = 0;     // Debug Level
+static int  printlev = 0;     // Print Level
+static POutPersist * so = NULL;  // Debug PPFOut file
 
 // -------------------------------------------------------------------------
@@ -57,76 +65,124 @@
 string msg;
 int rc = 0;
-POutPersist * so = NULL;
+RadSpectra * es = NULL;
+SpectralResponse * sr = NULL;
+double moy, sig;
+
+DataCards dc;
+so = NULL;
 
 try {
   string dcard = arg[1];
   cout << " Decoding parameters from file " << dcard << endl;
-  DataCards dc(dcard);
+  dc.ReadFile(dcard);
 
   rc = CheckCards(dc, msg);
-  if (rc) goto Fin;
-
-  cout << " skymix/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
-  cout << "  ... MapPath = " << (string)mapPath << "  DbgLev= " << debuglev << endl;
+  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 << " SphereGorski_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");
-
-  SphereGorski<float> outgs(hp_nside);
-  cout << " Output Gorski Map  created - NbPixels= " << outgs.NbPixels() << endl;
-  outgs.SetPixels(0.);
+if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
+
+SphereGorski<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", 1).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;
+  }
+  else {
+    cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
+    outgs.SetPixels(0.);
+  }
 
   // Decoding detection pass-band filter 
-  double nu0 = dc.DParam("FILTER", 0, 10.);
-  double s = dc.DParam("FILTER", 1, 1.);
-  double a = dc.DParam("FILTER", 2, 1.);
-  double numin = dc.DParam("FILTER", 3, 0.1);
-  double numax = dc.DParam("FILTER", 4, 9999);
-  GaussianFilter filt(nu0, s, a, numin, numax);
-  cout << " Filter decoded - Created " << endl;
-  cout << filt << endl;
-
-// FOR debug
-  if (debuglev > 0)  SpectralResponse2Nt(filt, *so, "filter");
-
+  sr = getSpectralResponse(dc);
   PrtTim(" After FilterCreation ");
 
-  SphereGorski<float> * ings = NULL;  // Our input map  
   FitsIoServer fios; // Our FITS IO Server
-  char * flnm, buff[64];
+  char * flnm, buff[90];
   string key;
 
+  SphereGorski<float> ings(hp_nside);  // The input map
+  double K = 1.;
   // Loop over sky component 
   int sk;
   for(sk = 0; sk<nskycomp; sk++) {
     cout << " Processing sky component No " << sk+1 << endl;
-    if (ings) { delete ings;  ings = NULL; }
-    ings = new SphereGorski<float>(hp_nside);
-    sprintf(buff, "MAPFITSFILE%d", sk+1);
-    key = buff;
+    sprintf(buff, "%d", sk+1);
+    key = (string)"MAPFITSFILE" + buff;
     flnm = BuildFITSFileName(dc.SParam(key, 0));
     cout << " Reading Input FITS map " << (string)flnm << endl;
-    fios.load(*ings, flnm, 1);
-
+    fios.load(ings, flnm, 2);
+    K = dc.DParam(key, 1, 1.);
     if (debuglev > 4) {  // Writing tne input map to the outppf
       FIO_SphereGorski<float> fiog(ings);
       fiog.Write(*so, key);
-      }  
-    }
-  }   // End of try block
+      }
+    // 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);
+    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 = 90;
+  rc = 50;
   }   
 
-
-Fin:
-if (so) delete so;  //  Closing the debug ppf file 
-if (rc == 0)  return(0);
-cerr << " Error condition -> Rc= " << rc << endl;
-cerr << " Msg= " << msg << endl;
-return(rc);
+// Saving the output map in FITS format 
+if (okout) {
+  FitsIoServer fios;
+  fios.save(outgs, arg[2]);
+  cout << "Output Map (SphereGorski<float>) written to FITS file " 
+       << (string)(arg[2]) << endl;
+  PrtTim("End of WriteFITS ");
+// Saving the output map in PPF format 
+  if (narg > 3) {
+   POutPersist s(arg[3]); 
+   FIO_SphereGorski<float> fiog(&outgs) ;
+   fiog.Write(s);
+   cout << "Output Map (SphereGorski<float>) written to POutPersist file "  
+        << (string)(arg[3]) << endl;
+   PrtTim("End of WritePPF ");
+   }
+  }
+  if (so) delete so;  //  Closing the debug ppf file 
+  return(rc);
 }
 
@@ -135,16 +191,37 @@
 //   Function to check datacards
 {
+rdmap = false;
 mapPath[0] = '\0';
 hp_nside = 32;
 nskycomp = 0;
 debuglev  = 0;
+printlev = 0;
 
 int rc = 0;
+string key, key2;
+
   //  Cheking datacards
-  if ( (!dc.HasKey("SKYMIX")) || (!dc.HasKey("FILTER")) ) {
+  if (dc.NbParam("SKYMIX") < 2)   {
      rc = 71; 
-     msg = "Invalid parameters  - NO @SKYMIX or @FILTER card ";
+     msg = "Invalid parameters  - Check @SKYMIX card ";
      return(rc);
   }
+  key = "READMAP";
+  if (dc.HasKey(key)) {
+    if (dc.NbParam(key) < 1) {
+     rc = 72; 
+     msg = "Invalid parameters  - Check @READMAP card ";
+     return(rc);
+    }
+   else rdmap = true;
+  }
+
+//  Checking detection filter specification 
+  key = "GAUSSFILTER";
+  key2 = "FILTERFITSFILE";
+  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) { 
+    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
+    rc = 73;  return(rc); 
+    }
 
   // Decoding number of component and pixelisation parameter
@@ -155,13 +232,14 @@
   if (ncomp < 1)  {
     msg = "Invalid parameters  - Check datacards @SKYMIX ";
-    rc = 72;
+    rc = 74;
     return(rc);
   }
 
+  // Checking detection filter specification
 //  Checking input FITS file specifications 
   int kc;
-  string key, key2;
   char buff[256];
   bool pb = false;
+  string key3;
   for(kc=0; kc<ncomp; kc++) {
     sprintf(buff, "MAPFITSFILE%d", kc+1);
@@ -173,8 +251,10 @@
     sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
     key = buff;
-    sprintf(buff, "SPECTRAFUNC%d", kc+1);
+    sprintf(buff, "BLACKBODY%d", kc+1);
     key2 = buff;
-    if ( (dc.NbParam(key) < 1) && (dc.NbParam(key2) < 2) )  {   
-      msg = "Missing card or invalid parameters : " + key + " or " + key2;
+    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;
       pb = true;  break;
     }
@@ -183,14 +263,8 @@
 
   if (pb)  {
-    rc = 72;
-    return(72);
-  }
-
-//  Checking detection filter specification 
-  key = "FILTER";
-  if (dc.NbParam(key) < 3) { 
-    msg = "Missing card or invalid parameters : " + key;
-    rc = 73;  return(rc); 
-    }
+    rc = 75;
+    return(75);
+  }
+
 
 // Initialiazing parameters
@@ -206,4 +280,6 @@
   key = "DEBUGLEVEL";
   debuglev =  dc.IParam(key, 0, 0);
+  key = "PRINTLEVEL";
+  printlev =  dc.IParam(key, 0, 0);
   return(rc);
 }
@@ -219,7 +295,106 @@
 
 /* Nouvelle-Fonction */
-// template <class T> 
-void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
-                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K)
+SpectralResponse * getSpectralResponse(DataCards & dc)
+{
+  SpectralResponse * filt = NULL;
+
+  string key = "FILTERFITSFILE";
+  string key2 = "GAUSSFILTER";
+  string ppfname = "filter";
+
+  if (dc.HasKey(key) ) {      // Reading FITS filter file
+    FitsIoServer fios;
+    char ifnm[256];
+    strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
+    ifnm[255] = '\0';
+    Matrix mtx(2,10);
+    fios.load(mtx, ifnm);
+    double numin = dc.DParam(key, 2, 1.);
+    double numax = dc.DParam(key, 3, 9999.);
+    Vector nu(mtx.NCols());
+    Vector fnu(mtx.NCols());
+    for(int k=0; k<mtx.NCols(); k++) {
+      nu(k) = mtx(0, k);
+      fnu(k) = mtx(1, k);
+    }
+  filt = new SpecRespVec(nu, fnu, numin, numax);
+  ppfname = key;
+  }
+  else if (dc.HasKey(key2) ) {   // creating GaussianFilter 
+    double nu0 = dc.DParam(key2, 0, 10.);
+    double s = dc.DParam(key2, 1, 1.);
+    double a = dc.DParam(key2, 2, 1.);
+    double numin = dc.DParam(key2, 3, 0.1);
+    double numax = dc.DParam(key2, 4, 9999);
+    filt = new GaussianFilter(nu0, s, a, numin, numax);
+    ppfname = key2;
+    }
+  if (filt == NULL) throw ParmError("datacard error ! No detection filter");
+  cout << " Filter decoded - Created " << endl;
+  cout << *filt << endl;
+
+// for debug
+  if (debuglev > 1)  SpectralResponse2Nt(*filt, *so, ppfname);
+  return(filt);
+}
+
+/* Nouvelle-Fonction */
+RadSpectra * getEmissionSpectra(DataCards & dc, int nc)
+{
+  char numb[16];
+  sprintf(numb, "%d", nc+1);
+
+  string key = (string)"SPECTRAFITSFILE" + numb;
+  string key2 = (string)"BLACKBODY" + numb;
+  string key3 = (string)"POWERLAWSPECTRA" + numb;
+  string ppfname = "espectra";
+
+  RadSpectra * rs = NULL;
+  if (dc.HasKey(key) ) {    // Reading emission spectra from file
+    char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
+    cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
+    FitsIoServer fios;
+    Matrix mtx(2,10);
+    fios.load(mtx, ifnm);
+    double numin = dc.DParam(key, 2, 1.);
+    double numax = dc.DParam(key, 3, 9999.);
+    Vector nu(mtx.NCols());
+    Vector tnu(mtx.NCols());
+    for(int k=0; k<mtx.NCols(); k++) {
+      nu(k) = mtx(0, k);
+      tnu(k) = mtx(1, k);
+    }
+  rs = new RadSpectraVec(nu, tnu, numin, numax);
+  ppfname = key;
+  }
+  else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
+    rs = new BlackBody(dc.DParam(key2, 0, 2.726));
+    ppfname = key2;
+  }
+  else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
+    double a = dc.DParam(key3, 0, 1.);
+    double nu0 = dc.DParam(key3, 1, 100.);
+    double dnu = dc.DParam(key3, 2, 10.);
+    double b = dc.DParam(key3, 3, 0.);
+    double numin = dc.DParam(key3, 4, 0.1);
+    double numax = dc.DParam(key3, 5, 9999);
+    rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
+    ppfname = key3;
+  }
+  if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
+  cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
+  cout << *rs << endl;
+// for debug
+  if (debuglev > 2)  RadSpec2Nt(*rs, *so, ppfname);
+  return(rs);
+}
+
+
+
+
+/* Nouvelle-Fonction */
+template <class T> 
+void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
+                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
 {
   // finalMap = finalMap + coeff* mapToAdd
@@ -229,4 +404,6 @@
     throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
   double coeff = rs.filteredIntegratedFlux(sr) * K;
+  if (printlev > 1) 
+    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl; 
   for(int i=0; i<finalMap.NbPixels(); i++)
     {
@@ -236,19 +413,19 @@
 
 /* Nouvelle-Fonction */
-void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
-                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K)
-{
-  // finalMap = finalMap + coeff* mapToAdd
-  // coeff    =  convolution of sr and rs 
-  // compute the coefficient corresponding to mapToAdd
-  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
-    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
-  double coeff = rs.filteredIntegratedFlux(sr) * K;
-  for(int i=0; i<finalMap.NbPixels(); i++)
-    {
-      finalMap(i) += coeff * mapToAdd(i);
-    }
-}
-
+template <class T>
+void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig)
+
+{
+  gmoy=0.;
+  gsig = 0.; 
+  double valok;
+  for(int k=0; k<(int)dbl.Size(); k++) {
+    valok = dbl(k);
+    gmoy += valok;  gsig += valok*valok; 
+  }
+  gmoy /= (double)dbl.Size();
+  gsig = gsig/(double)dbl.Size() - gmoy*gmoy;
+  if (gsig >= 0.) gsig = sqrt(gsig);
+}
 
 /* Nouvelle-Fonction */
Index: trunk/SophyaLib/SkyT/specrespvector.cc
===================================================================
--- trunk/SophyaLib/SkyT/specrespvector.cc	(revision 608)
+++ trunk/SophyaLib/SkyT/specrespvector.cc	(revision 610)
@@ -1,5 +1,5 @@
 //--------------------------------------------------------------------------
 // File and Version Information:
-//      $Id: specrespvector.cc,v 1.2 1999-11-20 21:00:53 ansari Exp $
+//      $Id: specrespvector.cc,v 1.3 1999-11-21 23:25:47 ansari Exp $
 //
 // Description:
@@ -57,13 +57,12 @@
 {
   if ( (nu < _numin) || (nu > _numax) ) return(0.);
-  double value = -99;
+  double value = 0.;
   int sizeVecOfNu = _vecOfNu.NElts();
-  double minVal = _vecOfNu(0);
-  if(nu <= minVal) return _vecOfFDeNu(0);
+  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
   if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
 
   for (int i=1; i<sizeVecOfNu; i++)
     {
-      if(nu>= minVal && nu < _vecOfNu(i))
+      if(nu < _vecOfNu(i))
 	{
 	  double up = _vecOfFDeNu(i) ;
@@ -75,8 +74,4 @@
 	  return value;
 	}
-      else 
-	{
-	  minVal = _vecOfNu(i);
-	}	
     }
   return value;
Index: trunk/SophyaLib/SkyT/tgsky.cc
===================================================================
--- trunk/SophyaLib/SkyT/tgsky.cc	(revision 610)
+++ trunk/SophyaLib/SkyT/tgsky.cc	(revision 610)
@@ -0,0 +1,116 @@
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+
+#include "sambainit.h"
+#include "spheregorski.h"
+#include "fitsioserver.h"
+#include "nbrandom.h"
+#include "bruit.h"
+#include "timing.h"
+
+// Test program to generate different type of random skies
+template <class T>
+void  MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig);
+
+// ------------- Main program --------------
+int main(int narg, char* arg[]) 
+{
+  double teta,phi;
+  double gmoy=0., gsig=0.;
+
+  PeidaInit();
+  InitTim();   // Initializing the CPU timer 
+  if ((narg < 5) || ((narg>1) && (strcmp(arg[1],"-h") == 0)) )  {
+    cout << " tgsky : Generation of random skies  " << endl;
+    cout << " Usage: tgsky NSide TypSky Params FitsFileName [PPFName] " << endl;
+    cout << "   - TypSky = 0  Params= m,sig  , Gaussian(m,sig) " << endl;
+    cout << "   - TypSky = 1  Params= m,sig  , OOFNoise(sig) + m" << endl;
+    cout << "   - TypSky = 2  Params= K,a,b,m  K*cos(a*teta)*sin(b*phi)+m" << endl;   
+    exit(0);
+    }
+
+  int  ns = atoi(arg[1]);
+  SphereGorski<float> sph(ns);
+
+  cout << "Filling SphereGorski<float> NSide= " << ns 
+       << "  NPixels= " << sph.NbPixels() << endl;
+
+  int typ = atoi(arg[2]);
+  double m,s,a,b,K;
+  m = 0.;
+  s = 1.;
+  a = b = 1.;
+  K = 1.;
+  if ((typ < 0) || (typ > 2)) typ = 0;
+  if (typ < 2)  { 
+    sscanf(arg[3],"%lg,%lg",&m,&s);
+    cout << "  TypSky= " << typ << " m= " << m << " sig= " << s << endl;
+    if (typ == 1) cout << "  ---> OOFNoise() " << endl;
+    else cout << "  ---> NoiseGenerator() " << endl;
+  }
+  else { 
+    sscanf(arg[3],"%lg,%lg,%lg,%lg",&K,&a,&b,&m);
+    cout << "  TypSky= " << typ << " K= " << m << " a= " << a << " b=" 
+         << b << " m= " << m << endl;
+    cout << "    ---> K * cos(a*teta) * sin(b*teta) + m" << endl;
+    }
+
+  PrtTim("End of Init ");
+
+  if (typ == 2) {
+    for (int j=0;j<sph.NbPixels();j++)   {
+      sph.PixThetaPhi(j,teta,phi);
+      sph(j)= K * cos(a*teta)*sin(b*phi) + m;
+    }
+  }
+  else {
+    NoiseGenerator * ng;
+    if (typ == 0) ng = new NoiseGenerator(s);
+    else ng = new OOFNoise(s);    
+    for (int j=0;j<sph.NbPixels();j++)   sph(j) = ng->Noise()+m;
+    delete ng;
+  }
+  PrtTim("End of Fill ");
+
+  // Computing mean and sigma on the sphere 
+  MeanSig(sph, gmoy, gsig);
+  cout << "SphereGorski<float> Mean= " << gmoy << "  Sigma = " << gsig << endl;
+  PrtTim("End of Mean-Sig ");
+
+  if (narg > 5) {
+   POutPersist s(arg[5]); 
+   FIO_SphereGorski<float> fiog(&sph) ;
+   fiog.Write(s);
+   cout << "SphereGorski<float> written to POutPersist file"  
+        << (string)(arg[5]) << endl;
+   PrtTim("End of WritePPF ");
+   }
+
+
+   
+  FitsIoServer fios;
+  fios.save(sph, arg[4]);
+  cout << "SphereGorski<float> written to FITS file " << (string)(arg[4]) << endl;
+  PrtTim("End of WriteFITS ");
+
+  cout <<  " ============ End of tgsky program ======== " << endl;
+  return 0;
+}
+
+/* Nouvelle-Methode */
+template <class T>
+void  MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig)
+
+{
+  gmoy=0.;
+  gsig = 0.; 
+  double valok;
+  for(int k=0; k<map.NbPixels(); k++) {
+    valok = map(k);
+    gmoy += valok;  gsig += valok*valok; 
+  }
+  gmoy /= (double)map.NbPixels();
+  gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
+  if (gsig >= 0.) gsig = sqrt(gsig);
+}
