Index: trunk/Cosmo/SimLSS/Makefile
===================================================================
--- trunk/Cosmo/SimLSS/Makefile	(revision 3194)
+++ trunk/Cosmo/SimLSS/Makefile	(revision 3196)
@@ -19,5 +19,6 @@
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
        $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
-       $(EXE)cmvtpoisson $(EXE)cmvconcherr $(EXE)cmvtinterp
+       $(EXE)cmvtpoisson $(EXE)cmvconcherr $(EXE)cmvtinterp \
+       $(EXE)cmvtstagn
 #$(EXE)cmvtluc
  
@@ -26,9 +27,10 @@
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o $(OBJ)cmvtinterp.o \
-          $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o
+          $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o \
+          $(OBJ)cmvtstagn.o
  
 LIBROBJ = \
           $(OBJ)cosmocalc.o $(OBJ)pkspectrum.o $(OBJ)schechter.o \
-          $(OBJ)planckspectra.o $(OBJ)integfunc.o $(OBJ)geneutils.o \
+          $(OBJ)planckspectra.o $(OBJ)geneutils.o $(OBJ)agnjackson.o \
           $(OBJ)genefluct3d.o
 
@@ -56,12 +58,12 @@
 $(OBJ)schechter.o: schechter.cc schechter.h constcosmo.h
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ schechter.cc
-$(OBJ)integfunc.o: integfunc.cc integfunc.h
-	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ integfunc.cc
 $(OBJ)planckspectra.o: planckspectra.cc planckspectra.h constcosmo.h
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ planckspectra.cc
-$(OBJ)pkspectrum.o: pkspectrum.cc pkspectrum.h constcosmo.h integfunc.h
+$(OBJ)pkspectrum.o: pkspectrum.cc pkspectrum.h constcosmo.h
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ pkspectrum.cc
 $(OBJ)geneutils.o: geneutils.cc geneutils.h
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ geneutils.cc
+$(OBJ)agnjackson.o: agnjackson.cc agnjackson.h
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ agnjackson.cc
 $(OBJ)genefluct3d.o: genefluct3d.cc genefluct3d.h
 	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ genefluct3d.cc
@@ -173,4 +175,12 @@
 
 ##############################################################################
+cmvtstagn: $(EXE)cmvtstagn
+	echo $@ " done"
+$(EXE)cmvtstagn: $(OBJ)cmvtstagn.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstagn.o $(MYLIB) 
+$(OBJ)cmvtstagn.o: cmvtstagn.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtstagn.cc
+
+##############################################################################
 cmvtluc: $(EXE)cmvtluc
 	echo $@ " done"
