Index: trunk/SophyaProg/PMixer/HowToUse
===================================================================
--- trunk/SophyaProg/PMixer/HowToUse	(revision 761)
+++ trunk/SophyaProg/PMixer/HowToUse	(revision 761)
@@ -0,0 +1,53 @@
+#!/bin/csh
+if ($#argv < 2)  then
+  echo 'Usage: HowToUse NSide_HealPix OutFITSfile'
+  exit 1
+endif
+
+
+echo '    *************************************************************';
+echo '		creating spectral response powerlaw';
+echo '    *************************************************************';
+echo '  ' 
+set cmd = "tgrsr 1 .2,120,50,-0.7,100,600 200 pl.fits"
+echo ' ' 
+echo 'Executing' $cmd
+time $cmd
+
+echo '	*************************************************************';
+echo '		creating gaussian radiation spectra sky';
+echo '	*************************************************************';
+set cmd = "tgsky $1 0 140,30 gauss.fits"
+echo ' ' 
+echo 'Executing' $cmd
+time $cmd
+
+echo '	*************************************************************';
+echo '		creating cosine/sine radiation spectra  sky';
+echo '	*************************************************************';
+set cmd = "tgsky $1 2 12,7,9,10 cossin.fits"
+echo ' ' 
+echo 'Executing' $cmd
+time $cmd
+
+echo '	*************************************************************';
+echo '		creating noise radiation spectra sky';
+echo '	*************************************************************';
+set cmd = "tgsky $1 1 160,44 noise.fits"
+echo ' ' 
+echo 'Executing' $cmd
+time $cmd
+
+echo '	*************************************************************'
+echo '		calling skymixer'
+echo '	*************************************************************'
+set xx = "@SKYMIX  3  $1"
+echo $xx > test.d
+cat howto.d >> test.d
+
+set cmd = "skymixer test.d $2"
+echo ' ' 
+echo 'Executing' $cmd
+time $cmd
+
+
Index: trunk/SophyaProg/PMixer/Makefile
===================================================================
--- trunk/SophyaProg/PMixer/Makefile	(revision 761)
+++ trunk/SophyaProg/PMixer/Makefile	(revision 761)
@@ -0,0 +1,71 @@
+MODULECXXREPNAME :=  PMixer
+include ../Mgr/Makefile.h
+ 
+ifdef NOSHLIB
+LIBF =  $(LIB)libSysTools.a $(LIB)libNTools.a $(LIB)libSamba.a $(LIB)libSkyMap.a $(LIB)libSkyT.a $(LIB)libFitsIOServer.a
+LIBS = -L$(LIB) -lSysTools -lNTools -lSamba -lSkyMap -lSkyT -lFitsIOServer  -lm
+ifeq ($(CXX),cxx)
+# librairie des instances de templates cxx
+LIBS := $(LIBS) -lertcxx
+endif
+else
+LIBF =  $(SLB)libsophya.so $(SLB)libextsophya.so
+LIBS = -L$(SLB) -lsophya -lextsophya  -lm
+endif
+LIBSORA := $(LIBS)
+ifdef ORACLE_HOME
+LIBSORA := -L$(LIB) -lDb++ $(LIBS) -lErosDb -L$(ORACLE_HOME)/lib -lclntsh
+endif
+LDLIBS := $(LIBS) $(LDLIBS)
+PROGS = $(EXE)skymixer $(EXE)extractRS $(EXE)tgrsr $(EXE)tgsky $(EXE)easyTest
+ 
+all: $(PROGS)
+
+clean:
+	rm -f $(PROGS)
+	if [ -d $(OBJ)cxxrep_PMixer ] ; then cd $(OBJ)cxxrep_PMixer ; rm -f *.o ; fi
+
+skymixer: $(EXE)skymixer
+	
+$(EXE)skymixer: $(OBJ)skymixer.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+  
+extractRS: $(EXE)extractRS
+	
+$(EXE)extractRS: $(OBJ)extractRS.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+  
+tgrsr: $(EXE)tgrsr
+	
+$(EXE)tgrsr: $(OBJ)tgrsr.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+  
+tgsky: $(EXE)tgsky
+	
+$(EXE)tgsky: $(OBJ)tgsky.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+  
+easyTest: $(EXE)easyTest
+	
+$(EXE)easyTest: $(OBJ)easyTest.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+  
+$(OBJ)easyTest.o: easyTest.cc $(INC)/blackbody.h \
+ $(INC)/radspec.h \
+ $(INC)/machdefs.h \
+ $(INC)/pclassids.h \
+ $(INC)/specresp.h \
+ $(INC)/anydataobj.h \
+ $(INC)/objfio.h \
+ $(INC)/ppersist.h \
+ $(INC)/pexceptions.h \
+ $(INC)/gnumd5.h \
+ $(INC)/convtools.h \
+ $(INC)/radspecvector.h \
+ $(INC)/tvector.h \
+ $(INC)/tmatrix.h \
+ $(INC)/ndatablock.h \
+ $(INC)/specrespvector.h \
+ $(INC)/squarefilt.h \
+ $(INC)/gaussfilt.h \
+ $(INC)/trianglefilt.h
Index: trunk/SophyaProg/PMixer/RSextract.d
===================================================================
--- trunk/SophyaProg/PMixer/RSextract.d	(revision 761)
+++ trunk/SophyaProg/PMixer/RSextract.d	(revision 761)
@@ -0,0 +1,88 @@
+#  EXT_RS  NbComponents  NSide_HealPix
+@EXT_RS  31 
+
+#  RD_EXT_PATH  PathForFITS   
+#  path for input Sky Component FITS files , from which the RadSpec
+#  will be extracted
+@RS_EXT_PATH /tmp_mnt/exp/planck/CartesRT/supertango/templates/
+
+#  ----    Sky map for different frequencies        ---------
+#  ----    From which will be extracted the radspectra   ----
+#  RADSPECMAPi FITSfilename freq(MHz)
+@RADSPECMAP1 	eran102_441_00015000mhz.fits
+@RADSPECMAP2 	eran102_441_00018884mhz.fits 
+@RADSPECMAP3 	eran102_441_00023773mhz.fits 
+@RADSPECMAP4	eran102_441_00029929mhz.fits 
+@RADSPECMAP5 	eran102_441_00037678mhz.fits 
+@RADSPECMAP6 	eran102_441_00047434mhz.fits 
+@RADSPECMAP7 	eran102_441_00059716mhz.fits 
+@RADSPECMAP8  	eran102_441_00075178mhz.fits 
+@RADSPECMAP9 	eran102_441_00094644mhz.fits 
+@RADSPECMAP10 	eran102_441_00119149mhz.fits 
+@RADSPECMAP11 	eran102_441_00150000mhz.fits
+@RADSPECMAP12 	eran102_441_00188839mhz.fits	
+@RADSPECMAP13 	eran102_441_00237734mhz.fits	00237734
+@RADSPECMAP14 	eran102_441_00299289mhz.fits 
+@RADSPECMAP15 	eran102_441_00376783mhz.fits
+@RADSPECMAP16 	eran102_441_00474342mhz.fits
+@RADSPECMAP17 	eran102_441_00597161mhz.fits
+@RADSPECMAP18 	eran102_441_00751781mhz.fits
+@RADSPECMAP19 	eran102_441_00946436mhz.fits
+@RADSPECMAP20 	eran102_441_01191492mhz.fits
+@RADSPECMAP21 	eran102_441_01500000mhz.fits	
+@RADSPECMAP22 	eran102_441_01888388mhz.fits	
+@RADSPECMAP23 	eran102_441_02377340mhz.fits
+@RADSPECMAP24 	eran102_441_02992893mhz.fits	
+@RADSPECMAP25 	eran102_441_03767830mhz.fits	
+@RADSPECMAP26 	eran102_441_04743416mhz.fits
+@RADSPECMAP27 	eran102_441_05971608mhz.fits	
+@RADSPECMAP28 	eran102_441_07517809mhz.fits	
+@RADSPECMAP29 	eran102_441_09464360mhz.fits
+@RADSPECMAP30 	eran102_441_11914924mhz.fits		
+@RADSPECMAP31 	eran102_441_15000000mhz.fits
+
+
+#  ----   Source point map for different frequencies        ---------
+#  ----    From which will be extracted the radspectra   ----
+#  SOURCEPTMAPi FITSfilename freq(MHz)
+@SOURCEPTMAP1 	eran102_441_00015000mhz.s_rare
+@SOURCEPTMAP2 	eran102_441_00018884mhz.s_rare
+@SOURCEPTMAP3 	eran102_441_00023773mhz.s_rare
+@SOURCEPTMAP4	eran102_441_00029929mhz.s_rare
+@SOURCEPTMAP5 	eran102_441_00037678mhz.s_rare
+@SOURCEPTMAP6 	eran102_441_00047434mhz.s_rare 
+@SOURCEPTMAP7 	eran102_441_00059716mhz.s_rare 	
+@SOURCEPTMAP8  	eran102_441_00075178mhz.s_rare 
+@SOURCEPTMAP9 	eran102_441_00094644mhz.s_rare 
+@SOURCEPTMAP10 	eran102_441_00119149mhz.s_rare
+@SOURCEPTMAP11 	eran102_441_00150000mhz.s_rare
+@SOURCEPTMAP12 	eran102_441_00188839mhz.s_rare	
+@SOURCEPTMAP13 	eran102_441_00237734mhz.s_rare
+@SOURCEPTMAP14 	eran102_441_00299289mhz.s_rare 
+@SOURCEPTMAP15 	eran102_441_00376783mhz.s_rare	
+@SOURCEPTMAP16 	eran102_441_00474342mhz.s_rare	
+@SOURCEPTMAP17 	eran102_441_00597161mhz.s_rare	
+@SOURCEPTMAP18 	eran102_441_00751781mhz.s_rare
+@SOURCEPTMAP19 	eran102_441_00946436mhz.s_rare
+@SOURCEPTMAP20 	eran102_441_01191492mhz.s_rare
+@SOURCEPTMAP21 	eran102_441_01500000mhz.s_rare
+@SOURCEPTMAP22 	eran102_441_01888388mhz.s_rare
+@SOURCEPTMAP23 	eran102_441_02377340mhz.s_rare
+@SOURCEPTMAP24 	eran102_441_02992893mhz.s_rare
+@SOURCEPTMAP25 	eran102_441_03767830mhz.s_rare
+@SOURCEPTMAP26 	eran102_441_04743416mhz.s_rare	
+@SOURCEPTMAP27 	eran102_441_05971608mhz.s_rare
+@SOURCEPTMAP28 	eran102_441_07517809mhz.s_rare
+@SOURCEPTMAP29 	eran102_441_09464360mhz.s_rare
+@SOURCEPTMAP30 	eran102_441_11914924mhz.s_rare		
+@SOURCEPTMAP31 	eran102_441_15000000mhz.s_rare
+
+
+# Filter definition now !!
+# define the detector (only one)
+@GAUSSFILTER 143 12 0.84 50. 300.
+
+#  ---- wanna print a lot of stuff..---------
+#  ----  let's try PRINTLEVEL 5          ----
+@PRINTLEVEL  0
+
Index: trunk/SophyaProg/PMixer/easyTest.cc
===================================================================
--- trunk/SophyaProg/PMixer/easyTest.cc	(revision 761)
+++ trunk/SophyaProg/PMixer/easyTest.cc	(revision 761)
@@ -0,0 +1,188 @@
+#include <iostream.h>
+#include "blackbody.h"
+#include "radspecvector.h"
+#include "specrespvector.h"
+#include "tvector.h"
+#include "squarefilt.h"
+#include "gaussfilt.h"
+#include "trianglefilt.h"
+
+
+int
+main() 
+{ 
+//   PInPersist in("out");
+//   SR_ReadSelf(in);
+//   cout << " " << endl;
+
+  cout << "BlackBody!" << endl;
+  BlackBody myBB(2.73);
+  double freqmin =  230.;
+  double freqmax = 500;
+  double logIF = myBB.logIntegratedFlux(freqmin,freqmax);
+  double IF    = myBB.integratedFlux(freqmin,freqmax);
+  cout << "log   " << logIF << endl;
+  cout << "nolog " << IF    << endl;
+  cout << "min et max " << myBB.minFreq() << " :: " << myBB.maxFreq() << endl;
+  cout << "full integration no log  " << myBB.integratedFlux(freqmin,freqmax) << endl;
+  cout << "full integration log     " << myBB.logIntegratedFlux(freqmin,freqmax) << endl;
+  cout << "full integration no log nolimit " << myBB.integratedFlux() << endl;
+  cout << "full integration log  nolimit  " << myBB.logIntegratedFlux() << endl;
+
+  cout << "squareFilter " << endl;
+  SquareFilter mySF(200,600);
+  double logSF = mySF.logIntegratedSpect(freqmin,freqmax);
+  double SF    = mySF.IntegratedSpect(freqmin,freqmax);
+  cout << "log   " << logSF << endl;
+  cout << "nolog " << SF    << endl;
+  cout << "min et max " << mySF.minFreq() << " :: " << mySF.maxFreq() << endl;
+  cout << "full integration no log  " << mySF.IntegratedSpect(200,700) << endl;
+  cout << "full integration log     " << mySF.logIntegratedSpect(200,700) << endl;
+  cout << "full integration no log nolimit  " << mySF.IntegratedSpect() << endl;
+  cout << "full integration log  nolimit   " << mySF.logIntegratedSpect() << endl;
+
+
+  cout << "triangleFilter" << endl;
+  TriangleFilter myTF(200,500,120,40);
+  double logTF = myTF.logIntegratedSpect(10.,freqmax);
+  double TF    = myTF.IntegratedSpect(10.,freqmax);
+  cout << "log   " << logTF << endl;
+  cout << "nolog " << TF    << endl;
+  cout << "min et max " << myTF.minFreq() << " :: " << myTF.maxFreq() << endl;
+  cout << "full integration no log  " << myTF.IntegratedSpect() << endl;
+  cout << "full integration log     " << myTF.logIntegratedSpect() << endl;
+  
+  cout << "....Testing RadSpectraVec..." << endl;
+  int maxVal = 10;
+  Vector vecOfNu(maxVal);
+  Vector vecOfFDeNu(maxVal);
+  for (int i=0; i<maxVal; i++)
+    {
+      vecOfNu(i) = (double)i+230.;
+      vecOfFDeNu(i) = 1.; //*(double)i+1;
+      cout << "entry values  " <<  vecOfNu(i) << " " << vecOfFDeNu(i) << endl;
+    }
+  freqmin = 200.;
+  freqmax = 259.;
+  RadSpectraVec myRSV(vecOfNu, vecOfFDeNu);
+  double RSV    = myRSV.integratedFlux(freqmin,freqmax);
+ 
+  double logRSV = myRSV.logIntegratedFlux(freqmin,freqmax);
+  cout << "nolog int " << RSV    << endl; 
+  cout << "log int   " << logRSV << endl;
+  cout << "min et max " << myRSV.minFreq() << " :: " << myRSV.maxFreq() << endl;
+  cout << "full integration no log  " << myRSV.integratedFlux() << endl;
+  cout << "full integration log     " << myRSV.logIntegratedFlux() << endl;
+
+  cout << "....Testing SpecRespVec..." << endl;
+  int SRVmaxVal = 100;
+  Vector vecOfNuRSV(SRVmaxVal);
+  Vector vecOfFDeNuRSV(SRVmaxVal);
+  for (int i=0; i<SRVmaxVal; i++)
+    {
+      vecOfNuRSV(i) = (double)i+130.;
+      vecOfFDeNuRSV(i) = 10*(double)i+1;
+      //      cout << "entry values  " <<  vecOfNuRSV(i) << " " << vecOfFDeNuRSV(i) << endl;
+    }
+  freqmin = vecOfNuRSV(0);
+  freqmax = vecOfNuRSV(SRVmaxVal-1);
+  SpecRespVec mySRV(vecOfNuRSV, vecOfFDeNuRSV);
+
+// Test ecriture en PPF 
+  {
+  cout << " Writing mySRV (RadSpecVec to PPF " << endl;
+  POutPersist s("srvec.ppf");
+  s << mySRV;
+  }
+  {
+  cout << " Reading  mySRV from PPF " << endl;
+  PInPersist s("srvec.ppf");
+  SpecRespVec srv;
+  s >> srv;
+  cout << srv;
+  }
+
+  double SRV    = mySRV.IntegratedSpect(freqmin,freqmax);
+  double logSRV = mySRV.logIntegratedSpect(freqmin,freqmax);
+  cout << "nolog " << SRV    << endl; 
+  cout << "log   " << logSRV << endl;
+  cout << "min et max " << mySRV.minFreq() << " :: " << mySRV.maxFreq() << endl;
+  cout << "full integration no log " << mySRV.IntegratedSpect() << endl;
+  cout << "full integration log    " << mySRV.logIntegratedSpect() << endl;
+  
+  cout << "----------------------------------------------------" << endl;
+  cout << "...Testing filter and integrate !" << endl;
+
+  cout << " BBSF nolog " << myBB.filteredIntegratedFlux(mySF,200.,400.) << endl;
+  cout << " BBSF   log " << myBB.filteredLogIntFlux(mySF,200.,400.) << endl;
+  cout << "full integration no log  " << myBB.integratedFlux(200.,400.) << endl;
+  cout << "full integration log     " << myBB.logIntegratedFlux(200.,400.) << endl;
+  cout << "SF" << endl;
+  cout << "full integration no log  " << mySF.IntegratedSpect(200.,400.) << endl;
+  cout << "full integration log     " << mySF.logIntegratedSpect(200.,400.) << endl;
+
+  cout << " BBTF nolog " << myBB.filteredIntegratedFlux(myTF,200.,400.) << endl;
+  cout << " BBTF   log " << myBB.filteredLogIntFlux(myTF,200.,400.) << endl;
+  cout << "TF" << endl;
+  cout << "full integration no log  " << myTF.IntegratedSpect(200.,400.) << endl;
+  cout << "full integration log     " << myTF.logIntegratedSpect(200.,400.) << endl;
+
+  /*
+  cout << " SZSF nolog " << mySZ.filteredIntegratedFlux(mySF,30.,100.) << endl;
+  cout << " SZSF  log " << mySZ.filteredLogIntFlux(mySF,30.,100.) << endl;
+  cout << "full integration no log  " << mySZ.integratedFlux(30.,100.) << endl;
+  cout << "full integration log     " << mySZ.logIntegratedFlux(30.,100.) << endl;
+
+  cout << " SZTF nolog " << mySZ.filteredIntegratedFlux(myTF,30.,100.) << endl;
+  cout << " SZTF   log " << mySZ.filteredLogIntFlux(myTF,30.,100.) << endl;
+			   */
+  cout << "********" << endl;
+  double min, max;
+  min=232;
+  max=235;
+
+  cout << " RSVSF nolog " << myRSV.filteredIntegratedFlux(mySF,min,max) << endl;
+  cout << " RSVSF   log " << myRSV.filteredLogIntFlux(mySF,min,max) << endl;
+  cout << "full integration no log  " << myRSV.integratedFlux(230.,500.) << endl;
+  cout << "full integration log     " << myRSV.logIntegratedFlux(230.,500.) << endl;
+  cout << "full integration no log min max " << myRSV.integratedFlux(min,max) << endl;
+  cout << "full integration log  min max   " << myRSV.logIntegratedFlux(min,max) << endl;
+
+  cout << " RSVTF nolog " << myRSV.filteredIntegratedFlux(myTF,230.,500.) << endl;
+  cout << " RSVTF   log " << myRSV.filteredLogIntFlux(myTF,230.,500.) << endl;
+  
+//   cout << "New TESTS USING IOBJS!" << endl;
+//   POutPersist out("out");
+//   TriangleFilter* myTF2 = new TriangleFilter(10,200,20,40);
+//   SR_WriteSelf(myTF2,out);
+
+//   POutPersist outV("outV");
+//   cout << "check before anything" << vecOfNu.NElts()<< endl;
+//   SpecRespVec* mySRV2 = new SpecRespVec(vecOfNu, vecOfFDeNu);
+//   SR_WriteSelf(mySRV2,outV);
+  
+//   Vector myVec(10);
+//   myVec.ReSize(10);
+//   cout << myVec.NElts();
+//   return 1;
+
+  cout << "read and write things in file...ppersist!" << endl;
+  TriangleFilter myTF2(10,200,20,40);
+//   PInPersist inTF("outTF");
+//   myTF2.ReadSelf(inTF);
+
+  {  // Test IO PPF GaussianFilter
+  cout << " Writing  GaussianFilter gf(276., 39.6, 0.794, 27., 690.) to PPF " << endl;
+  GaussianFilter gf(276., 39.6, 0.794, 27., 690.);
+  POutPersist s("fgauss.ppf");
+  s << gf ;
+  }
+  {
+  cout << " Reading  GaussianFilter from PPF " << endl;
+  PInPersist s("fgauss.ppf");
+  GaussianFilter gf2;
+  s >> gf2;
+  cout << gf2 << endl;
+  }
+
+} 
Index: trunk/SophyaProg/PMixer/extractRS.cc
===================================================================
--- trunk/SophyaProg/PMixer/extractRS.cc	(revision 761)
+++ trunk/SophyaProg/PMixer/extractRS.cc	(revision 761)
@@ -0,0 +1,508 @@
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <fstream.h>
+#include <iostream.h>
+#include <math.h>
+
+#include <string>
+#include <vector>
+
+#include "tmatrix.h"
+#include "tvector.h"
+#include "vector3d.h"
+
+#include "timing.h"
+#include "sambainit.h"
+#include "pexceptions.h"
+#include "datacards.h"
+#include "fitsioserver.h"
+
+#include "radspecvector.h"
+#include "blackbody.h"
+#include "nupower.h"
+
+#include "squarefilt.h"
+#include "trianglefilt.h"
+#include "specrespvector.h"
+#include "gaussfilt.h"
+
+// -----------------------------------------------------------------
+// ------------- Function declaration ------------------------------
+int CheckCards(DataCards & dc, string & msg);
+char * BuildFITSFileName(string const & fname);
+SpectralResponse * getSpectralResponse(DataCards & dc);
+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);
+float ComputeFrequency(DataCards &,string&);
+
+// -----------------------------------------------------------------
+
+//  ----- Global (static) variables ------------
+static bool rdmap = false;    // true -> Read map first 
+static char mapPath[256];     // Path for input maps
+static int  nskycomp = 0;     // Number of sky components
+static int  debuglev = 0;     // Debug Level
+static int  printlev = 0;     // Print Level
+static POutPersist * so = NULL;  // Debug PPFOut file
+
+// -------------------------------------------------------------------------
+//                             main program 
+// -------------------------------------------------------------------------
+int main(int narg, char * arg[])
+{
+  if ((narg < 4) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
+    cout << "  Usage: extractRS parameterFile outputfitsnameformaps  outputfitsnameforsource [outppfname]" << endl;
+    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];
+    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 << " !!!! extractRS/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
+    return(90);
+  }   
+  
+  
+  cout << " extractRS/Info : NComp = " <<  nskycomp << 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("extrRSdbg.ppf");
+  
+  
+  TMatrix<float> outgs;
+  bool okout = false;
+  
+  try {
+    
+    FitsIoServer fios; // Our FITS IO Server
+    char * flnm, buff[90];
+    string key;
+    string key2;
+   
+    double K = 1.;
+    
+    // Loop over sky frequencies components
+    int maxOfSkyComp = nskycomp;
+    int nbOfPix;
+    int sk;
+    Vector freq(maxOfSkyComp);
+    TMatrix<float> spectrum;
+    TMatrix<float> fdenu; //par freq et par pixel
+    int nb_source;
+    TMatrix<float> rareSource; //par freq et par numero de source
+	
+    for(sk = 0; sk<maxOfSkyComp; sk++) 
+      {
+	
+	TMatrix<float> ings;	
+	// Opening the file containing the map
+	
+	if (printlev > 2) cout << " Processing map's frequency component No " << sk+1 << endl;
+	sprintf(buff, "%d", sk+1);
+	key = (string)"RADSPECMAP" + buff;
+	
+	flnm = BuildFITSFileName(dc.SParam(key,0));
+	if (printlev > 2) cout << " Reading Input FITS map " << (string)flnm << endl;
+	
+	fios.load(ings, flnm) ;
+	if (printlev > 2) 
+	  {
+	    MeanSig(ings.DataBlock(), moy, sig );
+	    cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
+	  }
+
+	// Computing the frequency according to the name of the file
+	freq(sk) =  ComputeFrequency(dc,key);
+	if (printlev > 2)  cout << "filling spectrum for freq = " << freq(sk) << endl;
+	// Opening the file containing the point source
+	sprintf(buff, "%d", sk+1);
+	key2 = (string)"SOURCEPTMAP" + buff;
+	flnm = BuildFITSFileName(dc.SParam(key2,0));
+	const char* nameOfFile = flnm ;
+	if (printlev > 2) cout << " Reading Input FITS map " << nameOfFile << endl;
+	ifstream from(nameOfFile);
+	char ch;
+	string dum;
+	from >> nb_source >> dum;
+	if(sk==0) rareSource.ReSize(maxOfSkyComp,nb_source);
+	for (int ii=0; ii< nb_source; ii++)
+	  {
+	    double value = 0;
+	    from >> value;
+	    rareSource(sk,ii) = value*1.e-26; //for Jansky->W
+	  }
+	from.close();
+
+	nbOfPix = ings.NRows()*ings.NCols();
+	spectrum.ReSize(ings.NRows(),ings.NCols());
+	if(sk==0) fdenu.ReSize(maxOfSkyComp,nbOfPix);
+	int xy=0;
+	for(int x=0;x<ings.NRows(); x++)
+	  {
+	    for(int y=0;x<ings.NCols(); x++)
+	      {
+		spectrum(x,y) = ings(x,y)*1.e-26; //for Jansky->W
+		fdenu(sk,xy) = spectrum(x,y);  
+		xy += 1;
+	      }
+	  }
+	if (printlev > 2) 
+	  {
+	    MeanSig(spectrum.DataBlock(), moy, sig );
+	    cout << " MeanSig for spectrum  map - Mean= " << moy << " Sigma= " << sig << endl;
+	  }
+	
+	
+	K = dc.DParam(key, 1, 1.);
+	if (debuglev > 4) {  // Writing the input map to the outppf
+	  FIO_TMatrix<float> fiog;
+	  fiog.Write(*so, key);
+	}
+	
+	if(printlev>2) 
+	  {	
+	    sprintf(buff, "End of Proc. Comp. %d ", sk+1);
+	    cout << " --------------------------------------- \n" << endl;    
+	    PrtTim(buff);
+	  }
+      }   // End of sky component loop 
+    
+    // ---------------------------------------------------------------------
+    //      Integration in the detector band of the maps
+    // ---------------------------------------------------------------------
+    SpectralResponse* sr = NULL;
+    sr = getSpectralResponse(dc);
+  
+    if (printlev > 2)  cout  << "SpectralResponse!" << *sr << endl;
+    TVector<float> outCoeff(nbOfPix);
+
+    for(int pix=1;pix<nbOfPix;pix++)
+      {
+	Vector smallFDeNu(maxOfSkyComp);
+	for(sk = 0; sk<maxOfSkyComp; sk++) 
+	  {
+	    smallFDeNu(sk) =  fdenu(sk,pix);
+	  }
+	
+	RadSpectraVec radSV(freq,smallFDeNu);
+	if (printlev > 10)  cout << "RadiationSpectra" << radSV << endl;
+
+ 	TVector<float> smallFDeNuT(maxOfSkyComp);
+ 	for(sk = 0; sk<maxOfSkyComp; sk++) 
+ 	  {
+ 	    smallFDeNuT(sk) = smallFDeNu(sk);
+ 	  }
+	outCoeff(pix) = radSV.filteredIntegratedFlux(*sr);
+      }
+
+    FitsIoServer fiosrs;
+    fiosrs.save(outCoeff,arg[2]);
+    cout << "Output Map (TMatrix<float>) written to FITS file "  
+	 << arg[2]<< endl;
+    cout << endl;
+
+    // ---------------------------------------------------------------------
+    //      Integration in the detector band of the sources
+    // ---------------------------------------------------------------------
+    TVector<float> sourceCoeff(nb_source);
+    for(int nsource=0;nsource<nb_source;nsource++)
+      {
+	Vector SourceFDeNu(maxOfSkyComp);
+	for(sk = 0; sk<maxOfSkyComp; sk++) 
+	  {
+	    SourceFDeNu(sk) =  rareSource(sk,nsource);
+	  }
+	RadSpectraVec radSV_Source(freq,SourceFDeNu);
+	sourceCoeff(nsource) = radSV_Source.filteredIntegratedFlux(*sr);
+      }
+    FitsIoServer fiosrs_source;
+    fiosrs_source.save(sourceCoeff,arg[3]);
+    cout << "Output Map (TMatrix<float>) written to FITS file "  
+	 << arg[3]<< endl;
+    cout << endl;
+
+    // ---------------------------------------------------------------------
+    //     Check of write/read on FITS file
+    // ---------------------------------------------------------------------
+    if (printlev > 2)
+    {
+      cout << " Check of write/read on FITS file for skymap!" << endl;
+      TVector<float> outout;
+      FitsIoServer fiotest;
+      fiotest.load(outout, arg[2]);
+      cout << "Coeff(2)" << outout(2) << " ? " << outCoeff(2) << endl;
+      cout << endl;
+    }
+    if (printlev > 2)
+    {
+      cout << " Check of write/read on FITS file for sources!" << endl;
+      FitsIoServer fiotest;
+      TVector<float> outout2;
+      fiotest.load(outout2, arg[3]);
+      cout << "sourceCoeff(2)" << outout2(2) << " ? " << sourceCoeff(2) << endl;
+    }
+
+
+  }   // End of try block
+  
+  
+  catch (PException exc) {
+    msg = exc.Msg();
+    cerr << " !!!! extractRS - Catched PException  - Msg= " << exc.Msg() << endl;
+    rc = 50;
+  }   
+  catch (PThrowable exc) {
+    msg = exc.Msg();
+    cerr << " !!!! extractRS - Catched PThrowable  - Msg= " << exc.Msg() << endl;
+    rc = 50;
+  }   
+  catch (...) {
+    cerr << " Une exception de type inconnue !!! " << endl;
+    exit(99);
+  }
+  
+  // Saving the output map in FITS format 
+  if (okout) {
+    FitsIoServer fios;
+    fios.save(outgs, arg[2]);
+    cout << "Output Map (TMatrix<float>) written to FITS file " 
+	 << (string)(arg[2]) << endl;
+    if(printlev>2) PrtTim("End of WriteFITS ");
+    // Saving the output map in PPF format 
+    if (narg > 3) {
+      POutPersist s(arg[3]); 
+      FIO_TMatrix<float> fiog ;
+      fiog.Write(s);
+      cout << "Output Map (TMatrix<float>) written to POutPersist file "  
+	   << (string)(arg[3]) << endl;
+      if(printlev>2) PrtTim("End of WritePPF ");
+    }
+  }
+  if (so) delete so;  //  Closing the debug ppf file 
+  return(rc);
+}
+
+
+
+/* Nouvelle-Fonction */
+int CheckCards(DataCards & dc, string & msg)
+  //   Function to check datacards
+{
+  rdmap = false;
+  mapPath[0] = '\0';
+  nskycomp = 0;
+  debuglev  = 0;
+  printlev = 0;
+  
+  int rc = 0;
+  string key, key2;
+  
+  //  Checking datacards
+  if (dc.NbParam("EXT_RS") < 1)   {
+     rc = 71; 
+     msg = "Invalid parameters  - Check @EXT_RS card ";
+     return(rc);
+  }
+  int kc;
+  char buff[256];
+  bool pb = false;
+  string key3;
+  int ncomp = 0;
+  ncomp = dc.IParam("EXT_RS", 0);
+  
+  for(kc=0; kc<ncomp; kc++) 
+    {
+      key = "RS_EXT_PATH";
+      if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
+      mapPath[255] = '\0';
+      sprintf(buff, "RADSPECMAP%d", kc+1);
+      key = buff;
+      
+      if (dc.NbParam(key) < 1) {   
+	msg = "Missing or invalid card : " + key;
+	pb = true;  break;
+      }
+    }
+  
+  //  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); 
+  }
+  
+  if (pb)  {
+    rc = 75;
+    return(75);
+  }
+  
+  
+  // Initialiazing parameters
+  rc = 0;
+  msg = "OK";
+  nskycomp = ncomp;
+  
+  // Checking for PATH definition card
+  key = "RS_EXT_PATH";
+  if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
+  mapPath[255] = '\0';
+  key = "DEBUGLEVEL";
+  debuglev =  dc.IParam(key, 0, 0);
+  key = "PRINTLEVEL";
+  printlev =  dc.IParam(key, 0, 0);
+  return(rc);
+}
+
+static char buff_flnm[1024];   // Mal protege !
+/* Nouvelle-Fonction */
+char* BuildFITSFileName(string const & fname)
+{
+if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
+else sprintf(buff_flnm, "%s", fname.c_str());
+return(buff_flnm);
+}
+
+
+/* Nouvelle-Fonction */
+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 */
+template <class T> 
+void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
+                  PixelMap<T>& 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;
+  if (printlev > 1) 
+    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl; 
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      finalMap(i) += coeff * mapToAdd(i);
+    }
+}
+
+
+/* Nouvelle-Fonction */
+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);
+}
+
+
+float ComputeFrequency(DataCards & dc,string& key)
+{
+  ofstream filename("temp_filename");
+  filename << dc.SParam(key,0) << '\n';
+  filename.close();
+  
+  ifstream fromtempfile("temp_filename");
+  char chtemp;
+  string temp_str;
+  int count = 0;
+  do { 
+    fromtempfile.get(chtemp);    
+    if(count> 12 && count < 20)
+      {
+	temp_str+= chtemp; 
+      }
+    count+=1;
+  } while(chtemp!='\n'); 
+  fromtempfile.close();
+  
+  ofstream filename2("temp_filename");
+  filename2 << temp_str << '\n';
+  filename2.close();
+  
+  int frequency;
+  ifstream fromtempfile2("temp_filename");
+  fromtempfile2 >> frequency;
+  fromtempfile2.close();
+  return ((float)frequency/1000.);
+}
Index: trunk/SophyaProg/PMixer/howto.d
===================================================================
--- trunk/SophyaProg/PMixer/howto.d	(revision 761)
+++ trunk/SophyaProg/PMixer/howto.d	(revision 761)
@@ -0,0 +1,47 @@
+#  SKYMIX  NbComponents  NSide_HealPix
+SKYMIX  3  128
+
+#  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 /data/opus/Rz/SOPHYA/
+
+
+#  ----        Sky components  ---------
+#      ----     HealPix maps   ----
+#  MAPFITSFILEi FITSfilename Normalisation
+@MAPFITSFILE1  gauss.fits  1.0 
+
+#      ---   Emission Spectra definition -----
+#   SPECTRAFITSFILEi  FITSfilename  Fmin Fmax 
+@SPECTRAFITSFILE1  pl.fits 100 600  
+
+#   Other possible definition of emission spectra
+#   BLACKBODYi  temperature
+#  The blackbody object returns flux in units of watt/m2/sr/hz
+#   We thus scale the map (x 1.e17) 
+@MAPFITSFILE2  cossin.fits  1.e17 
+@BLACKBODY2 2.726
+
+#   For power-law spectra  f = a ((nu-nu0)/dnu)^b)
+#   POWERLAWSPECTRAi  a nu0 dnu b Fmin Fmax
+@MAPFITSFILE3  noise.fits    0.1 
+@POWERLAWSPECTRA3 1. 150. 50. -0.5 100. 500.
+
+#  Define the Debug level
+@DEBUGLEVEL  0
+#  Define the Print level
+@PRINTLEVEL  5
Index: trunk/SophyaProg/PMixer/makeprog
===================================================================
--- trunk/SophyaProg/PMixer/makeprog	(revision 761)
+++ trunk/SophyaProg/PMixer/makeprog	(revision 761)
@@ -0,0 +1,49 @@
+MODULECXXREPNAME := SkyT
+# MODULEDECCXXFLAGS := -msg_quiet
+include $(DPCDEVREP)/Include/Makefile.h
+
+ifdef NOSHLIB
+LIBF =  $(LIB)libFitsIOServer.a $(LIB)libSkyT.a $(LIB)libSamba.a $(LIB)libNTools.a $(LIB)libSysTools.a 
+LIBS = -L$(LIB) -lSkyT -lSamba -lNTools -lSysTools -lm
+ifeq ($(CXX),cxx)
+# librairie des instances de templates cxx
+LIBS := $(LIBS) -lertcxx
+endif
+else
+LIBF =  $(SLB)/libtools.so 
+LIBS = -L$(SLB) -ltools -lm 
+ifeq ($(MACHEROS),Linux)
+LIBS := $(LIBS) -ldl
+endif
+endif
+
+
+myexe = $(EXE)easyTest $(EXE)skymixer $(EXE)extractRS $(EXE)tgsky $(EXE)tgrsr 
+
+all : $(myexe)
+
+
+$(EXE)easyTest :  $(OBJ)easyTest.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)skymixer :  $(OBJ)skymixer.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)extractRS :  $(OBJ)extractRS.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)tgsky :  $(OBJ)tgsky.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)tgrsr :  $(OBJ)tgrsr.o 
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(OBJ)easyTest.o : easyTest.cc
+
+$(OBJ)tgsky.o : tgsky.cc
+
+$(OBJ)tgrsr.o : tgrsr.cc
+
+$(OBJ)skymixer.o : skymixer.cc
+
+$(OBJ)extractRS.o : extractRS.cc
Index: trunk/SophyaProg/PMixer/programs
===================================================================
--- trunk/SophyaProg/PMixer/programs	(revision 761)
+++ trunk/SophyaProg/PMixer/programs	(revision 761)
@@ -0,0 +1,5 @@
+skymixer skymixer.o
+extractRS extractRS.o
+tgrsr tgrsr.o
+tgsky tgsky.o
+easyTest easyTest.o
Index: trunk/SophyaProg/PMixer/skm.d
===================================================================
--- trunk/SophyaProg/PMixer/skm.d	(revision 761)
+++ trunk/SophyaProg/PMixer/skm.d	(revision 761)
@@ -0,0 +1,43 @@
+#  SKYMIX  NbComponents  NSide_HealPix
+@SKYMIX  2  64
+
+#  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  
+
+#   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  0
+#  Define the Print level
+@PRINTLEVEL  0
Index: trunk/SophyaProg/PMixer/skymixer.cc
===================================================================
--- trunk/SophyaProg/PMixer/skymixer.cc	(revision 761)
+++ trunk/SophyaProg/PMixer/skymixer.cc	(revision 761)
@@ -0,0 +1,478 @@
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+
+#include <iostream.h>
+#include <math.h>
+
+#include <string>
+#include <vector>
+
+#include "timing.h"
+#include "sambainit.h"
+#include "pexceptions.h"
+#include "datacards.h"
+#include "fitsioserver.h"
+
+#include "spheregorski.h"
+
+#include "radspecvector.h"
+#include "blackbody.h"
+#include "nupower.h"
+
+#include "squarefilt.h"
+#include "trianglefilt.h"
+#include "specrespvector.h"
+#include "gaussfilt.h"
+
+// -----------------------------------------------------------------
+// ------------- 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);
+
+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;     // Debug Level
+static int  printlev = 0;     // Print Level
+static POutPersist * so = NULL;  // Debug PPFOut file
+
+// -------------------------------------------------------------------------
+//                             main program 
+// -------------------------------------------------------------------------
+int main(int narg, char * arg[])
+{
+  if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
+    cout << "  Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
+    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];
+  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 << " 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);
+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;
+
+  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;
+    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, 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);
+      }
+    // 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 = 50;
+  }   
+
+// 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);
+}
+
+/* Nouvelle-Fonction */
+int CheckCards(DataCards & dc, string & msg)
+//   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.NbParam("SKYMIX") < 2)   {
+     rc = 71; 
+     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
+  int mg = 32;
+  int ncomp = 0;
+  ncomp = dc.IParam("SKYMIX", 0, 0);
+  mg = dc.IParam("SKYMIX", 1, 32);
+  if (ncomp < 1)  {
+    msg = "Invalid parameters  - Check datacards @SKYMIX ";
+    rc = 74;
+    return(rc);
+  }
+
+  // Checking detection filter specification
+//  Checking input FITS file specifications 
+  int kc;
+  char buff[256];
+  bool pb = false;
+  string key3;
+  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;
+      pb = true;  break;
+    }
+    sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
+    key = buff;
+    sprintf(buff, "BLACKBODY%d", kc+1);
+    key2 = buff;
+    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;
+    }
+
+    }
+
+  if (pb)  {
+    rc = 75;
+    return(75);
+  }
+
+
+// Initialiazing parameters
+  rc = 0;
+  msg = "OK";
+  nskycomp = ncomp;
+  hp_nside = mg;  
+
+// Checking for PATH definition card
+  key = "MAPPATH";
+  if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
+  mapPath[255] = '\0';
+  key = "DEBUGLEVEL";
+  debuglev =  dc.IParam(key, 0, 0);
+  key = "PRINTLEVEL";
+  printlev =  dc.IParam(key, 0, 0);
+  return(rc);
+}
+
+static char buff_flnm[1024];   // Mal protege !
+/* Nouvelle-Fonction */
+char* BuildFITSFileName(string const & fname)
+{
+if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
+else sprintf(buff_flnm, "%s", fname.c_str());
+return(buff_flnm);
+}
+
+/* Nouvelle-Fonction */
+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
+  // 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;
+  if (printlev > 1) 
+    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl; 
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      finalMap(i) += coeff * mapToAdd(i);
+    }
+}
+
+/* Nouvelle-Fonction */
+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 */
+void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
+{
+  char *ntn[2] = {"nu","fnu"};
+  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
+  float xnt[2];
+  double nu;
+  double numin = rs.minFreq();
+  double numax = rs.maxFreq();
+  int nmax = 500;
+  double dnu = (numax-numin)/nmax;
+  for(int k=0; k<nmax; k++) {
+    nu = numin+k*dnu;
+    xnt[0] = nu;
+    xnt[1] = rs.flux(nu);
+    nt.Fill(xnt);
+  }
+  ObjFileIO<NTuple> oiont(nt);
+  oiont.Write(so, name);
+  return;
+}
+
+/* Nouvelle-Fonction */
+void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
+{
+  char *ntn[2] = {"nu","tnu"};
+  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
+  float xnt[2];
+  double nu;
+  double numin = sr.minFreq();
+  double numax = sr.maxFreq();
+  int nmax = 500;
+  double dnu = (numax-numin)/nmax;
+  for(int k=0; k<nmax; k++) {
+    nu = numin+k*dnu;
+    xnt[0] = nu;
+    xnt[1] = sr.transmission(nu);
+    nt.Fill(xnt);
+  }
+  ObjFileIO<NTuple> oiont(nt);
+  oiont.Write(so, name);
+  return;
+}
+
Index: trunk/SophyaProg/PMixer/tgrsr.cc
===================================================================
--- trunk/SophyaProg/PMixer/tgrsr.cc	(revision 761)
+++ trunk/SophyaProg/PMixer/tgrsr.cc	(revision 761)
@@ -0,0 +1,151 @@
+#include "machdefs.h"
+#include <iostream.h>
+#include <typeinfo>
+#include <math.h>
+
+#include "tvector.h"
+#include "cimage.h"
+#include "sambainit.h"
+
+#include "radspecvector.h"
+#include "specrespvector.h"
+#include "nupower.h"
+#include "gaussfilt.h"
+
+#include "fitsioserver.h"
+
+#include "timing.h"
+
+// Test program to generate different type of RadSpectra
+// and SpectralResponse in the form of RadSpectraVec and SpecRespVec
+// in FITS format
+
+// ------------- Main program --------------
+int main(int narg, char* arg[]) 
+{
+  PeidaInit();
+  InitTim();   // Initializing the CPU timer 
+
+  if ((narg < 5) || ((narg>1) && (strcmp(arg[1],"-h") == 0)) )  {
+    cout << " tgrsr : Generation of RadSpectraVec and SpecRespVec in FITS  " << endl;
+    cout << " Usage: tgrsr 0/1 Params NPoints FitsFileName [ImageR4FITSName] [ppfname] " << endl;
+    cout << " 0 -> GaussianFilter  T = A exp(((nu-nu0)/dnu)^2" << endl;
+    cout << " Params : A,Nu0,DNu,MinFreq,MaxFreq " << endl;
+    cout << " 1 -> PowerLawSpectra f = a ((nu-nu0)/dnu)^b " << endl;
+    cout << " Params : A,Nu0,dnu,b,MinFreq,MaxFreq " << endl;
+    exit(0);
+    }
+
+  int typ = 0; 
+  typ = atoi(arg[1]);
+  int k, npt = 100;
+  npt = atoi(arg[3]);
+  if (npt < 5) npt = 5;
+  cout << " ---- tgrsr/Info: Type= " << typ << " NPoints= " << npt << endl;
+  Matrix mtx(2, npt);
+  ImageR4 img(npt, 2);
+
+  char buff[256];
+  if (typ == 0) {  // Generation de SpecRespVec from GaussianFilter
+    double a,nu0,dnu,fmin,fmax;
+    a = 1.0;
+    nu0 = 143.;
+    dnu = 10.;
+    fmin = 10.;
+    fmax = 1000.;
+    sscanf(arg[2],"%lg,%lg,%lg,%lg,%lg",&a,&nu0,&dnu,&fmin,&fmax);
+    cout << " GaussianFilter  T = A exp(((nu-nu0)/dnu)^2" << endl;
+    sprintf(buff,"A=%g Nu0=%g DNu=%g FMin=%g FMax=%g",
+            a,nu0,dnu,fmin,fmax);
+    cout << "Decoded params: " << (string)buff << endl;
+    GaussianFilter sr(nu0, dnu, a, fmin, fmax);
+    cout << sr << endl;
+    cout << " Creating SpecRespVec() ... " << endl;
+    Vector v1(10);
+    Vector v2(10);
+    SpecRespVec srv(v1, v2, fmin, fmax);
+    Vector & nu = srv.getNuVec();
+    Vector & tnu = srv.getTNuVec();
+    nu.ReSize(npt);
+    tnu.ReSize(npt);
+    double freq;
+    double dfreq = (fmax-fmin)/(npt-3);
+    for(k=0; k<npt; k++) {
+      freq = fmin+(k-1)*dfreq;
+      nu(k) = freq;
+      tnu(k) = sr.transmission(freq);
+      //      if ((k>10) && (k<25))  cout << " DBG - K = " << k <<
+      //        " Nu= " << nu(k) << " Tnu= " << tnu(k) << endl; 
+      mtx(0,k) = nu(k);
+      mtx(1,k) = tnu(k);
+      img(k,0) = nu(k);
+      img(k,1) = tnu(k);
+    }
+    if ((typ == 0) && (narg > 6)) { // Writing to ppf  file
+      cout << "writing SpecRespVec(" << npt << ") to PPF file: " 
+           << (string)(arg[6]) << endl;
+      cout << sr ; 
+      cout << srv ; 
+      {
+        POutPersist pos(arg[6]);
+        pos << srv;
+      }
+      PrtTim("End of ImageR4->PPF ");
+    }
+
+  }
+  else { // Generation de RadSpectraVec from PowerLawSpectra
+    double a,nu0,dnu,b,fmin,fmax;
+    a = 1.0;
+    nu0 = 143.;
+    b = 0.;
+    dnu = 10.;
+    fmin = 10.;
+    fmax = 1000.;
+    sscanf(arg[2],"%lg,%lg,%lg,%lg,%lg,%lg",&a,&nu0,&dnu,&b,&fmin,&fmax);
+    cout << " PowerLawSpectra f = a ((nu-nu0)/dnu)^b " << endl;
+    sprintf(buff,"A=%g Nu0=%g DNu=%g B=%g FMin=%g FMax=%g",
+            a,nu0,dnu,b,fmin,fmax);
+    cout << "Decoded params: " << (string)buff << endl;
+    PowerLawSpectra rs(a,b, nu0, dnu, fmin, fmax);
+    cout << rs << endl;
+    cout << " Creating RadSpectraVec() ... " << endl;
+    Vector v1(10);
+    Vector v2(10);
+    RadSpectraVec rsv(v1, v2, fmin, fmax);
+    Vector & nu = rsv.getNuVec();
+    Vector & fnu = rsv.getFNuVec();
+    nu.ReSize(npt);
+    fnu.ReSize(npt);
+    double freq;
+    double dfreq = (fmax-fmin)/(npt-3);
+    for(k=0; k<npt; k++) {
+      freq = fmin+(k-1)*dfreq;
+      nu(k) = freq;
+      fnu(k) = rs.flux(freq);
+      //      if ((k>10) && (k<25))  cout << " DBG2 - K = " << k <<
+      //        " Nu= " << nu(k) << " Fnu= " << fnu(k) << endl; 
+      mtx(0,k) = nu(k);
+      mtx(1,k) = fnu(k);
+      img(k,0) = nu(k);
+      img(k,1) = fnu(k);
+    }
+  }
+
+  PrtTim("End of filling ");
+
+  FitsIoServer fios;
+  cout << "\n Writing Matrix NRows= " << mtx.NRows() << "  NCols= " 
+       << mtx.NCols() << " to FITS file: " << (string)(arg[4]) << endl;
+  fios.save(mtx, arg[4]);
+  PrtTim("End of Matrix->FITS ");
+
+  if (narg > 5) { // Writing ImageR4 fo FITS file
+    cout << "wrting ImageR4(" << npt << ", 2) to FITS file: " 
+         << (string)(arg[5]) << endl;
+    fios.save(img, arg[5]);
+    PrtTim("End of ImageR4->FITS ");
+  }
+  cout <<  " ============ End of tgrsr program ======== " << endl;
+  return 0;
+}
Index: trunk/SophyaProg/PMixer/tgsky.cc
===================================================================
--- trunk/SophyaProg/PMixer/tgsky.cc	(revision 761)
+++ trunk/SophyaProg/PMixer/tgsky.cc	(revision 761)
@@ -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);
+}
Index: trunk/SophyaProg/PMixer/usedlibs
===================================================================
--- trunk/SophyaProg/PMixer/usedlibs	(revision 761)
+++ trunk/SophyaProg/PMixer/usedlibs	(revision 761)
@@ -0,0 +1,6 @@
+SysTools
+NTools
+Samba
+SkyMap
+SkyT
+FitsIOServer
Index: trunk/SophyaProg/PMixer/usedslibs
===================================================================
--- trunk/SophyaProg/PMixer/usedslibs	(revision 761)
+++ trunk/SophyaProg/PMixer/usedslibs	(revision 761)
@@ -0,0 +1,2 @@
+sophya
+extsophya