Index: trunk/Cosmo/SimLSS/agnjackson.cc
===================================================================
--- trunk/Cosmo/SimLSS/agnjackson.cc	(revision 3196)
+++ trunk/Cosmo/SimLSS/agnjackson.cc	(revision 3196)
@@ -0,0 +1,96 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <vector>
+
+#include "pexceptions.h"
+
+#include "histos.h"
+#include "srandgen.h"
+
+#include "geneutils.h"
+#include "agnjackson.h"
+
+
+AGNJackson::AGNJackson(void)
+  : nobjang_(0.) , fluxang_(0.) , dndls_(NULL), tirls_(NULL)
+{
+
+  //--- Les valeurs lues sur la courbe de Jackson (from JIM)
+  const int nval = 10;
+  // log10(S en Jy)
+  double x[nval] = {-7.,-6.,-5.,-4.,-3.,-2.,-1.,0.,1.,2.};
+  // S^2.5*dN/dS (en Jy^1.5 /sr)
+  double y[nval] = {0.0035,0.14,0.95,2.5,7.,50.,260.,300.,150.,100.};
+
+  if(sizeof(x)!=sizeof(y) || sizeof(y)!=8*nval) {
+    cout<<"AGNJackson::AGNJackson_Error: incompatible x,y,nval sizes"<<endl;
+    throw SzMismatchError("AGNJackson::AGNJackson_Error: incompatible x,y,nval sizes");
+  }
+
+  xjack_.resize(0); yjack_.resize(0);
+  for(int i=0;i<nval;i++)
+    {xjack_.push_back(x[i]); yjack_.push_back(y[i]);}
+
+  // On cree l'histo: dndls = dN/dlog10(S) en fct de  log10(S en Jy)
+  int nbin = nval*20;
+  double xmin = x[0], xmax=x[nval-1];
+  dndls_ = new Histo(xmin,xmax,nbin);
+  dndls_->Zero();
+
+  // On remplit les histos
+  // yjack_: S^2.5 * dN/dS (en Jy^1.5 /sr) en fct de log10(S en Jy)
+  // dndls: dN/dlog10(S)  (en 1/log10(Jy)/sr) en fct de log10(Jy en Jy)
+  //   avec dN/dlog10(S) = dN/dS * dS/dlog10(S)
+  //                     = (S^2.5*dN/dS)*S^(-2.5) * (ln(10)*S)
+  // On va interpoler en log(dN/dlog10(S)) versus log10(S)
+  vector<double> Y;
+  for(int i=0;i<nval;i++) {
+    double v = yjack_[i]*pow(10.,-1.5*xjack_[i])*M_LN10;
+    Y.push_back(log(v));
+  }
+  for(int i=0;i<nbin;i++) {
+    double lf = dndls_->BinCenter(i);
+    // Interpolation lineaire (parabolic not OK!)
+    double v = InterpTab(lf,xjack_,Y,1);
+    (*dndls_)(i) = exp(v);
+  }
+
+  // On remplit la fonction de tirage
+  tirls_ = new FunRan(*dndls_,true);
+
+  // On calcul le nombre moyen d'AGN et le flux moyen par sr
+  // ATTENTION: dn/dlog10(S) en fct de log10(s)
+  nobjang_ = fluxang_ = 0.;
+  for(int i=0;i<nbin;i++) {
+    double lf = dndls_->BinCenter(i);
+    nobjang_ += (*dndls_)(i);
+    fluxang_ += (*dndls_)(i)*pow(10.,lf);
+  }
+  nobjang_ *= dndls_->BinWidth();
+  fluxang_ *= dndls_->BinWidth();
+
+}
+
+AGNJackson::~AGNJackson(void)
+{
+  if(dndls_!=NULL) delete dndls_;
+  if(tirls_!=NULL) delete tirls_;
+}
+
+void AGNJackson::Print(void)
+{
+  cout<<"AGNJackson::Print: nobj="<<nobjang_
+      <<" /sr, flux="<<fluxang_<<" Jy/sr"<<endl;
+  dndls_->Show();
+}
+
+void AGNJackson::OrigJack(vector<double>& xjack,vector<double>& yjack)
+{
+  xjack = xjack_;
+  yjack = yjack_;
+}
Index: trunk/Cosmo/SimLSS/agnjackson.h
===================================================================
--- trunk/Cosmo/SimLSS/agnjackson.h	(revision 3196)
+++ trunk/Cosmo/SimLSS/agnjackson.h	(revision 3196)
@@ -0,0 +1,39 @@
+#ifndef AGNJACKSON_SEEN
+#define AGNJACKSON_SEEN
+
+#include "machdefs.h"
+#include <math.h>
+#include "genericfunc.h"
+#include "histos.h"
+#include "perandom.h"
+
+namespace SOPHYA {
+
+//----------------------------------------------------
+// C.A.Jackson (New Astronomy Review 40 (2004) 1187-1193)
+class AGNJackson {
+public:
+  AGNJackson(void);
+  virtual ~AGNJackson(void);
+
+  Histo& dNdlS(void) {return *dndls_;}
+  FunRan& TirL10S(void) {return *tirls_;}
+  double NObjAng(void) {return nobjang_;}
+  double FluxAng(void) {return fluxang_;}
+  double L10FluxJY(void) {return tirls_->Random();}
+  double FluxJY(void) {return pow(10.,L10FluxJY());}
+  void Print(void);
+
+  void OrigJack(vector<double>& xjack,vector<double>& yjack);
+
+protected:
+  vector<double> xjack_,yjack_;
+  double nobjang_,fluxang_;
+  Histo *dndls_;
+  FunRan *tirls_;
+
+};
+
+}
+
+#endif
Index: trunk/Cosmo/SimLSS/cmvdefsurv.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3196)
@@ -11,10 +11,9 @@
 #include "cosmocalc.h"
 #include "geneutils.h"
-#include "integfunc.h"
 #include "schechter.h"
 #include "planckspectra.h"
 
 /* --- Check Peterson at al. astro-ph/0606104 v1
-cmvdefsurv -z 0.0025 -x 1 -U 0.75,0.3,0.7,-1,1 -V 300 -O 400000,6000 -A 75 -M 6.156e9 -F 3 -2 1.5
+cmvdefsurv -z 0.0025 -x 1 -U 0.75,0.3,0.7,-1,1 -V 300 -O 400000,6000 -N 75 -M 6.156e9 -F 3 -2 1.5
  --- */
 
@@ -35,5 +34,5 @@
       <<" -L lobewidth : taille du lobe d\'observation (FWHM) en arcmin (def= 1\')"<<endl
       <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
-      <<" -A Tsys : temperature du system (K)"<<endl
+      <<" -N Tsys : temperature du system (K)"<<endl
       <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
       <<"                   (indnu==0 no evolution with freq.)"<<endl
@@ -43,4 +42,5 @@
       <<" -U h100,om0,ol0,w0,or0,flat : cosmology"<<endl
       <<" -2 : two polarisations measured"<<endl
+      <<" -A <log10(S_agn)> : moyenne du flux AGN en Jy dans le pixel"<<endl
       <<" redshift : redshift moyen du survey"<<endl
       <<endl;
@@ -59,5 +59,4 @@
  double alpha = -1.37;
  cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
-
 
  // --- Arguments
@@ -76,8 +75,9 @@
  double vrot = 300.; // largeur en vitesse (km/s) pour elargissement doppler
  double facpolar = 0.5; // si on ne mesure les 2 polars -> 1.0
+ double lflux_agn = -3.;
 
  // --- Decodage arguments
  char c;
-  while((c = getopt(narg,arg,"hP2x:y:z:A:S:O:M:F:V:U:L:")) != -1) {
+  while((c = getopt(narg,arg,"hP2x:y:z:N:S:O:M:F:V:U:L:A:")) != -1) {
   switch (c) {
   case 'P' :
@@ -99,5 +99,5 @@
     sscanf(optarg,"%lf",&lobewidth);
     break;
-  case 'A' :
+  case 'N' :
     sscanf(optarg,"%lf",&Tsys);
     break;
@@ -119,4 +119,7 @@
   case '2' :
     facpolar = 1.0;
+    break;
+  case 'A' :
+    sscanf(optarg,"%lf",&lflux_agn);
     break;
   case 'h' :
@@ -281,7 +284,8 @@
  double lobecyl = sqrt(8.)*slobe; // diametre du lobe cylindrique equiv en arcmin
  double lobearea = M_PI*lobecyl*lobecyl/4.;  // en arcmin^2 (hypothese lobe gaussien)
+ double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
+ if(lobewidth<=0.) nlobes = 1.;
  cout<<"\nBeam FWHM = "<<lobewidth<<"\' -> sigma = "<<slobe<<"\' -> "
      <<" Dcyl = "<<lobecyl<<"\' -> area = "<<lobearea<<" arcmin^2"<<endl;
- double nlobes = rad2min(adtx)*rad2min(adty)/lobearea;
  cout<<"Number of beams in one transversal pixel = "<<nlobes<<endl;
 
@@ -290,9 +294,9 @@
  cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
 
- double fhi = hifactor*FluxHI(mhiref,dlum);
+ double fhi = hifactor*Msol2FluxHI(mhiref,dlum);
  cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
      <<"  Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
-     <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
-     <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
+     <<"      in ["<<hifactor*Msol2FluxHI(mhiref,dlumlim[0])
+     <<","<<hifactor*Msol2FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
  double sfhi = fhi / (dnuhiz*1e9) / Jansky2Watt_cst;
  cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
@@ -342,4 +346,13 @@
  cout<<"                       flux density = "<<scmb<<" Jy"<<endl;
 
+ // AGN
+ double flux_agn = pow(10.,lflux_agn);
+ cout<<"AGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="<<flux_agn<<" Jy"<<endl;
+ double flux_agn_pix = flux_agn*(dnuhiz*1e9);
+ double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);
+ double lmass_agn_pix = log10(mass_agn_pix);
+ cout<<"...pixel: f="<<flux_agn_pix<<" 10^-26 W/m^2"
+     <<" -> "<<mass_agn_pix<<" Msol -> log10 = "<<lmass_agn_pix<<endl;
+
  // --- Noise analysis
  cout<<"\n--- Noise analysis"<<endl;
Index: trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3196)
@@ -15,5 +15,4 @@
 #include "schechter.h"
 #include "geneutils.h"
-#include "integfunc.h"
 #include "genefluct3d.h"
 
@@ -32,4 +31,5 @@
      <<" -2 : compute 2D spectrum"<<endl
      <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
+     <<" -A <log10(S_agn en Msol)>,sigma : distribution des AGN par pixel"<<endl
      <<" -W : write cube in FITS format"<<endl
      <<" -P : write cube in PPF format"<<endl
@@ -77,4 +77,8 @@
  double snoise= 0.;   // en equivalent MSol
 
+ // *** AGN
+ bool do_agn = false;
+ double lmsol_agn=-99., lsigma_agn=0.;   // en equivalent MSol
+
  // *** type de generation
  bool computefourier0=false;
@@ -91,5 +95,5 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
   switch (c) {
   case 'a' :
@@ -122,4 +126,8 @@
   case 'M' :
     sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
+    break;
+  case 'A' :
+    do_agn = true;
+    sscanf(optarg,"%lf,%lf",&lmsol_agn,&lsigma_agn);
     break;
   case 'V' :
@@ -153,4 +161,5 @@
      <<", schnpt="<<schnpt<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
+ if(do_agn) cout<<"AGN: <log10(Msol)>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
 
  //-----------------------------------------------------------------
@@ -476,4 +485,13 @@
  }
 
+ if(do_agn) {
+   cout<<"\n--- Add AGN: <mass>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
+   fluct3d.AddAGN(lmsol_agn,lsigma_agn);
+   nm = fluct3d.MeanSigma2(rm,rs2);
+   cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+       <<rs2<<" -> "<<sqrt(rs2)<<endl;
+   PrtTim(">>>> End Add AGN");
+ }
+
  if(snoise>0.) {
    cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
Index: trunk/Cosmo/SimLSS/cmvtintfun.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3196)
@@ -10,5 +10,5 @@
 #include "ntuple.h"
 
-#include "integfunc.h"
+#include "geneutils.h"
 
 // cmvtintfun [n] [xmin,xmax,npt]
Index: trunk/Cosmo/SimLSS/cmvtstagn.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstagn.cc	(revision 3196)
+++ trunk/Cosmo/SimLSS/cmvtstagn.cc	(revision 3196)
@@ -0,0 +1,198 @@
+// Distribution des flux des AGN dans un pixel de taille donne
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include "timing.h"
+#include "ntuple.h"
+
+#include "geneutils.h"
+#include "agnjackson.h"
+
+void usage(void);
+void usage(void)
+{
+cout<<"cmvtstagn "<<endl
+ <<" -a angsol: taille pixel (en sr si >0 ou arcmin^2 si <0)"<<endl
+ <<" -x nbin,xmin,xmax: def histo des log10(S Jy)/pixel"<<endl
+ <<"              si nbin<0 on a de <S>+xmin , <S>-xmax"<<endl
+ <<" -n npix: nombre de tirages de pixel"<<endl;
+}
+
+
+int main(int narg,char *arg[])
+{
+ double angsol = 1.;  // sr
+ double lxmin=-3.,lxmax=4.; int nbinlx = -100;
+ int npix = 100;
+
+ // --- Decodage arguments
+ char c;
+ while((c = getopt(narg,arg,"ha:x:n:")) != -1) {
+  switch (c) {
+  case 'a' :
+    sscanf(optarg,"%lf",&angsol);
+    break;
+  case 'x' :
+    sscanf(optarg,"%d,%lf,%lf",&nbinlx,&lxmin,&lxmax);
+    break;
+  case 'n' :
+    sscanf(optarg,"%d",&npix);
+    break;
+  case 'h' :
+  default :
+    usage();
+    return -1;
+  }
+ }
+
+ // --- Class de distribution d'AGN
+ double sr2amin = AngSol(M_PI/180./60./2.,M_PI/180./60./2.);
+ if(angsol<0) angsol *= -sr2amin;
+ cout<<"angsol = "<<angsol<<" sr -> "<<angsol/sr2amin<<" arcmin^2"<<endl;
+ AGNJackson agn;
+ agn.Print();
+
+ cout<<endl;
+ vector<double> xjack, yjack; agn.OrigJack(xjack,yjack);
+ Histo hdndls0 = agn.dNdlS();
+ Histo hdndls(hdndls0); hdndls *= angsol;
+ Histo htirls = agn.TirL10S();
+ double nobjang = agn.NObjAng()*angsol;
+ double fluxang = agn.FluxAng()*angsol;
+ cout<<"Nombre moyen d\'objets dans notre pixel: "<<nobjang<<endl;
+ cout<<"Flux moyen dans notre pixel: "<<fluxang<<" Jy"<<endl;
+
+ // --- Distribution sur angsol
+ // Histo du nombre d'objets / angsol
+ cout<<endl;
+ double mu = nobjang; if(mu<=0.) mu=1.;
+ double xmin = mu-5.*sqrt(mu); if(xmin<0.) xmin=0.; xmin=floor(xmin);
+ double xmax = mu+7.*sqrt(mu); xmax=floor(xmax+1.);
+ long nb = long(xmax-xmin+1.01); while(nb>250) nb /= 2;
+ // Histo du nombre d'objets/pixel
+ Histo hnobj(xmin,xmax,nb); hnobj.ReCenterBin();
+ cout<<"hnobj: "; hnobj.Show();
+ // Histo de check du tirage
+ Histo hlflux(htirls); hlflux.Zero();
+ // Histo du flux/pixel
+ if(nbinlx==0 || nbinlx==1) nbinlx = -100;
+ if(nbinlx<0) {
+   nbinlx = -nbinlx;
+   lxmin = log10(fluxang) + lxmin;
+   lxmax = log10(fluxang) + lxmax;
+   cout<<"Recomputing: lxmin="<<lxmin<<" lxmax="<<lxmax<<" nbin="<<nbinlx<<endl;
+ }
+ Histo hlfang(lxmin,lxmax,nbinlx); hlfang.Zero();
+ cout<<"hlfang: "; hlfang.Show();
+
+ // Tirage aleatoire
+ cout<<endl;
+ if(npix>1) {
+   long l = npix/10; if(l<=0) l=1;
+   cout<<"Making "<<npix<<" trials"<<endl;
+   unsigned long nzero=0, ntir=0;
+   double sumpix=0., sumpix2=0., sumn=0., sumn2=0.;
+   for(long i=0;i<npix;i++) {
+     if(i%l==0) cout<<"...trial: "<<i<<endl;
+     unsigned long n = PoissRandLimit(mu,10.);
+     hnobj.Add((double)n);
+     sumn += n; sumn2 += n*n;
+     double fsum = 0.; unsigned long nf = 0;
+     for(long j=0;j<n;j++) {
+       double f = agn.FluxJY();
+       hlflux.Add(log10(f));
+       fsum += f;
+       nf++;
+     }
+     ntir += nf;
+     if(nf==0) {
+       nzero++;
+     } else {
+       sumpix += fsum; sumpix2 += fsum*fsum;
+       hlfang.Add(log10(fsum));
+     }
+   }
+   sumn /= (double)npix;
+   sumn2 = sumn2/(double)npix - sumn*sumn;
+   sumpix /= (double)npix;
+   sumpix2 = sumpix2/(double)npix - sumpix*sumpix;
+   cout<<"Number of try: "<<ntir<<endl;
+   cout<<"Number of pixels with no flux: "<<nzero<<" / "<<npix<<endl;
+   cout<<"Mean object number/pix: "<<sumn<<" (check "<<(double)ntir/(double)npix
+       <<") s^2="<<sumn2<<" s="<<sqrt(fabs(sumn2))<<endl;
+   cout<<"Mean flux/pix: "<<sumpix<<" s^2="<<sumpix2<<" s="<<sqrt(fabs(sumpix2))<<endl;
+   // Info de l'histo des flux
+   cout<<"hlfang: "; hlfang.Show();
+   // normalisation de hlflux
+   hlflux *= hdndls.BinWidth()/hlflux.BinWidth()*hdndls.Sum()/hlflux.Sum();
+ }
+
+ // --- Ecriture sur fichier ppf
+ string tag = "cmvtstagn.ppf";
+ POutPersist pos(tag);
+
+ TVector<r_8> vxjack(xjack), vyjack(yjack);
+ tag = "xjack"; pos.PutObject(vxjack,tag);
+ tag = "yjack"; pos.PutObject(vyjack,tag);
+
+ tag = "dndls0"; pos.PutObject(hdndls0,tag);
+ tag = "dndls"; pos.PutObject(hdndls,tag);
+
+ tag = "tirls"; pos.PutObject(htirls,tag);
+ tag = "hlflux"; pos.PutObject(hlflux,tag);
+
+ tag = "nobj"; pos.PutObject(hnobj,tag);
+ tag = "hlfang"; pos.PutObject(hlfang,tag);
+
+ return 0;
+}
+
+
+/*
+openppf cmvtstagn.ppf
+
+del ntjack
+c++exec char *vname[2] = {"x","y"}; double xnt[2]; \
+NTuple ntjack(2,vname); \
+for(int i=0;i<xjack.Size();i++) \
+  {xnt[0]=xjack(i); xnt[1]=yjack(i); ntjack.Fill(xnt);} \
+KeepObj(ntjack); cout<<"OK";
+
+# Jackson initial et son interpolation (Jy^1.5/sr)
+n/plot ntjack.y%x ! ! "nsta connectpoints"
+vecplot xjack yjack "same red fcirclemarker7"
+
+# Check de ce qu'on a extrapole
+set yxt "y*pow(10.,-1.5*x)*M_LN10"
+n/plot dndls0.log(val)%x ! ! "nsta"
+n/plot ntjack.log($yxt)%x ! ! "nsta red fcirclemarker7 same"
+
+# la pdf : histo en 1/log10(Jy)/sr
+disp dndls0 "nsta"
+n/plot ntjack.$yxt%x ! ! "nsta red fcirclemarker7 same"
+
+# la pdf : histo en 1/log10(Jy)/angsol
+disp dndls "nsta"
+
+# la pdf en flux: dn/dlog10(S)*S en microJy/log10(Jy)/sr
+n/plot dndls.val*pow(10.,x)*1e6%x ! ! "nsta connectpoints"
+
+# la fonction de tirage en 1/Jy/(angsol sr)
+disp tirls "nsta"
+
+# Check tirages
+disp dndls "nsta"
+disp hlflux "nsta same red"
+
+# Resultat: nombre d'objets (poisson) par pixel
+disp nobj
+
+# Resultat: distribution du flux par pixel
+disp hlfang
+
+ */
Index: trunk/Cosmo/SimLSS/cmvtstblack.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstblack.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cmvtstblack.cc	(revision 3196)
@@ -10,6 +10,6 @@
 #include "ntuple.h"
 
-#include "integfunc.h"
 #include "constcosmo.h"
+#include "geneutils.h"
 #include "planckspectra.h"
 
Index: trunk/Cosmo/SimLSS/cmvtstsch.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3196)
@@ -14,5 +14,4 @@
 
 #include "geneutils.h"
-#include "integfunc.h"
 #include "cosmocalc.h"
 #include "schechter.h"
@@ -113,6 +112,9 @@
  cout<<endl;
  double d = 1000.;
+ double f = Msol2FluxHI(mstar,d);
  cout<<"FluxHI: d="<<d<<" Mpc:"<<"  M= "<<mstar<<" Msol  "
-     <<"-> Flux= "<<FluxHI(mstar,d)<<" W/m^2"<<endl;
+     <<"-> Flux= "<<f<<" W/m^2"<<endl;
+ cout<<"Check... MassHI: d="<<d<<" Mpc:"<<"  f= "<<f<<" W/m^2  "
+     <<"-> Mass= "<<FluxHI2Msol(f,d)<<" Msol"<<endl;
 
  //------- Random
Index: trunk/Cosmo/SimLSS/cosmocalc.cc
===================================================================
--- trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3196)
@@ -13,5 +13,5 @@
 #include "constcosmo.h"
 #include "cosmocalc.h"
-#include "integfunc.h"
+#include "geneutils.h"
 
 
Index: trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3196)
@@ -20,5 +20,4 @@
 
 #include "constcosmo.h"
-#include "integfunc.h"
 #include "geneutils.h"
 
@@ -970,4 +969,63 @@
 }
 
+void GeneFluct3D::AddAGN(double lmsol,double lsigma)
+// Add AGN flux into simulation:
+// --- Procedure:
+// 1. lancer "cmvdefsurv" avec les parametres du survey
+//    et recuperer l'angle solide "angsol sr" du pixel elementaire
+//    au centre du cube.
+// 2. lancer "cmvtstagn" pour cet angle solide -> cmvtstagn.ppf
+// 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
+//    cad une moyenne <log10(S)> et un sigma "sig"
+//    Attention: la distribution n'etant pas gaussienne
+// 4. re-lancer "cmvdefsurv" en ajoutant l'info sur les AGN
+//    "cmvdefsurv ... -G <log10(S)> ..."
+//    et recuperer le log10(masse solaire equivalente)
+// --- Limitations actuelle du code:
+// a. l'angle solide du pixel est pris au centre du cube
+//    et ne varie pas avec la distance a l'interieur du cube
+// b. la taille en dNu des pixels (selon z) est supposee constante
+//    et egale a celle calculee pour le centre du cube
+// c. les AGN sont supposes plats en flux
+// d. le flux des AGN est mis dans une colonne Oz (indice k)
+//    et pas sur la ligne de visee
+// e. la distribution est approximee a une gaussienne
+// .. C'est une approximation pour un observateur loin du centre du cube
+//    et pour un cube peu epais / distance observateur 
+// --- Parametres de la routine:
+// lmsol : c'est le <log10(Msol)> qui est la conversion en masse solaire
+//         du flux depose dans un pixel par les AGN
+// lsigma : c'est le sigma de la distribution
+// - Comme on est en echelle log10():
+//   on tire log10(Msol) + X
+//   ou X est une realisation sur une gaussienne de variance "sig^2"
+//   La masse realisee est donc: Msol*10^X
+// - Pas de probleme de pixel negatif car on a une multiplication!
+{
+  if(lp_>0) cout<<"--- AddAGN: lmsol = "<<lmsol<<" , sigma = "<<lsigma<<endl;
+  check_array_alloc();
+
+  double m = pow(10.,lmsol);
+
+  double sum=0., sum2=0., nsum=0.;
+  for(long l=0;l<Nz_;l++) {
+    double a = lsigma*NorRand();
+    a = m*pow(10.,a);
+    // On met le meme tirage le long de Oz (indice k)
+    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
+      int_8 ip = IndexR(i,j,l);
+      data_[ip] += a;
+    }
+    sum += a; sum2 += a*a; nsum += 1.;
+  }
+
+ if(nsum>1.) {
+   sum /= nsum;
+   sum2 = sum2/nsum - sum*sum;
+   cout<<"...Mean mass="<<sum<<" Msol , s^2="<<sum2<<" s="<<sqrt(fabs(sum2))<<endl;
+ }
+ 
+}
+
 void GeneFluct3D::AddNoise2Real(double snoise)
 // add noise to every pixels (meme les <=0 !)
Index: trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3194)
+++ trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3196)
@@ -83,4 +83,5 @@
   double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
   double TurnMass2Flux(void);
+  void AddAGN(double lmsol,double lsigma);
   void AddNoise2Real(double snoise);
 
Index: trunk/Cosmo/SimLSS/geneutils.cc
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/geneutils.cc	(revision 3196)
@@ -150,4 +150,49 @@
 }
 
+//----------------------------------------------------
+double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint)
+// Interpole in x0 the table Y = f(X)
+//           X doit etre ordonne par ordre croissant (strictement)
+// typint = 0 : nearest value
+//          1 : linear interpolation
+//          2 : parabolique interpolation
+{
+ long n = X.size();
+ if(Y.size()<n || n<2)
+   throw ParmError("InterpTab_Error :  incompatible size between X and Y tables!");
+
+ if(x0<X[0] || x0>X[n-1]) return 0.;
+ if(typint>2) typint = 0;
+
+ // Recherche des indices encadrants par dichotomie
+ long k, klo=0, khi=n-1;
+ while (khi-klo > 1) {
+   k = (khi+klo) >> 1;
+   if (X[k] > x0) khi=k; else klo=k;
+ }
+
+ // Quel est le plus proche?
+ k = (x0-X[klo]<X[khi]-x0) ? klo: khi;
+
+ // On retourne le plus proche
+ if(typint==0) return Y[k];
+
+ // On retourne l'extrapolation lineaire
+ if(typint==1 || n<3)
+   return Y[klo] + (Y[khi]-Y[klo])/(X[khi]-X[klo])*(x0-X[klo]);
+
+ // On retourne l'extrapolation parabolique
+ if(k==0) k++; else if(k==n-1) k--;
+ klo = k-1; khi = k+1;
+ double x1 = X[klo]-X[k], x2 = X[khi]-X[k];
+ double y1 = Y[klo]-Y[k], y2 = Y[khi]-Y[k];
+ double den = x1*x2*(x1-x2);
+ double a = (y1*x2-y2*x1)/den;
+ double b = (y2*x1*x1-y1*x2*x2)/den;
+ x0 -= X[k];
+ return Y[k] + (a*x0+b)*x0;;
+
+}
+
 //-------------------------------------------------------------------
 int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex)
@@ -252,2 +297,163 @@
   return n;
 }
+
+//-------------------------------------------------------------------
+
+static unsigned short IntegrateFunc_GlOrder = 0;
+static vector<double> IntegrateFunc_x;
+static vector<double> IntegrateFunc_w;
+
+///////////////////////////////////////////////////////////
+//************** Integration of Functions ***************//
+///////////////////////////////////////////////////////////
+
+
+////////////////////////////////////////////////////////////////////////////////////
+double IntegrateFunc(GenericFunc& func,double xmin,double xmax
+       ,double perc,double dxinc,double dxmax,unsigned short glorder)
+// --- Integration adaptative ---
+//     Integrate[func[x], {x,xmin,xmax}]
+// ..xmin,xmax are the integration limits
+// ..dxinc is the searching increment x = xmin+i*dxinc
+//         if <0  npt = int(|dxinc|)   (si<2 alors npt=100)
+//                et dxinc = (xmax-xmin)/npt
+// ..dxmax is the maximum possible increment (if dxmax<=0 no test)
+// ---
+// Partition de [xmin,xmax] en intervalles [x(i),x(i+1)]:
+// on parcourt [xmin,xmax] par pas de "dxinc" : x(i) = xmin + i*dxinc
+// on cree un intervalle [x(i),x(i+1)]
+//     - si |f[x(i+1)] - f[x(i)]| > perc*|f[[x(i)]|
+//     - si |x(i+1)-x(i)| >= dxmax  (si dxmax>0.)
+// Dans un intervalle [x(i),x(i+1)] la fonction est integree
+//   avec une methode de gauss-legendre d'ordre "glorder"
+{
+  double signe = 1.;
+  if(xmin>xmax) {double tmp=xmax; xmax=xmin; xmin=tmp; signe=-1.;}
+
+  if(glorder==0) glorder = 4;
+  if(perc<=0.) perc=0.1;
+  if(dxinc<=0.) {int n=int(-dxinc); if(n<2) n=100; dxinc=(xmax-xmin)/n;}
+  if(glorder != IntegrateFunc_GlOrder) {
+    IntegrateFunc_GlOrder = glorder;
+    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
+  }
+
+  // Recherche des intervalles [x(i),x(i+1)]
+  int_4 ninter = 0;
+  double sum = 0., xbas=xmin, fbas=func(xbas), fabsfbas=fabs(fbas);
+  for(double x=xmin+dxinc; x<xmax+dxinc/2.; x += dxinc) {
+    double f = func(x);
+    double dx = x-xbas;
+    bool doit = false;
+    if( x>xmax ) {doit = true; x=xmax;}
+    else if( dxmax>0. && dx>dxmax ) doit = true;
+    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
+    if( ! doit ) continue;
+    double s = 0.;
+    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++)
+      s += IntegrateFunc_w[i]*func(xbas+IntegrateFunc_x[i]*dx);
+    sum += s*dx;
+    xbas = x; fbas = f; fabsfbas=fabs(fbas); ninter++;
+  }
+  //cout<<"Ninter="<<ninter<<endl;
+
+  return sum*signe;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
+         ,double perc,double dlxinc,double dlxmax,unsigned short  glorder)
+// --- Integration adaptative ---
+// Idem IntegrateFunc but integrate on logarithmic (base 10) intervals:
+//    Int[ f(x) dx ] = Int[ x*f(x) dlog10(x) ] * log(10)
+// ..lxmin,lxmax are the log10() of the integration limits
+// ..dlxinc is the searching logarithmic (base 10) increment lx = lxmin+i*dlxinc
+// ..dlxmax is the maximum possible logarithmic (base 10) increment (if dlxmax<=0 no test)
+// Remarque: to be use if "x*f(x) versus log10(x)" looks like a polynomial
+//            better than "f(x) versus x"
+// ATTENTION: la fonction func que l'on passe en argument
+//            est "func(x)" et non pas "func(log10(x))"
+{
+  double signe = 1.;
+  if(lxmin>lxmax) {double tmp=lxmax; lxmax=lxmin; lxmin=tmp; signe=-1.;}
+
+  if(glorder==0) glorder = 4;
+  if(perc<=0.) perc=0.1;
+  if(dlxinc<=0.) {int n=int(-dlxinc); if(n<2) n=100; dlxinc=(lxmax-lxmin)/n;}
+  if(glorder != IntegrateFunc_GlOrder) {
+    IntegrateFunc_GlOrder = glorder;
+    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
+  }
+
+  // Recherche des intervalles [lx(i),lx(i+1)]
+  int_4 ninter = 0;
+  double sum = 0., lxbas=lxmin, fbas=func(pow(10.,lxbas)), fabsfbas=fabs(fbas);
+  for(double lx=lxmin+dlxinc; lx<lxmax+dlxinc/2.; lx += dlxinc) {
+    double f = func(pow(10.,lx));
+    double dlx = lx-lxbas;
+    bool doit = false;
+    if( lx>lxmax ) {doit = true; lx=lxmax;}
+    else if( dlxmax>0. && dlx>dlxmax ) doit = true;
+    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
+    if( ! doit ) continue;
+    double s = 0.;
+    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++) {
+      double y = pow(10.,lxbas+IntegrateFunc_x[i]*dlx);
+      s += IntegrateFunc_w[i]*y*func(y);
+    }
+    sum += s*dlx;
+    lxbas = lx; fbas = f; fabsfbas=fabs(fbas); ninter++;
+  }
+  //cout<<"Ninter="<<ninter<<endl;
+
+  return M_LN10*sum*signe;
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+/*
+Integration des fonctions a une dimension y=f(x) par la Methode de Gauss-Legendre.
+  --> Calcul des coefficients du developpement pour [x1,x2]
+| INPUT:
+|  x1,x2 : bornes de l'intervalle (dans nbinteg.h -> x1=-0.5 x2=0.5)
+|  glorder = degre n du developpement de Gauss-Legendre
+| OUTPUT:
+|  x[] = valeur des abscisses ou l'on calcule (dim=n)
+|  w[] = valeur des coefficients associes
+| REMARQUES:
+|  - x et w seront dimensionnes a n.
+|  - l'integration est rigoureuse si sur l'intervalle d'integration
+|    la fonction f(x) peut etre approximee par un polynome
+|    de degre 2*m (monome le + haut x**(2*n-1) )
+*/
+void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1,double x2)
+{
+  if(glorder==0) return;
+  int n = (int)glorder;
+  x.resize(n,0.); w.resize(n,0.);
+
+   int m,j,i;
+   double z1,z,xm,xl,pp,p3,p2,p1;
+
+   m=(n+1)/2;
+   xm=0.5*(x2+x1);
+   xl=0.5*(x2-x1);
+   for (i=1;i<=m;i++)  {
+      z=cos(M_PI*(i-0.25)/(n+0.5));
+      do {
+         p1=1.0;
+         p2=0.0;
+         for (j=1;j<=n;j++) {
+            p3=p2;
+            p2=p1;
+            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
+         }
+         pp=n*(z*p1-p2)/(z*z-1.0);
+         z1=z;
+         z=z1-p1/pp;
+      } while (fabs(z-z1) > 3.0e-11);  // epsilon
+      x[i-1]=xm-xl*z;
+      x[n-i]=xm+xl*z;
+      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
+      w[n-i]=w[i-1];
+   }
+}
Index: trunk/Cosmo/SimLSS/geneutils.h
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.h	(revision 3194)
+++ trunk/Cosmo/SimLSS/geneutils.h	(revision 3196)
@@ -79,4 +79,7 @@
 
 //----------------------------------------------------
+double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint=0);
+
+//----------------------------------------------------
 int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex=false);
 int FuncToVec(GenericFunc& func,TVector<r_8>& h,double xmin,double xmax,bool logaxex=false);
@@ -89,3 +92,12 @@
 unsigned long PoissRandLimit(double mu,double mumax=10.);
 
+//----------------------------------------------------
+double IntegrateFunc(GenericFunc& func,double xmin,double xmax
+         ,double perc=0.1,double dxinc=-1.,double dxmax=-1.,unsigned short glorder=4);
+
+double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
+         ,double perc=0.1,double dlxinc=-1.,double dlxmax=-1.,unsigned short glorder=4);
+
+void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1=0.,double x2=1.);
+
 #endif
Index: trunk/Cosmo/SimLSS/pkspectrum.cc
===================================================================
--- trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3196)
@@ -11,5 +11,5 @@
 
 #include "constcosmo.h"
-#include "integfunc.h"
+#include "geneutils.h"
 #include "pkspectrum.h"
 
@@ -390,5 +390,5 @@
 void VarianceSpectrum::SetInteg(double dperc,double dlogkinc,double dlogkmax,unsigned short glorder)
 // ATTENTION: on n'integre pas f(k)*dk mais k*f(k)*d(log10(k))
-// see argument details in function IntegrateFuncLog (integfunc.cc)
+// see argument details in function IntegrateFuncLog (geneutils.cc)
 {
   dperc_ = dperc;  if(dperc_<=0.) dperc_ = 0.1;
Index: trunk/Cosmo/SimLSS/schechter.cc
===================================================================
--- trunk/Cosmo/SimLSS/schechter.cc	(revision 3194)
+++ trunk/Cosmo/SimLSS/schechter.cc	(revision 3196)
@@ -68,5 +68,5 @@
 ///////////////////////////////////////////////////////////
 
-double FluxHI(double m,double d)
+double Msol2FluxHI(double m,double d)
 // Input:
 //    m : masse de HI en "Msol"
@@ -88,2 +88,12 @@
   return  2.0e-28 * m / (d*d);
 }
+
+double FluxHI2Msol(double f,double d)
+// Input:
+//    f : flux total emis a 21 cm en W/m^2
+//    d : distance en "Mpc"  (si cosmo d=DLum(z))
+// Return:
+//    m : masse de HI en "Msol"
+{
+  return f *d*d / 2.0e-28;
+}
Index: trunk/Cosmo/SimLSS/schechter.h
===================================================================
--- trunk/Cosmo/SimLSS/schechter.h	(revision 3194)
+++ trunk/Cosmo/SimLSS/schechter.h	(revision 3196)
@@ -25,5 +25,6 @@
 
 //-----------------------------------------------------------------------------------
-double FluxHI(double m,double d);
+double Msol2FluxHI(double m,double d);
+double FluxHI2Msol(double f,double d);
 
 #endif
