Index: trunk/Cosmo/SimLSS/Makefile
===================================================================
--- trunk/Cosmo/SimLSS/Makefile	(revision 3115)
+++ trunk/Cosmo/SimLSS/Makefile	(revision 3115)
@@ -0,0 +1,171 @@
+include $(SOPHYABASE)/include/sophyamake.inc
+#include ../BuildMgr/sophyamake.inc
+
+OBJ = ${CMVPROG}/obj/
+EXE = ${CMVPROG}/exe/
+LIB = ${CMVPROG}/lib/
+MACH = $(shell echo `uname`)
+ifeq ($(MACH),OSF1)
+CXXREP = -ptr $(OBJ)/CmvBAO_cxxrep/
+endif
+
+MYEXTINC = ${EXTLIBDIR}/Include
+MYLIB = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lfftw3 -lm
+#MYLIB = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lfftw3_threads -lfftw3 -lm
+
+#---- Les programmes utilitaires de calcul de cartes
+PROGS = \
+       $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
+       $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
+       $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
+       $(EXE)cmvtluc $(EXE)cmvtpoisson $(EXE)cmvconchprof
+ 
+PROGSOBJ = \
+          $(OBJ)cmvtuniv.o $(OBJ)cmvtransf.o $(OBJ)cmvtgrowth.o $(OBJ)cmvtstpk.o \
+          $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
+          $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o \
+          $(OBJ)cmvtluc.o $(OBJ)cmvtpoisson.o $(OBJ)cmvconchprof.o
+ 
+LIBROBJ = \
+          $(OBJ)cosmocalc.o $(OBJ)pkspectrum.o $(OBJ)schechter.o \
+          $(OBJ)planckspectra.o $(OBJ)integfunc.o $(OBJ)geneutils.o \
+          $(OBJ)genefluct3d.o
+
+LIBR = $(LIB)libcmvsimbao.a
+
+#----
+all: lib prog
+
+lib: $(LIBR)
+
+prog: $(PROGS)
+
+clean:
+	rm -f *.o
+	rm -f $(PROGS) $(LIBR)
+	rm -f $(PROGSOBJ) $(LIBROBJ)
+	rm -rf $(OBJ)/CmvBAO_cxxrep/
+
+##############################################################################
+$(LIB)libcmvsimbao.a : $(LIBROBJ)
+	$(AR) $(ARFLAGS) $@ $? 
+
+$(OBJ)cosmocalc.o: cosmocalc.cc cosmocalc.h constcosmo.h
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cosmocalc.cc
+$(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
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ pkspectrum.cc
+$(OBJ)geneutils.o: geneutils.cc geneutils.h
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ geneutils.cc
+$(OBJ)genefluct3d.o: genefluct3d.cc genefluct3d.h
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ genefluct3d.cc
+
+##############################################################################
+cmvtuniv: $(EXE)cmvtuniv
+	echo $@ " done"
+$(EXE)cmvtuniv: $(OBJ)cmvtuniv.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtuniv.o $(MYLIB) 
+$(OBJ)cmvtuniv.o: cmvtuniv.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtuniv.cc
+
+##############################################################################
+cmvtransf: $(EXE)cmvtransf
+	echo $@ " done"
+$(EXE)cmvtransf: $(OBJ)cmvtransf.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtransf.o $(MYLIB) 
+$(OBJ)cmvtransf.o: cmvtransf.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtransf.cc
+
+##############################################################################
+cmvtgrowth: $(EXE)cmvtgrowth
+	echo $@ " done"
+$(EXE)cmvtgrowth: $(OBJ)cmvtgrowth.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtgrowth.o $(MYLIB) 
+$(OBJ)cmvtgrowth.o: cmvtgrowth.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtgrowth.cc
+
+##############################################################################
+cmvtstpk: $(EXE)cmvtstpk
+	echo $@ " done"
+$(EXE)cmvtstpk: $(OBJ)cmvtstpk.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstpk.o $(MYLIB) 
+$(OBJ)cmvtstpk.o: cmvtstpk.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtstpk.cc
+
+##############################################################################
+cmvtstsch: $(EXE)cmvtstsch
+	echo $@ " done"
+$(EXE)cmvtstsch: $(OBJ)cmvtstsch.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstsch.o $(MYLIB) 
+$(OBJ)cmvtstsch.o: cmvtstsch.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtstsch.cc
+
+##############################################################################
+cmvtstblack: $(EXE)cmvtstblack
+	echo $@ " done"
+$(EXE)cmvtstblack: $(OBJ)cmvtstblack.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstblack.o $(MYLIB) 
+$(OBJ)cmvtstblack.o: cmvtstblack.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtstblack.cc
+
+##############################################################################
+cmvtvarspec: $(EXE)cmvtvarspec
+	echo $@ " done"
+$(EXE)cmvtvarspec: $(OBJ)cmvtvarspec.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtvarspec.o $(MYLIB) 
+$(OBJ)cmvtvarspec.o: cmvtvarspec.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtvarspec.cc
+
+##############################################################################
+cmvdefsurv: $(EXE)cmvdefsurv
+	echo $@ " done"
+$(EXE)cmvdefsurv: $(OBJ)cmvdefsurv.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvdefsurv.o $(MYLIB) 
+$(OBJ)cmvdefsurv.o: cmvdefsurv.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvdefsurv.cc
+
+##############################################################################
+cmvobserv3d: $(EXE)cmvobserv3d
+	echo $@ " done"
+$(EXE)cmvobserv3d: $(OBJ)cmvobserv3d.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvobserv3d.o $(MYLIB) 
+$(OBJ)cmvobserv3d.o: cmvobserv3d.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvobserv3d.cc
+
+##############################################################################
+cmvtintfun: $(EXE)cmvtintfun
+	echo $@ " done"
+$(EXE)cmvtintfun: $(OBJ)cmvtintfun.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtintfun.o $(MYLIB) 
+$(OBJ)cmvtintfun.o: cmvtintfun.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtintfun.cc
+
+##############################################################################
+cmvconchprof: $(EXE)cmvconchprof
+	echo $@ " done"
+$(EXE)cmvconchprof: $(OBJ)cmvconchprof.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvconchprof.o $(MYLIB) 
+$(OBJ)cmvconchprof.o: cmvconchprof.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvconchprof.cc
+
+
+##############################################################################
+cmvtpoisson: $(EXE)cmvtpoisson
+	echo $@ " done"
+$(EXE)cmvtpoisson: $(OBJ)cmvtpoisson.o $(LIB)libcmvsimbao.a
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtpoisson.o $(MYLIB) 
+$(OBJ)cmvtpoisson.o: cmvtpoisson.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtpoisson.cc
+
+##############################################################################
+cmvtluc: $(EXE)cmvtluc
+	echo $@ " done"
+$(EXE)cmvtluc: $(OBJ)cmvtluc.o
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtluc.o $(SOPHYASLBLIST) -L$(SOPHYALIBP) -lLUC
+$(OBJ)cmvtluc.o: cmvtluc.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtluc.cc
Index: trunk/Cosmo/SimLSS/arrctcast.h
===================================================================
--- trunk/Cosmo/SimLSS/arrctcast.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/arrctcast.h	(revision 3115)
@@ -0,0 +1,105 @@
+#ifndef ARRCTCAST_SEEN
+#define ARRCTCAST_SEEN
+
+#include "tarray.h"
+#include "sopnamsp.h"
+
+//-----------------------------------------------------------------------
+// Classe de Bridge et fonctions utilitaires pour permettre des 
+// cast de type de donnees - sans conversion - pour des tableaux (TArray)
+//                    R. Ansari - Oct 2006
+//-----------------------------------------------------------------------
+
+//******************************************
+// Classe NDCBridge : Bridge sur NDataBlock, utile pour cast - sans conversion - de type 
+
+template <class T>
+class NDCBridge : public Bridge {
+ public:
+  NDCBridge(NDataBlock< T > & db) : dbcs_(db) { 
+  //DBG cout << "+++ NDCBridge() this= " << hex << (long int)this << dec << endl; 
+  }
+  virtual ~NDCBridge() {
+  //DBG cout << "--- ~NDCBridge() this= " << hex << (long int)this << dec << endl; 
+  }
+ protected:
+  NDataBlock< T > dbcs_;
+};
+
+//******************************************
+// Fonction de cast de type de conteny - sans conversion - de tableaux
+// pour tableau packe en memoire
+template <class T1, class T2>								    
+TArray< T2 > ArrayCast( TArray< T1 > & a, T2 ct, sa_size_t doff=0, sa_size_t dstep=1) 
+{
+  NDCBridge< T1 > * db = new NDCBridge< T1 >(a.DataBlock());
+  //      if (a.AvgStep() != 1)
+  if (a.IsPacked() == false)
+    throw ParmError("ArrayCast - Pb: Not a packed array");
+  float rs = (float)sizeof(T1)/(float)sizeof(T2);
+  //DBG cout << "--DBG-ArrayCast SizeOf T1= " << sizeof(T1)
+  //DBG     << " T2= " << sizeof(T2) << " rs=T1/T2= " << rs << endl;
+  
+  int_4 ndim = a.Rank();
+  sa_size_t sz[BASEARRAY_MAXNDIMS];
+  for(int_4 k=0; k<BASEARRAY_MAXNDIMS; k++)   sz[k] = 0;
+  for(int_4 k=0; k<ndim; k++)  sz[k] = a.Size(k);
+  // Pour gerer correctement les vecteurs colonnes
+  int kszsc = 0;
+  if ((sz[0] == 1) && (sz[1] > 1))   kszsc = 1;
+  if (dstep > 1) 
+    sz[kszsc] = sz[kszsc]*sizeof(T1)/sizeof(T2)/dstep;
+  else 
+    sz[kszsc] = sz[kszsc]*sizeof(T1)/sizeof(T2);
+  //DBG  cout << " --DBG-ArrayCast SzX=" << a.SizeX() 
+  //     << " rs=" << rs << " dstep=" << dstep << " --> SzX=" << sz[0] << endl;
+  
+  sa_size_t step = dstep; 
+  sa_size_t off = doff; 
+  TArray< T2 > ac(ndim, sz, (T2 *)(a.Data()), step, off, db);
+  ac.SetMemoryMapping(a.GetMemoryMapping());
+  ac.SetTemp(true);
+  return ac;
+}
+
+
+//******************************************
+// cast de contenu - sans conversion - tableau complex vers float
+template <class T>
+TArray<T> ArrCastC2R(TArray< complex<T> > & a)  
+{
+  T x;
+  return ArrayCast(a, x);
+}
+
+//******************************************
+// cast de contenu - sans conversion - tableau float vers complex
+template <class T>
+TArray< complex<T> > ArrCastR2C(TArray< T > & a)  
+{
+  complex<T> x;
+  //  return ArrayCast< TArray< T > , complex<T> > (a, x);
+  return ArrayCast(a, x);
+}
+
+
+//******************************************
+// Acces a la partie reelle d'un tableau complex - en shared data
+template <class T>
+TArray<T> SDRealPart(TArray< complex<T> > & a)  
+{
+  T x;
+  return ArrayCast(a, x, 0, 2);
+}
+
+//******************************************
+// Acces a la partie imaginaire d'un tableau complex - en shared data
+template <class T>
+TArray<T> SDImagPart(TArray< complex<T> > & a)  
+{
+  T x;
+  return ArrayCast(a, x, 1, 2);
+}
+
+
+#endif
Index: trunk/Cosmo/SimLSS/cmvconchprof.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvconchprof.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvconchprof.cc	(revision 3115)
@@ -0,0 +1,91 @@
+#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 "hisprof.h"
+
+void usage(void);
+void usage(void)
+{
+ cout<<"cmvconchprof -n namehprof file1.ppf file2.ppf ..."<<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ string nameh = ""; // "hpkgen" "hpkgenf" "hpkrec"
+ vector<string> ppfname;
+
+ // --- Decodage arguments
+ char c;
+ while((c = getopt(narg,arg,"hn:")) != -1) {
+  switch (c) {
+  case 'n' :
+    nameh = optarg;
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+ if(nameh.size()<=0 || optind>=narg) {usage(); return -1;}
+ for (int i=optind;i<narg;i++) ppfname.push_back(arg[i]);
+ cout<<"nameh = "<<nameh<<" , number of file to be read = "<<ppfname.size()<<endl;
+
+ // --- lecture et moyenne des HProf
+ HProf *hpconc = NULL;
+ int nhp = 0;
+ int itest = 0;
+ double sum=0., sum2=0., nsum=0;
+ for (int ifile=0;ifile<ppfname.size();ifile++) {
+   PInPersist pis(ppfname[ifile].c_str());
+   HProf hp;
+   pis.GetObject(hp,nameh);
+   int n = hp.NBins();
+   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<n<<endl;
+   if(ifile==0) { 
+     hpconc = new HProf(hp.XMin(),hp.XMax(),hp.NBins());
+     hpconc->SetErrOpt(true);
+     nhp = n;
+     itest = int(hpconc->NBins()/5.+0.5);
+   }
+   if(n!=nhp) {cout<<"BAD NUMBER OF BINS"<<endl; continue;}
+   for(int i=0;i<nhp;i++) hpconc->AddBin(i,hp(i));
+   sum+=hp(itest); sum2+=hp(itest)*hp(itest), nsum++;
+ }
+ hpconc->UpdateHisto(true);
+ if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
+ cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
+ cout<<"         "<<" mean="<<(*hpconc)(itest)<<" sigma^2="<<hpconc->Error2(itest)<<endl;
+
+ // --- ecriture ascii
+ if(1) {
+   FILE * fdata = fopen("cmvconchprof.data","w");
+   fprintf(fdata,"# k(Mpc^-1) <pkz_reconstruit> variance_pkz\n");
+   for(int i=0;i<hpconc->NBins();i++) {
+     double k = hpconc->BinCenter(i);
+     fprintf(fdata,"%e %e %e\n",k,(*hpconc)(i),hpconc->Error2(i));
+   }
+   fclose(fdata);
+ }
+
+ // --- ecriture ppf
+ string tagobs = "cmvconchprof.ppf";
+ POutPersist posobs(tagobs);
+ tagobs = "hpconc"; posobs.PutObject(*hpconc,tagobs);
+
+ delete hpconc;
+
+ return 0;
+}
+
+/*
+openppf cmvconchprof.ppf
+
+disp hpconc
+n/plot hpconc.val%log10(x) x>0 ! "connectpoints"
+*/
Index: trunk/Cosmo/SimLSS/cmvdefsurv.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvdefsurv.cc	(revision 3115)
@@ -0,0 +1,334 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "constcosmo.h"
+#include "cosmocalc.h"
+#include "geneutils.h"
+#include "integfunc.h"
+#include "schechter.h"
+#include "planckspectra.h"
+
+inline double rad2deg(double trad) {return trad/M_PI*180.;}
+inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
+inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
+inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
+inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
+inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}
+
+void usage(void);
+void usage(void) {
+  cout<<"cmvdefsurv [-r] -x adtx,atxlarg [-y adty,atylarg] -z dred,redlarg redshift"<<endl
+      <<" -x adtx,atxlarg : resolution en Theta_x (arcmin), largeur (degre)"<<endl
+      <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
+      <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
+      <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
+      <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
+      <<" -A Tbruit : temperature de bruit (K)"<<endl
+      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
+      <<"                   (indnu==0 no evolution with freq.)"<<endl
+      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
+      <<" -F  : HI flux factor to be applied for our redshift"<<endl
+      <<" -1 : on ne mesure qu\'une polarisation"<<endl
+      <<" redshift : redshift moyen du survey"<<endl
+      <<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ // --- Valeurs fixes
+ // WMAP
+ unsigned short flat = 0;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
+ // Schechter
+ double h75 = h100 / 0.75;
+ double nstar = 0.006*pow(h75,3.);  //
+ double mstar = pow(10.,9.8/(h75*h75));  // MSol
+ double alpha = -1.37;
+ cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
+
+
+ // --- Arguments
+ bool inmpc = false;
+ double adtx=1., atxlarg=90., adty=-1., atylarg=-1.;
+ double dx=1.,txlarg=1000., dy=-1.,tylarg=1000., dz=1.,tzlarg=100.;
+ int nx,ny,nz;
+ double redshift = 1., dred=0.01, redlarg=0.3;
+ double tobs = 3600., surfeff = 400000.;
+ double Tbruit=75.;
+ // a 408 MHz (Haslam) + evol index a -2.6
+ double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
+ double mhiref = -1.; // reference Mass en HI (def integ schechter)
+ double hifactor = 1.;
+ double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5
+
+ // --- Decodage arguments
+ char c;
+  while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
+  switch (c) {
+  case 'P' :
+    inmpc = true;
+    break;
+  case 'x' :
+    sscanf(optarg,"%lf,%lf",&adtx,&atxlarg);
+    break;
+  case 'y' :
+    sscanf(optarg,"%lf,%lf",&adty,&atylarg);
+    break;
+  case 'z' :
+    sscanf(optarg,"%lf,%lf",&dred,&redlarg);
+    break;
+  case 'O' :
+    sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
+    break;
+  case 'A' :
+    sscanf(optarg,"%lf",&Tbruit);
+    break;
+  case 'S' :
+    sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
+    break;
+  case 'M' :
+    sscanf(optarg,"%lf",&mhiref);
+    break;
+  case 'F' :
+    sscanf(optarg,"%lf",&hifactor);
+    break;
+  case '1' :
+    facpolar = 0.5;
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+ if(optind>=narg) {usage(); return-1;}
+ sscanf(arg[optind],"%lf",&redshift);
+ if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}
+
+ // --- Initialisation de la Cosmologie
+ cout<<"\n--- Cosmology for z = "<<redshift<<endl;
+ CosmoCalc univ(flat,true,2.*redshift);
+ double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
+ univ.SetInteg(perc,dzinc,dzmax,glorder);
+ univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print(redshift);
+
+ double dang    = univ.Dang(redshift);
+ double dtrcom  = univ.Dtrcom(redshift);
+ double dlum    = univ.Dlum(redshift);
+ double dloscom = univ.Dloscom(redshift);
+ double dlosdz  = univ.Dhubble()/univ.E(redshift);
+ cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
+     <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;
+
+ cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
+ cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
+ cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;
+
+ cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;
+
+ cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
+ cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
+     <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;
+
+ // --- Mise en forme dans les unites appropriees
+ if(adty<=0.) adty=adtx;
+ if(atylarg<=0.) atylarg=atxlarg;
+ if(inmpc) {
+   dx = adtx; txlarg = atxlarg; nx = int(txlarg/dx+0.5);
+   dy = adty; txlarg = atxlarg; ny = int(tylarg/dy+0.5);
+   dz = dred; tzlarg = redlarg; nz = int(tzlarg/dz+0.5);
+   adtx = dx/dtrcom; atxlarg = adtx*nx;
+   adty = dy/dtrcom; atylarg = adty*ny;
+   dred = dz/dlosdz; redlarg = dred*nz;
+ } else {
+   adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg); nx = int(atxlarg/adtx+0.5);
+   adty = min2rad(adty); atylarg = deg2rad(atylarg); ny = int(atylarg/adty+0.5);
+   nz = int(redlarg/dred+0.5);
+   dx = adtx*dtrcom; txlarg = dx*nx;
+   dy = adty*dtrcom; tylarg = dy*ny;
+   dz = dred*dlosdz; tzlarg = dz*nz;
+ }
+ double Npix = (double)nx*(double)ny*(double)nz;
+
+ double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
+ if(redlim[0]<=0.)
+   {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
+ double dtrlim[2] = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
+ double loslim[2] = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
+ double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};
+
+ cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;
+ cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
+     <<"dred = "<<dred<<"  redlarg = "<<redlarg<<endl
+     <<"  dz = "<<dz<<" Mpc  redlarg = "<<tzlarg<<" Mpc, nz = "<<nz<<" pix"<<endl;
+ cout<<"---- Transverse X:"<<endl
+     <<"adtx = "<<rad2min(adtx)<<"\',  atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
+     <<"  dx = "<<dx<<" Mpc,  txlarg = "<<txlarg<<" Mpc, nx = "<<nx<<" pix"<<endl;
+ cout<<"---- Transverse Y:"<<endl
+     <<"adty = "<<rad2min(adty)<<"\',  atylarg = "<<rad2deg(atylarg)<<" d"<<endl
+     <<"  dy = "<<dy<<" Mpc,  tylarg = "<<tylarg<<" Mpc, ny = "<<ny<<" pix"<<endl;
+ cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;
+
+ // --- Cosmolographie Transverse
+ cout<<"\n--- Transverse"<<endl;
+ cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
+                         <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;
+
+ cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
+     <<"       with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
+ cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
+     <<"          with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;
+
+ cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
+     <<"       with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
+ cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
+     <<"          with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;
+
+ // --- Cosmolographie Line of sight
+ cout<<"\n--- Line of Sight"<<endl;
+ cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
+     <<loslim[0]<<","<<loslim[1]<<"]"<<endl
+     <<"                                    diff = "
+     <<loslim[1]-loslim[0]<<" Mpc"<<endl;
+
+ cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
+ cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;
+
+ // --- Solid Angle & Volume
+ cout<<"\n--- Solid angle"<<endl;
+ double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
+ cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
+ double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
+ cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;
+
+ cout<<"\n--- Volume"<<endl;
+ double dvol = dx*dy*dz;
+ cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
+ double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
+ cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
+     <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;
+
+ // --- Fourier space: k = omega = 2*Pi*Nu
+ cout<<"\n--- Fourier space"<<endl;
+ cout<<"Array size: nx = "<<nx<<",  ny = "<<ny<<",  nz = "<<nz<<endl;
+ double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
+ double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
+ double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
+ cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1  (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
+     <<"            dk_y = "<<dk_y<<" Mpc^-1  (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
+ cout<<"Nyquist:    kx = "<<knyq_x<<" Mpc^-1  (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
+     <<"            ky = "<<knyq_y<<" Mpc^-1  (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
+ cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1  (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
+ cout<<"Nyquist:    kz = "<<knyq_z<<" Mpc^-1  (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;
+
+ // --- Masse de HI
+ cout<<"\n--- Mass HI"<<endl;
+ Schechter sch(nstar,mstar,alpha);
+ sch.SetOutValue(1);
+ cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
+ cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
+ int npt = 10000;
+ double lnx1=log10(1e-6), lnx2=log10(1e+14), dlnx=(lnx2-lnx1)/npt;
+ double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
+ cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;
+
+ double masshipix = masshimpc3*dvol;
+ double masshitot = masshimpc3*vol;
+ cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
+     <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
+ if(mhiref<=0.)  mhiref = masshipix;
+
+ // --- Survey values
+ double unplusz = 1.+redshift;
+ double nuhiz = Fr_HyperFin_Par / unplusz;  // GHz
+ // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
+ //     = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
+ double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
+               / (1.- (dred/.2/unplusz)*(dred/.2/unplusz));
+ cout<<"\n--- Observation:"<<endl
+     <<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
+     <<"      nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
+ cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;
+
+ // --- Power emitted by HI
+ cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
+ cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;
+
+ double fhi = hifactor*FluxHI(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;
+ double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;
+
+ // --- Signal analysis
+ cout<<"\n--- Signal analysis"<<endl;
+ PlanckSpectra planck(T_CMB_Par);
+ planck.SetApprox(1);  // Rayleigh spectra
+ planck.SetVar(0); // frequency
+ planck.SetUnitOut(0); // output en W/....
+ planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz
+
+ double hnu = h_Planck_Cst*(nuhiz*1.e+9);
+ cout<<"Facteur polar = "<<facpolar<<",   h.nu="<<hnu<<" J"<<endl;
+
+ double psig = facpolar * fhi * surfeff;
+ double esig = psig * tobs;
+ double nsig = esig / hnu;
+ double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
+ double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
+ cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
+ cout<<"     flux density = "<<ssig<<" Jy"<<endl;
+
+ double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);
+
+ double tsynch = Tsynch408;
+        if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
+ planck.SetTemperature(tsynch);
+ double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
+ double esynch = psynch * tobs;
+ double nsynch = esynch / hnu;
+ double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
+     <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
+     <<"             P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
+ cout<<"     flux density = "<<ssynch<<" Jy"<<endl;
+
+ double tcmb = T_CMB_Par;
+ planck.SetTemperature(tcmb);
+ double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
+ double ecmb = pcmb * tobs;
+ double ncmb = ecmb / hnu;
+ double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
+ cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
+ cout<<"     flux density = "<<scmb<<" Jy"<<endl;
+
+ // --- Noise analysis
+ cout<<"\n--- Noise analysis"<<endl;
+ double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
+ double ebruit = pbruit * tobs;
+ cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;
+
+ double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
+ cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;
+
+ double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
+ cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;
+
+ double SsN = ssig/slim;
+ cout<<"\nSignal to noise ratio = "<<SsN<<endl;
+
+ double smass = mhiref/ssig*slim;
+ cout<<"\nSigma noise equivalent = "<<smass<<" Msol"<<endl;
+
+ return 0;
+}
Index: trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3115)
@@ -0,0 +1,381 @@
+#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 "matharr.h"
+#include "srandgen.h"
+#include "perandom.h"
+
+#include "arrctcast.h"
+
+#include "constcosmo.h"
+#include "schechter.h"
+#include "geneutils.h"
+#include "integfunc.h"
+#include "genefluct3d.h"
+
+void usage(void);
+void usage(void)
+{
+ cout<<"cmvobserv3d [-a] [-0]"<<endl
+     <<" -a : init auto de l aleatoire"<<endl
+     <<" -0 : use methode ComputeFourier0"<<endl
+     <<" -x nx,dx : taille en x (npix,Mpc)"<<endl
+     <<" -y ny,dy : taille en y (npix,Mpc)"<<endl
+     <<" -z nz,dz : taille en z (npix,Mpc)"<<endl
+     <<" -Z zref : redshift median"<<endl
+     <<" -s scale : sigma du bruit en Msol"<<endl
+     <<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ InitTim();
+
+ //-----------------------------------------------------------------
+ // *** Survey definition
+ int nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
+ //int nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
+ //int nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
+
+ // *** Cosmography definition   (WMAP)
+ double ob0 = 0.0444356;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double zref = 0.5;
+
+ // *** Spectrum and variance definition
+ double ns = 1., as = 1.;
+ double R=8./h100, Rg=R/sqrt(5.);
+ double sigmaR = 1.;
+
+ double kmin=1e-5,kmax=1000.;
+ int npt = 10000;
+ double lkmin=log10(kmin), lkmax=log10(kmax);
+ double eps=1.e-3;
+
+ // *** Schechter mass function definition
+ double h75 = 0.71 / 0.75;
+ double nstar = 0.006*pow(h75,3.); 
+ double mstar = pow(10.,9.8/(h75*h75));  // MSol
+ double alpha = -1.37;
+
+ double schmin=1e8, schmax=1e12;
+ int schnpt = 1000;
+ double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
+
+ // *** Niveau de bruit
+ double snoise= 0.;   // en equivalent MSol
+
+ // *** type de generation
+ bool computefourier0=false;
+ unsigned short nthread=4;
+
+ // --- Decodage arguments
+
+ char c;
+ while((c = getopt(narg,arg,"ha0x:y:z:s:Z:")) != -1) {
+  switch (c) {
+  case 'a' :
+    Auto_Ini_Ranf(5);
+    break;
+  case '0' :
+    computefourier0 = true;
+    break;
+  case 'x' :
+    sscanf(optarg,"%d,%lf",&nx,&dx);
+    break;
+  case 'y' :
+    sscanf(optarg,"%d,%lf",&ny,&dy);
+    break;
+  case 'z' :
+    sscanf(optarg,"%d,%lf",&nz,&dz);
+    break;
+  case 's' :
+    sscanf(optarg,"%lf",&snoise);
+    break;
+  case 'Z' :
+    sscanf(optarg,"%lf",&zref);
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+
+ string tagobs = "cmvobserv3d.ppf";
+ POutPersist posobs(tagobs);
+
+ cout<<"zref="<<zref<<endl;
+ cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
+ cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"
+     <<", npt="<<npt<<endl;
+ cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
+ cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
+ cout<<"schmin="<<schmin<<" ("<<lschmin
+     <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
+     <<", schnpt="<<schnpt<<endl;
+ cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
+
+ InitialSpectrum pkini(ns,as);
+
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+ //tf.SetNoOscEnv(2);
+
+ GrowthFactor d1(om0,ol0);
+
+ PkSpectrum0 pk0(pkini,tf);
+
+ PkSpectrumZ pkz(pk0,d1,zref);
+ 
+ Schechter sch(nstar,mstar,alpha);
+
+ //-----------------------------------------------------------------
+ pkz.SetZ(0.);
+ cout<<endl<<"\n--- Compute variance for top-hat R="<<R
+     <<" at z="<<pkz.GetZ()<<endl;
+ VarianceSpectrum varpk_th(pkz,0);
+ double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
+ double pkmax_th = varpk_th(kfind_th);
+ cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
+ double k1=kmin, k2=kmax;
+ int rc = varpk_th.FindLimits(R,pkmax_th/1.e4,k1,k2,eps);
+ cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
+     <<k2<<" ("<<log10(k2)<<")"<<endl;
+
+ double ldlk = (log10(k2)-log10(k1))/npt;
+ varpk_th.SetInteg(0.01,ldlk,-1.,4);
+ double sr2 = varpk_th.Variance(R,k1,k2);
+ cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
+
+ double normpkz = sigmaR*sigmaR/sr2;
+ pkz.SetScale(normpkz);
+ cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
+
+ pkz.SetZ(zref);
+
+ Histo hpkz (lkmin,lkmax,npt); hpkz.ReCenterBin();
+ FuncToHisto(pkz,hpkz,true);
+ {
+ tagobs = "hpkz"; posobs.PutObject(hpkz,tagobs);
+ }
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
+ VarianceSpectrum varpk_int(pkz,2);
+
+ double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
+ double pkmax_int = varpk_int(kfind_int);
+ cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
+ double k1int=kmin, k2int=kmax;
+ int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
+ cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
+     <<k2int<<" ("<<log10(k2int)<<")"<<endl;
+
+ double ldlkint = (log10(k2int)-log10(k1int))/npt;
+ varpk_int.SetInteg(0.01,ldlkint,-1.,4);
+ double sr2int = varpk_int.Variance(R,k1int,k2int);
+ cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
+
+ PrtTim("End of definition");
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Compute galaxies density number"<<endl;
+
+ sch.SetOutValue(0);
+ cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
+ double ngal_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
+ cout<<"Galaxy density number = "<<ngal_by_mpc3<<" /Mpc^3 between limits"<<endl;
+
+ sch.SetOutValue(1);
+ cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
+ double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
+ cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
+
+ //-----------------------------------------------------------------
+ // FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f  with e+f=0 ou 1
+ cout<<endl<<"\n--- Compute Spectrum Realization"<<endl;
+
+ pkz.SetZ(zref);
+ TArray< complex<r_8> > pkgen;
+ GeneFluct3D fluct3d(pkgen,pkz);
+ fluct3d.SetNThread(nthread);
+ fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ fluct3d.Print();
+ double knyqmax = fluct3d.GetKmax();
+ double dkmin = fluct3d.GetKinc()[0];
+
+ cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
+ // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
+ // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
+ // Vc/Vs = 0.675   ->  keff = kmax * (0.675)^(1/3) = kmax * 0.877
+ knyqmax *= 0.877;
+ ldlkint = (log10(knyqmax)-log10(k1int))/npt;
+ varpk_int.SetInteg(0.01,ldlkint,-1.,4);
+ double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax);
+ cout<<"varpk_int(<"<<knyqmax<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
+
+ cout<<"\n--- Computing a realization in Fourier space"<<endl;
+ if(computefourier0) fluct3d.ComputeFourier0();
+   else              fluct3d.ComputeFourier();
+
+ cout<<"\n--- Checking realization spectra"<<endl;
+ int nhprof = int(fluct3d.GetKmax()/dkmin+0.5);
+ HProf hpkgen(0.,fluct3d.GetKmax(),nhprof);
+ hpkgen.ReCenterBin();
+ fluct3d.ComputeSpectrum(hpkgen);
+ {
+ tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
+ }
+
+ if(1) {
+   cout<<"\n--- Computing convolution by pixel shape"<<endl;
+   fluct3d.FilterByPixel();
+
+   cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
+   HProf hpkgenf(hpkgen); hpkgenf.Zero();
+   fluct3d.ComputeSpectrum(hpkgenf);
+   {
+   tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs);
+   }
+ }
+
+ if(0) {
+   cout<<"\n--- Writing cmvobserv3dk.ppf"<<endl;
+   string tag = "cmvobserv3dk.ppf";
+   POutPersist pos(tag);
+   tag = "pkgen"; pos.PutObject(pkgen,tag);
+ }
+
+ cout<<"\n--- Computing a realization in real space"<<endl;
+ fluct3d.ComputeReal();
+ r_8 undouble=0.;
+ TArray<r_8> rgen = ArrayCast(pkgen,undouble);
+ double rmin,rmax; rgen.MinMax(rmin,rmax);
+ cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+
+ cout<<"\n--- Check mean and variance in real space"<<endl;
+ size_t nlowone = fluct3d.NumberOfBad(-1.,1e+200);
+ double rm,rs2; size_t nm;
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
+
+ if(1) {
+   cout<<"\n--- Check variance sigmaR in real space"<<endl;
+   double varr;
+   size_t nvarr = fluct3d.VarianceFrReal(R,varr);
+   cout<<"R="<<R<<" : sigmaR^2="<<varr<<" -> "<<sqrt(varr)<<",   n="<<nvarr<<endl;
+ }
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Converting fluctuations into mass"<<endl;
+ fluct3d.TurnFluct2Mass();
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ cout<<"1+rgen: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ cout<<"\n--- Converting mass into galaxy number"<<endl;
+ rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
+ cout<<rm<<" galaxies put into survey"<<endl;
+ nm = fluct3d.MeanSigma2(rm,rs2,0.);
+ cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ cout<<"\n--- Set negative pixels to BAD"<<endl;
+ nm = fluct3d.SetToVal(0.,1e+200,-999.);
+ cout<<nm<<" negative in survey set to BAD"<<endl;
+ nm = fluct3d.MeanSigma2(rm,rs2,-998.);
+ cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ cout<<"\n--- Apply poisson on galaxie number"<<endl;
+ nm = fluct3d.ApplyPoisson();
+ cout<<nm<<" galaxies into survey after poisson"<<endl;
+ nm = fluct3d.MeanSigma2(rm,rs2,-998.);
+ cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ cout<<"\n--- Convert Galaxies number to HI mass"<<endl;
+ int nhmdndm = int( (lschmax-lschmin+1.)*100. + 0.5);
+ Histo hmdndm(lschmin,lschmax,nhmdndm);
+ sch.SetOutValue(1);
+ FuncToHisto(sch,hmdndm,true);
+ FunRan tirhmdndm(hmdndm,true);
+ {
+ tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
+ Histo hdum1(tirhmdndm);
+ tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
+ }
+ double mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
+ cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
+ nm = fluct3d.MeanSigma2(rm,rs2,0.);
+ cout<<"HI mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+ cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
+
+ cout<<"\n--- Set BAD pixels to Zero"<<endl;
+ nm = fluct3d.SetToVal(-998.,1e+200,0.);
+ cout<<nm<<" BAD in survey set to zero"<<endl;
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
+ fluct3d.AddNoise2Real(snoise);
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+
+ if(0) {
+   cout<<"\n--- Writing cmvobserv3dr.ppf"<<endl;
+   string tag = "cmvobserv3dr.ppf";
+   POutPersist pos(tag);
+   tag = "rgen"; pos.PutObject(rgen,tag);
+ }
+
+ //-----------------------------------------------------------------
+ // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
+ if(1) {
+   cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
+   fluct3d.ReComputeFourier();
+   HProf hpkrec(0.,fluct3d.GetKmax(),nhprof);
+   hpkrec.ReCenterBin();
+   fluct3d.ComputeSpectrum(hpkrec);
+   tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
+ }
+
+ return 0;
+}
+
+/*
+openppf cmvobserv3dk.ppf
+openppf cmvobserv3dr.ppf
+openppf cmvobserv3d.ppf
+
+objaoper pkgen sliceyz 0
+
+n/plot hpkz.val%x x>0 ! "connectpoints"
+n/plot hpkgen.val%log10(x) x>0 ! "same red connectpoints"
+n/plot hpkgenf.val%log10(x) x>0 ! "same pink connectpoints"
+n/plot hpkrec.val%log10(x) x>0 ! "same blue connectpoints"
+
+set k pow(10.,x)
+n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
+
+defscript rgensl
+ objaoper rgen sliceyz $1
+ disp sliceyz_$1 
+endscript 
+
+rgensl 0
+
+ */
Index: trunk/Cosmo/SimLSS/cmvtgrowth.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3115)
@@ -0,0 +1,63 @@
+#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 "cosmocalc.h"
+#include "pkspectrum.h"
+
+void usage(void);
+void usage(void) {cout<<"cmvtgrowth z1,z2,dz [Omatter0,Lambda0]"<<endl;}
+
+void tstprint(CosmoCalc& univ,double z1,double z2,double dz);
+void tstspeed(CosmoCalc& univ,double z1,double z2,double dz);
+void tstntuple(CosmoCalc& univ,double z1,double z2,double dz);
+
+int main(int narg,char *arg[])
+{
+
+ double z1=0., z2=10., dz=1.;
+ if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);
+ cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
+
+ double om0=0.267804, ol0=0.73;
+ if(narg>2) sscanf(arg[2],"%lf,%lf",&om0,&ol0);
+ cout<<"Om0="<<om0<<"  Ol0="<<ol0<<endl;
+
+ GrowthFactor growth(om0,ol0);
+ cout<<"D1(z=0) = "<<growth(0.)<<endl;
+
+ const int n = 2;
+ char *vname[n] = {"z","d1"};
+ NTuple nt(n,vname);
+ double xnt[n];
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   xnt[0]  = z;
+   xnt[1]  = growth(z);
+   nt.Fill(xnt);
+ }
+
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtgrowth.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+
+ return 0;
+}
+
+
+/*
+openppf cmvtgrowth.ppf
+
+set cut 1
+set cut z<5
+
+n/plot nt.d1%z $cut ! "nsta connectpoints"
+n/plot nt.1./(1.+z)%z $cut ! "nsta connectpoints red same"
+*/
Index: trunk/Cosmo/SimLSS/cmvtintfun.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtintfun.cc	(revision 3115)
@@ -0,0 +1,61 @@
+#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 "integfunc.h"
+
+// cmvtintfun [n] [xmin,xmax,npt]
+
+class MyFunc : public GenericFunc {
+public:
+  MyFunc(double n) {n_=n;}
+  virtual ~MyFunc(void) {};
+  virtual double operator() (double x)
+          {return (n_>=0.) ? pow(x,n_) : exp(n_*x);}
+  double integ(double xmin,double xmax) {
+    return (n_>=0.) ? (pow(xmax,n_+1.)-pow(xmin,n_+1.))/(n_+1.)
+                    : (exp(n_*xmax)-exp(n_*xmin))/n_;
+  }
+protected:
+  double n_;
+};
+
+
+int main(int narg,char *arg[])
+{
+
+  double n=3.;
+  if(narg>1) sscanf(arg[1],"%lf",&n);
+  double xmin=1., xmax=10.;
+  int npt = 100;
+  if(narg>2) sscanf(arg[2],"%lf,%lf,%d",&xmin,&xmax,&npt);
+
+  double lxmin=log10(xmin), lxmax=log10(xmax);
+  double perc=0.01;
+  double dxinc=(xmax-xmin)/npt, dxmax=-1.;
+  double dlxinc=(lxmax-lxmin)/npt, dlxmax=-1.;
+  unsigned short glorder=4;
+
+  if(narg>1) sscanf(arg[1],"%lf",&n);
+  if(narg>2) sscanf(arg[2],"%lf,%lf,%d",&xmin,&xmax,&npt);
+  cout<<"n="<<n<<" xmin="<<xmin<<" xmax="<<xmax<<" npt="<<npt<<endl;
+
+  MyFunc myfunc(n);
+
+  double integ = IntegrateFunc(myfunc,xmin,xmax,perc,dxinc,dxmax,glorder);
+  cout<<"IntegrateFunc() = "<<integ<<endl;
+
+  double integlog = IntegrateFuncLog(myfunc,lxmin,lxmax,perc,dlxinc,dlxmax,glorder);
+  cout<<"IntegrateFuncLog() = "<<integlog<<endl;
+
+  cout<<"Integrate() = "<<myfunc.integ(xmin,xmax)<<endl;
+
+ return 0;
+}
Index: trunk/Cosmo/SimLSS/cmvtluc.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtluc.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtluc.cc	(revision 3115)
@@ -0,0 +1,103 @@
+#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 "luc.h"
+
+void usage(void);
+void usage(void)
+{
+  cout<<"cmvtuniv z1,z2,dz"<<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ // -- WMAP
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1., ob0=0.045, oph0=5e-5;
+ // -- ouvert matter only
+ //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1., ob0=0.045, oph0=5e-5;
+ // -- plat matter only
+ //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1., ob0=0.045, oph0=5e-5;
+
+ double z1=0., z2=10., dz=1.;
+ if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);
+
+ SimpleUniverse univ;
+ univ.SetPrtDbgLevel(1);
+
+ univ.Seth(h100);
+ univ.SetOmegaMatter(om0);
+ univ.SetOmegaLambda(ol0);
+ //univ.SetDarkEnergy(ol0,w0);
+ univ.SetOmegaRadiation(or0);
+ univ.SetOmegaBaryon(ob0);
+ univ.SetOmegaPhoton(oph0);
+
+ //univ.SetFlatUniverse_OmegaMatter();
+ //univ.SetFlatUniverse_OmegaLambda();
+
+ /*
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   univ.SetEmissionRedShift(z,0.001);
+   cout<<endl;
+   univ.Print();
+ }
+ */
+
+ /**/
+ double norm = univ.HubbleLengthMpc(); // 1.;
+ double norm3 = pow(norm,3.);
+ const int n = 11;
+ char *vname[n] = {
+   "z","hz","om","or","ol","ok","ot",
+   "da","dl","dvc","vc0"
+ };
+ NTuple nt(n,vname);
+ double xnt[n];
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   univ.SetEmissionRedShift(z,0.001);
+   xnt[0]  = univ.ZE();
+   xnt[1]  = univ.HZE();
+   xnt[2]  = univ.OmegaMatterZE();
+   xnt[3]  = univ.OmegaRadiationZE();
+   xnt[4]  = 0.; //univ.OmegaLambdaZE()+univ.OmegaDarkZE();
+   xnt[5]  = univ.OmegaCurvZE();
+   xnt[6]  = 0.; //univ.OmegaTotalZE();
+   xnt[7]  = univ.AngularDiameterDistanceMpc()/norm;
+   xnt[8]  = univ.LuminosityDistanceMpc()/norm;
+   if(xnt[8]<0.) cout<<xnt[8]<<endl;
+   xnt[9]  = 0.; //univ.dVolZE()/norm3;
+   xnt[10] = 0.; //univ.Vol4PiZE()/norm3;
+   nt.Fill(xnt);
+ }
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtluc.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+ /**/
+
+ return 0;
+}
+
+
+/*
+openppf cmvtluc.ppf
+
+set cut 1
+set cut z<100.
+
+zone
+n/plot nt.hz%z $cut ! "nsta"
+
+zone
+n/plot nt.dl%z $cut ! "nsta"
+n/plot nt.da%z $cut ! "nsta same red"
+
+*/
Index: trunk/Cosmo/SimLSS/cmvtpoisson.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtpoisson.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtpoisson.cc	(revision 3115)
@@ -0,0 +1,46 @@
+#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 "srandgen.h"
+
+#include "geneutils.h"
+
+int main(int narg,char *arg[])
+{
+  unsigned long nalea=10000000;
+  // Poisson de moyenne mu  -> mean=mu,  variance=sigma^2=mu
+  //   mumax=10  -> on est a 10*sigma de zero
+  double mu = 5, mumax=10.;
+
+  cout<<"mu="<<mu<<" mumax="<<mumax<<" nalea="<<nalea<<endl;
+
+  Auto_Ini_Ranf(5);
+  InitTim();
+  double sum=0., sum2=0.;
+  double s=0., s2=0.;
+  unsigned long pr = nalea/25;
+  for(unsigned long i=0;i<nalea;i++) {
+    double m = mu*drand01(); // on fait varier la moyenne
+    sum += m; sum2 += m*m;
+    //int a = PoissRand(m);
+    unsigned long a = PoissRandLimit(m,mumax);
+    if(i%pr==0) cout<<i<<endl;
+    s += (double)a;  s2 += (double)a*(double)a;
+  }
+  s /= (double)nalea;
+  s2 = s2/(double)nalea -s*s;
+  cout<<"Poisson: mean="<<s<<"  sigma^2="<<s2<<" -> "<<sqrt(s2)<<"  n="<<nalea<<endl;
+
+  sum /= (double)nalea;
+  sum2 = sum2/(double)nalea -sum*sum;
+  cout<<"Moyenne generee: m="<<sum<<"  s^2="<<sum2<<" -> "<<sqrt(fabs(sum2))<<endl;
+  PrtTim("Fin");
+
+  return 0;
+}
Index: trunk/Cosmo/SimLSS/cmvtransf.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtransf.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtransf.cc	(revision 3115)
@@ -0,0 +1,82 @@
+#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 "constcosmo.h"
+#include "pkspectrum.h"
+
+void usage(void);
+void usage(void) {cout<<"cmvtuniv [k1,k2,npt] [h100,Om0,Ob0,tcmb]"<<endl;}
+
+int main(int narg,char *arg[])
+{
+ double k1 = 1e-6, k2 = 10.; int_4 npt = 100;
+ double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
+
+ if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&k1,&k2,&npt);
+ if(k1<=0.) {cout<<"k1 must be >0"<<endl; return -1;}
+ if(npt<=0) npt = 100;
+ cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
+ if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
+ cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
+
+ cout<<"TransfertEisenstein with baryon"<<endl;
+ TransfertEisenstein tf(h100,Om0-Ob0,Ob0,tcmb);
+ cout<<"KPeak = "<<tf.KPeak()<<"  1/Mpc"<<endl;
+
+ cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.1"<<endl;
+ TransfertEisenstein tfnosc1(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.1
+ tfnosc1.SetNoOscEnv(1);
+
+ cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.2"<<endl;
+ TransfertEisenstein tfnosc2(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.2
+ tfnosc2.SetNoOscEnv(2);
+
+ cout<<endl<<"TransfertEisenstein without baryon"<<endl;
+ TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons
+
+ const int n = 5;
+ char *vname[n] = {"k","t","tnosc1","tnosc2","tnob"};
+ NTuple nt(n,vname);
+ double xnt[n];
+
+ double lnk1 = log10(k1), lnk2=log10(k2), dlnk=(lnk2-lnk1)/npt;
+ for(double lnk=lnk1;lnk<lnk2+dlnk/2.;lnk+=dlnk) {
+   double k = pow(10.,lnk);
+   xnt[0]  = k;
+   xnt[1]  = tf(k);
+   xnt[2]  = tfnosc1(k);
+   xnt[3]  = tfnosc2(k);
+   xnt[4]  = tfnob(k);
+   nt.Fill(xnt);
+ }
+
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtransf.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+
+ return 0;
+}
+
+/*
+openppf cmvtransf.ppf
+
+set k k
+set k log10(k)
+
+n/plot nt.t%$k ! ! "nsta crossmarker3"
+n/plot nt.tnob%$k ! ! "nsta red crossmarker3 same"
+n/plot nt.tnosc1%$k ! ! "nsta blue circlemarker3 same"
+n/plot nt.tnosc2%$k ! ! "nsta green trianglemarker3 same"
+
+n/plot nt.t/tnosc1%$k ! ! "nsta crossmarker3"
+n/plot nt.t/tnosc2%$k ! ! "nsta crossmarker3 red same"
+*/
Index: trunk/Cosmo/SimLSS/cmvtstblack.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstblack.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtstblack.cc	(revision 3115)
@@ -0,0 +1,314 @@
+#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 "integfunc.h"
+#include "constcosmo.h"
+#include "planckspectra.h"
+
+void compute_xroot(void);
+
+// cmvtstblack [T,nutest(Hz)] [deriv 0/1]
+
+int main(int narg,char *arg[])
+{
+  //{compute_xroot(); return -41;}
+ double T = 2.718, nutest=-1.;
+ unsigned short deriv = 0;
+ double x0;
+
+ if(narg>1) sscanf(arg[1],"%lf,%lf",&T,&nutest);
+ if(T<=0.) T = 2.718;
+ if(narg>2) sscanf(arg[2],"%u",&deriv);
+ if(deriv!=0) deriv = 1;
+
+ cout<<"Temperature "<<T<<" K , deriv="<<deriv<<endl;
+
+ //--------------------------------- Radiance Frequence
+ // frequence planck radiance W/..
+ PlanckSpectra pradf(T);
+ pradf.SetDeriv(deriv); pradf.SetVar(0);
+ pradf.SetApprox(0); pradf.SetTypSpectra(0); pradf.SetUnitOut(0);
+ PlanckSpectra phpradf(pradf); phpradf.SetUnitOut(1);
+
+ // frequence rayleigh radiance W/..
+ PlanckSpectra rradf(T);
+ rradf.SetDeriv(deriv); rradf.SetVar(0);
+ rradf.SetApprox(1); rradf.SetTypSpectra(0); rradf.SetUnitOut(0);
+ PlanckSpectra phrradf(rradf); phrradf.SetUnitOut(1);
+
+ // frequence wien radiance W/..
+ PlanckSpectra wradf(T);
+ wradf.SetDeriv(deriv); wradf.SetVar(0);
+ wradf.SetApprox(2); wradf.SetTypSpectra(0); wradf.SetUnitOut(0);
+ PlanckSpectra phwradf(wradf); phwradf.SetUnitOut(1);
+
+ // Check
+ if(nutest<=0.) nutest = k_Boltzman_Cst*T/h_Planck_Cst;
+ cout<<"Test for t="<<T<<" K,  nu="<<nutest<<" Hz"<<endl
+     <<"    planck = "<<pradf(nutest)<<endl
+     <<"      wien = "<<wradf(nutest)<<endl
+     <<"  rayleigh = "<<rradf(nutest)<<endl;
+
+ //--------------------------------- Densite Frequence
+ // frequence planck density W/..
+ PlanckSpectra pdensf(T);
+ pdensf.SetDeriv(deriv); pdensf.SetVar(0);
+ pdensf.SetApprox(0); pdensf.SetTypSpectra(2); pdensf.SetUnitOut(0);
+ PlanckSpectra phpdensf(pdensf); phpdensf.SetUnitOut(1);
+
+ // frequence rayleigh density W/..
+ PlanckSpectra rdensf(T);
+ rdensf.SetDeriv(deriv); rdensf.SetVar(0);
+ rdensf.SetApprox(1); rdensf.SetTypSpectra(2); rdensf.SetUnitOut(0);
+ PlanckSpectra phrdensf(rdensf); phrdensf.SetUnitOut(1);
+
+ // frequence wien density W/..
+ PlanckSpectra wdensf(T);
+ wdensf.SetDeriv(deriv); wdensf.SetVar(0);
+ wdensf.SetApprox(2); wdensf.SetTypSpectra(2); wdensf.SetUnitOut(0);
+ PlanckSpectra phwdensf(wdensf); phwdensf.SetUnitOut(1);
+
+ //--------------------------------- Radiance Longeur d'onde
+ // longueur d'onde planck radiance W/..
+ PlanckSpectra pradl(T);
+ pradl.SetDeriv(deriv); pradl.SetVar(1);
+ pradl.SetApprox(0); pradl.SetTypSpectra(0); pradl.SetUnitOut(0);
+ PlanckSpectra phpradl(pradl); phpradl.SetUnitOut(1);
+
+ // longueur d'onde rayleigh radiance W/..
+ PlanckSpectra rradl(T);
+ rradl.SetDeriv(deriv); rradl.SetVar(1);
+ rradl.SetApprox(1); rradl.SetTypSpectra(0); rradl.SetUnitOut(0);
+ PlanckSpectra phrradl(rradl); phrradl.SetUnitOut(1);
+
+ // longueur d'onde wien radiance W/..
+ PlanckSpectra wradl(T);
+ wradl.SetDeriv(deriv); wradl.SetVar(1);
+ wradl.SetApprox(2); wradl.SetTypSpectra(0); wradl.SetUnitOut(0);
+ PlanckSpectra phwradl(wradl); phwradl.SetUnitOut(1);
+
+ //--------------------------------- Densite Longeur d'onde
+ // longueur d'onde planck density W/..
+ PlanckSpectra pdensl(T);
+ pdensl.SetDeriv(deriv); pdensl.SetVar(1);
+ pdensl.SetApprox(0); pdensl.SetTypSpectra(2); pdensl.SetUnitOut(0);
+ PlanckSpectra phpdensl(pdensl); phpdensl.SetUnitOut(1);
+
+ // longueur d'onde rayleigh density W/..
+ PlanckSpectra rdensl(T);
+ rdensl.SetDeriv(deriv); rdensl.SetVar(1);
+ rdensl.SetApprox(1); rdensl.SetTypSpectra(2); rdensl.SetUnitOut(0);
+ PlanckSpectra phrdensl(rdensl); phrdensl.SetUnitOut(1);
+
+ // longueur d'onde wien density W/..
+ PlanckSpectra wdensl(T);
+ wdensl.SetDeriv(deriv); wdensl.SetVar(1);
+ wdensl.SetApprox(2); wdensl.SetTypSpectra(2); wdensl.SetUnitOut(0);
+ PlanckSpectra phwdensl(wdensl); phwdensl.SetUnitOut(1);
+
+ //--------------------------------- Le maximum des spectres de planck
+ cout<<endl;
+ double eps = 0.001;
+ x0 = pradf.WienLaw();
+ printf("Planck maximum for radiance energy: f = %e Hz  -> %e W/m^2/sr/Hz\n",x0,pradf(x0));
+ printf("                             check:   = %e\n",pradf.FindMaximum(eps));
+ x0 = phpradf.WienLaw();
+ printf("Planck maximum for radiance photon: f = %e Hz  -> %e ph/s/m^2/sr/Hz\n",x0,phpradf(x0));
+ printf("                             check:   = %e\n",phpradf.FindMaximum(eps));
+
+ x0 = pdensf.WienLaw();
+ printf("Planck maximum for density energy:  f = %e Hz  -> %e J/m^3/Hz\n",x0,pdensf(x0));
+ printf("                             check:   = %e\n",pdensf.FindMaximum(eps));
+ x0 = phpdensf.WienLaw();
+ printf("Planck maximum for density photon:  f = %e Hz  -> %e ph/m^3/Hz\n",x0,phpdensf(x0));
+ printf("                             check:   = %e\n",phpdensf.FindMaximum(eps));
+
+ x0 = pradl.WienLaw();
+ printf("Planck maximum for radiance energy: l = %e m  -> %e W/m^2/sr/m\n",x0,pradl(x0));
+ printf("                             check:   = %e\n",pradl.FindMaximum(eps));
+ x0 = phpradl.WienLaw();
+ printf("Planck maximum for radiance photon: l = %e m  -> %e ph/s/m^2/sr/m\n",x0,phpradl(x0));
+ printf("                             check:   = %e\n",phpradl.FindMaximum(eps));
+
+ x0 = pdensl.WienLaw();
+ printf("Planck maximum for density energy:  l = %e m  -> %e J/m^3/m\n",x0,pdensl(x0));
+ printf("                             check:   = %e\n",pdensl.FindMaximum(eps));
+ x0 = phpdensl.WienLaw();
+ printf("Planck maximum for density photon:  l = %e m  -> %e ph/m^3/m\n",x0,phpdensl(x0));
+ printf("                             check:   = %e\n",phpdensl.FindMaximum(eps));
+
+ //--------------------------------- La loi de Wien
+ cout<<endl;
+ x0 = wradf.WienLaw();
+ printf("Wien law for radiance energy: f = %e Hz  -> %e W/m^2/sr/Hz\n",x0,wradf(x0));
+ printf("                       check:   = %e\n",wradf.FindMaximum(eps));
+ x0 = phwradf.WienLaw();
+ printf("Wien law for radiance photon: f = %e Hz  -> %e ph/s/m^2/sr/Hz\n",x0,phwradf(x0));
+ printf("                       check:   = %e\n",phwradf.FindMaximum(eps));
+
+ x0 = wdensf.WienLaw();
+ printf("Wien law for density energy:  f = %e Hz  -> %e J/m^3/Hz\n",x0,wdensf(x0));
+ printf("                      check:    = %e\n",wdensf.FindMaximum(eps));
+ x0 = phwdensf.WienLaw();
+ printf("Wien law for density photon:  f = %e Hz  -> %e  ph/m^3/Hz\n",x0,phwdensf(x0));
+ printf("                      check:    = %e\n",phwdensf.FindMaximum(eps));
+
+ x0 = wradl.WienLaw();
+ printf("Wien law for radiance energy: l = %e m  -> %e W/m^2/sr/m\n",x0,wradl(x0));
+ printf("                       check:   = %e\n",wradl.FindMaximum(eps));
+ x0 = phwradl.WienLaw();
+ printf("Wien law for radiance photon: l = %e m  -> %e ph/s/m^2/sr/m\n",x0,phwradl(x0));
+ printf("                       check:   = %e\n",phwradl.FindMaximum(eps));
+
+ x0 = wdensl.WienLaw();
+ printf("Wien law for density energy:  l = %e m  -> %e J/m^3/m\n",x0,wdensl(x0));
+ printf("                      check:    = %e\n",wdensl.FindMaximum(eps));
+ x0 = phwdensl.WienLaw();
+ printf("Wien law for density photon:  l = %e m  -> %e ph/m^3/m\n",x0,phwdensl(x0));
+ printf("                      check:    = %e\n",phwdensl.FindMaximum(eps));
+
+ //--------------------------------- Densite totale d'energie
+ cout<<endl;
+ cout<<"Densite totale d'energie: "<<pdensf.PlanckEnergie()<<" J/m^3"<<endl;
+ if(! deriv) {
+   double x1,x2,perc,dxinc;
+   x0 = pdensf.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
+   cout<<"   check by integration: "<<IntegrateFunc(pdensf,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
+   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
+   cout<<"   check by integration: "<<IntegrateFuncLog(pdensf,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
+
+   x0 = pdensl.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
+   cout<<"   check by integration: "<<IntegrateFunc(pdensl,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
+   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
+   cout<<"   check by integration: "<<IntegrateFuncLog(pdensl,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
+ }
+
+ //--------------------------------- Densite totale de photons
+ cout<<endl;
+ cout<<"Densite totale de photons: "<<pdensf.PlanckPhoton()<<" ph/m^3"<<endl;
+ if(! deriv) {
+   double x1,x2,perc,dxinc;
+   x0 = phpdensf.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
+   cout<<"   check by integration: "<<IntegrateFunc(phpdensf,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
+   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
+   cout<<"   check by integration: "<<IntegrateFuncLog(phpdensf,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
+
+   x0 = phpdensl.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
+   cout<<"   check by integration: "<<IntegrateFunc(phpdensl,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
+   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
+   cout<<"   check by integration: "<<IntegrateFuncLog(phpdensl,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
+ }
+
+ //--------------------------------- NTuple
+ cout<<endl;
+ int npt = 1000;
+ double frac = 10000.;
+ x0 = wradf.WienLaw();
+ double xmin = x0/frac, xmax = x0*frac;
+ cout<<"xmin="<<xmin<<" Hz,  xmax="<<xmax<<" Hz"<<endl;
+ double lnx1 = log10(xmin), lnx2 = log10(xmax), dlnx = (lnx2-lnx1)/npt;
+ cout<<"lnx1="<<lnx1<<",  lnx2="<<lnx2<<", dlnx="<<dlnx<<endl;
+
+ const int n = 14;
+ char *vname[n] = {"f","l","prf","prl","pdf","pdl","rrf","rrl","rdf","rdl","wrf","wrl","wdf","wdl"};
+ NTuple nt(n,vname);
+ double xnt[n];
+ for(double lx=lnx1;lx<lnx2+dlnx/2.;lx+=dlnx) {
+   double f = pow(10.,lx);
+   double l = wradf.F2L(f);
+   xnt[0] = f;
+   xnt[1] = l;
+   xnt[2] = pradf(f);
+   xnt[3] = pradl(l);
+   xnt[4] = pdensf(f);
+   xnt[5] = pdensl(l);
+   xnt[6] = rradf(f);
+   xnt[7] = rradl(l);
+   xnt[8] = rdensf(f);
+   xnt[9] = rdensl(l);
+   xnt[10] = wradf(f);
+   xnt[11] = wradl(l);
+   xnt[12] = wdensf(f);
+   xnt[13] = wdensl(l);
+   nt.Fill(xnt);
+ }
+
+ cout<<"\n>>>> Ecriture"<<endl;
+ string tag = "cmvtstblack.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+
+ return 0;
+}
+
+//---------------------------------------
+void compute_xroot(void)
+// Pour calculer les x qui donnent le maximum de f(x) pour:
+// f(x) = x^n * exp(-x)
+// ou
+// f(x) = x^n / (exp(x)-1)
+// ou
+// f(x) = x^n * exp(x) / (exp(x)-1)^2
+{
+  double n = 3.;
+  cout<<"n="<<n<<endl;
+  int npt=10;
+  double xmin=0.1,xmax=10., dx=(xmax-xmin)/npt;
+  for(int itry=0;itry<100;itry++) {
+    double x0=xmin, v0=-1e10;
+    for(double x=xmin;x<xmax+dx/2.; x+=dx) {
+      double f = pow(x,n) * exp(-x);  // Wien ou dWein/dT
+      //double f = pow(x,n) / (exp(x)-1);  // Planck
+      //double f = pow(x,n)  * exp(x) / ((exp(x)-1)*(exp(x)-1));  // dPlanck/dT
+      if(f>v0) {v0=f; x0=x;}
+    }
+    printf("x0 = %.18f , v0=%.18f , dx=%e\n",x0,v0,dx);
+    xmin = x0-dx; xmax = x0+dx; dx=(xmax-xmin)/npt;
+    if(dx<1e-15) break;
+  }
+}
+
+
+/*
+openppf cmvtstblack.ppf
+
+# Energie
+set fac 1.
+# Photon
+set fac (6.6260693e-34*f)
+
+n/plot nt.log10(f)%log10(l)
+
+# radiance frequence
+n/plot nt.prf/${fac}%log10(f) ! ! "nsta connectpoints"
+n/plot nt.rrf/${fac}%log10(f) ! ! "nsta connectpoints red same"
+n/plot nt.wrf/${fac}%log10(f) ! ! "nsta connectpoints blue same"
+
+# radiance longeur d'onde
+n/plot nt.prl/${fac}%log10(l) ! ! "nsta connectpoints"
+n/plot nt.rrl/${fac}%log10(l) ! ! "nsta connectpoints red same"
+n/plot nt.wrl/${fac}%log10(l) ! ! "nsta connectpoints blue same"
+
+# densite frequence
+n/plot nt.pdf/${fac}%log10(f) ! ! "nsta connectpoints"
+n/plot nt.rdf/${fac}%log10(f) ! ! "nsta connectpoints red same"
+n/plot nt.wdf/${fac}%log10(f) ! ! "nsta connectpoints blue same"
+
+# densite longeur d'onde
+n/plot nt.pdl/${fac}%log10(l) ! ! "nsta connectpoints"
+n/plot nt.rdl/${fac}%log10(l) ! ! "nsta connectpoints red same"
+n/plot nt.wdl/${fac}%log10(l) ! ! "nsta connectpoints blue same"
+
+# Check:  p(l)dl = p(f)df avec dl = c/f^2 df = l^2/c df
+n/plot nt.prf%log10(f) ! ! "nsta connectpoints"
+n/plot nt.prl*2.998e8/(f*f)%log10(f) ! ! "nsta connectpoints red same"
+
+ */
Index: trunk/Cosmo/SimLSS/cmvtstpk.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3115)
@@ -0,0 +1,231 @@
+#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 "histos.h"
+#include "tvector.h"
+#include "ntuple.h"
+#include "perandom.h"
+
+#include "constcosmo.h"
+#include "pkspectrum.h"
+#include "geneutils.h"
+
+void usage(void);
+void usage(void) {
+ cout<<"cmvtstpk -k npt,lkmin,lkmax -s scale zval"<<endl
+     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
+     <<" -s scale : on multiplie pkz par scale"<<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ double ob0 = 0.0444356;
+ // -- WMAP
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ // -- ouvert matter only
+ //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ // -- plat matter only
+ //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;
+
+ double ns = 1., as = 1.;
+
+ int npt = 1000;
+ double lkmin = -4., lkmax=2.;
+ double scale = 1.;
+
+ char c;
+  while((c = getopt(narg,arg,"hk:s:")) != -1) {
+  switch (c) {
+  case 'k' :
+    sscanf(optarg,"%d,%lf,%lf",&npt,&lkmin,&lkmax);
+    if(npt<=0) npt=1000;
+    if(lkmax<lkmin) {lkmin=-4.; lkmax=2.;}
+    break;
+  case 's' :
+    sscanf(optarg,"%lf",&scale);
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+ double dlk=(lkmax-lkmin)/npt;
+
+ // Fill z value into list
+ double zval = 0.;
+ if(optind<narg) zval = atof(arg[optind]);
+
+ cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
+ cout<<"scale="<<scale<<endl;
+ cout<<"zval="<<zval<<endl;
+
+ //--------------------------
+ InitialSpectrum pkini(ns,as);
+
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+ cout<<"kpeak="<<tf.KPeak()<<endl;
+ TransfertEisenstein tfnosc2(tf); tfnosc2.SetNoOscEnv(2);
+ TransfertEisenstein tfnosc1(tf); tfnosc1.SetNoOscEnv(1);
+ TransfertEisenstein tfnob(h100,om0,0.,T_CMB_Par,true);
+
+ GrowthFactor d1(om0,ol0);
+ cout<<"GrowthFactor: "<<d1(zval)<<endl;
+
+ PkSpectrum0 pk0(pkini,tf);
+ PkSpectrum0 pk0nosc2(pkini,tfnosc2);
+ PkSpectrum0 pk0nosc1(pkini,tfnosc1);
+ PkSpectrum0 pk0nob(pkini,tfnob);
+
+ PkSpectrumZ pkz(pk0,d1,zval);
+ PkSpectrumZ pkznosc2(pk0nosc2,d1,zval);
+ PkSpectrumZ pkznosc1(pk0nosc1,d1,zval);
+ PkSpectrumZ pkznob(pk0nob,d1,zval);
+
+ //--------------------------
+ Histo hd1(0.,1000.,10000); hd1.ReCenterBin();
+ FuncToHisto(d1,hd1,false);
+
+ Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
+ FuncToHisto(pkz,hpkz,true);
+ TVector<r_8> vpkz(npt);
+ FuncToVec(pkz,vpkz,lkmin,lkmax,true);
+
+ FunRan talea(hpkz,true);
+ Histo halea(hpkz); halea.Zero();
+ int nalea = 1000000;
+ for(int i=0;i<nalea;i++) halea.Add(talea.Random());
+ halea *= hpkz.Sum()/halea.Sum();
+
+ //--------------------------
+ const int n = 14;
+ char *vname[n] = {
+   "k","pkini",
+   "tf","pk0","pk",
+   "tfnosc2","pk0nosc2","pknosc2",
+   "tfnosc1","pk0nosc1","pknosc1",
+   "tfnob","pk0nob","pknob"
+ };
+ NTuple nt(n,vname);
+ double xnt[n];
+
+ for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
+   double k = pow(10.,lk);
+   xnt[0] = k;
+   xnt[1] = pkini(k);
+   xnt[2] = tf(k);
+   xnt[3] = pk0(k);
+   xnt[4] = pkz(k,zval);
+   xnt[5] = tfnosc2(k);
+   xnt[6] = pk0nosc2(k);
+   xnt[7] = pkznosc2(k,zval);
+   xnt[8] = tfnosc1(k);
+   xnt[9] = pk0nosc1(k);
+   xnt[10] = pkznosc1(k,zval);
+   xnt[11] = tfnob(k);
+   xnt[12] = pk0nob(k);
+   xnt[13] = pkznob(k,zval);
+   nt.Fill(xnt);
+ }
+
+ //--------------------------
+ cout<<">>>> ASCII"<<endl;
+ if(1) {
+   FILE * fdata = fopen("cmvtstpk.data","w");
+   fprintf(fdata,"# z= %g : k(Mpc^-1) pkini, tf pk0 pkz, tfnosc2 pk0nosc2 pkznosc2, ",zval);
+   fprintf(fdata,"tfnosc1 pk0nosc1 pkznosc1, tfnob pk0nob pkznob\n");
+   for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
+     double k = pow(10.,lk);
+     fprintf(fdata,"%e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
+	     k,pkini(k),
+	     tf(k),pk0(k),scale*pkz(k,zval),
+	     tfnosc2(k),pk0nosc2(k),scale*pkznosc2(k,zval),
+	     tfnosc1(k),pk0nosc1(k),scale*pkznosc1(k,zval),
+	     tfnob(k),pk0nob(k),scale*pkznob(k,zval));
+   }
+   fclose(fdata);
+ }
+
+ //--------------------------
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtstpk.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+ tag = "hd1"; pos.PutObject(hd1,tag);
+ tag = "hpkz"; pos.PutObject(hpkz,tag);
+ tag = "vpkz"; pos.PutObject(vpkz,tag);
+ tag = "halea"; pos.PutObject(halea,tag);
+ return 0;
+}
+
+/*
+openppf cmvtstpk.ppf
+
+#### growth-factor
+zone
+n/plot hd1.log10(val)%log10(x) val>0.&&x>0. ! "nsta connectpoints"
+
+#### Spectre initial
+zone
+n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta crossmarker3"
+
+#### fct de transfert
+zone
+n/plot nt.tf%log10(k) ! ! "nsta connectpoints"
+n/plot nt.tfnosc2%log10(k) ! ! "nsta connectpoints same red"
+n/plot nt.tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
+n/plot nt.tfnob%log10(k) ! ! "nsta connectpoints same green"
+
+n/plot nt.tf/tfnosc2%log10(k) ! ! "nsta connectpoints red"
+n/plot nt.tf/tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
+n/plot nt.tf/tfnob%log10(k) ! ! "nsta connectpoints same red"
+addline -10 1 10 1
+
+#### Spectre a z=0
+zone
+n/plot nt.pk0%log10(k) ! ! "nsta connectpoints"
+n/plot nt.pk0nosc2%log10(k) ! ! "nsta connectpoints same red"
+n/plot nt.pk0nosc1%log10(k) ! ! "nsta connectpoints same blue"
+n/plot nt.pk0nob%log10(k) ! ! "nsta connectpoints same green"
+
+# Check
+zone 2 2
+n/plot nt.pk0/pkini-tf*tf%log10(k) pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%log10(k) pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%log10(k) pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nob/pkini-tfnob*tfnob%log10(k) pkini>0 ! "nsta crossmarker3"
+
+#### Spectre a z
+zone
+n/plot nt.pk%log10(k) ! ! "nsta connectpoints"
+n/plot nt.pknosc2%log10(k) ! ! "nsta connectpoints same red"
+n/plot nt.pknosc1%log10(k) ! ! "nsta connectpoints same blue"
+n/plot nt.pknob%log10(k) ! ! "nsta connectpoints same green"
+
+n/plot nt.pk/pknosc2%log10(k) ! ! "nsta connectpoints red"
+n/plot nt.pk/pknosc1%log10(k) ! ! "nsta connectpoints same blue"
+n/plot nt.pk/pknob%log10(k) ! ! "nsta connectpoints same green"
+addline -10 1 10 1
+
+#### Le spectre version Delta^2
+set D2 k*k*k*pk/(2*M_PI*M_PI)
+n/plot nt.$D2%log10(k)  !  ! "nsta crossmarker3 connectpoints"
+
+#### Test des transferts dans Histo et TVector
+zone 1 2
+n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"
+disp hpkz "same red"
+
+disp vpkz
+c++exec cout<<hpkz(0)<<" =?= "<<vpkz(0)<<endl;
+
+zone
+disp hpkz
+disp halea "same red"
+
+*/
Index: trunk/Cosmo/SimLSS/cmvtstsch.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3115)
@@ -0,0 +1,211 @@
+#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 "histos.h"
+#include "ntuple.h"
+#include "srandgen.h"
+#include "perandom.h"
+
+#include "geneutils.h"
+#include "integfunc.h"
+#include "cosmocalc.h"
+#include "schechter.h"
+
+void usage(void);
+void usage(void)
+{
+ cout<<"cmvtstsch -m xmin,xmax -a -i perc,dlxinc,dlxmax,glorder -n nbin -N nalea"<<endl
+     <<"  xmin,xmax : limites en masse"<<endl
+     <<"  perc : see IntegrateFunc"<<endl
+     <<"  dlnx : see IntegrateFunc"<<endl
+     <<"         si <0 nombre de point dlx = (log10(xmax)-log10(xmin))/(-dlx)"<<endl
+     <<"  dlxmax : see IntegrateFunc, si <0 -dlxmax*dlx, si ==0 no check"<<endl
+     <<"  glorder : odre gauss-legendre"<<endl
+     <<"  nbin : nombre de bins"<<endl
+     <<"  nalea : nombre de tirages aleatoires"<<endl
+     <<"  -a : init auto de l\'aleatoire"<<endl;
+}
+
+
+int main(int narg,char *arg[])
+{
+ double h75 = 0.71 / 0.75;
+ double nstar = 0.006*pow(h75,3.);  //  
+ double mstar = pow(10.,9.8/(h75*h75));  // MSol
+ double alpha = -1.37;
+ cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
+
+ double xmin=1e6, xmax=1e12;
+ double perc=0.01,dlxinc=-1000.,dlxmax=0.; unsigned short glorder=4;
+ int npt = 1000;
+ int nalea = 100000;
+
+ char c;
+ while((c = getopt(narg,arg,"ham:i:n:N:")) != -1) {
+   switch (c) {
+   case 'm' :
+     sscanf(optarg,"%lf,%lf",&xmin,&xmax);
+     break;
+   case 'i' :
+     int ig;
+     sscanf(optarg,"%lf,%lf,%lf,%u",&perc,&dlxinc,&dlxmax,&ig);
+     if(perc<=0.) perc = 0.01;
+     glorder = (ig<=0) ? 4: ig;
+     break;
+   case 'n' :
+     sscanf(optarg,"%d",&npt);
+     if(npt<=0) npt = 1000;
+     break;
+   case 'N' :
+     sscanf(optarg,"%d",&nalea);
+     break;
+   case 'a' :
+     Auto_Ini_Ranf(5);
+     break;
+   case 'h' :
+   default :
+     usage(); return -1;
+   }
+ }
+
+ double lnx1=log10(xmin), lnx2=log10(xmax), dlnx = (lnx2-lnx1)/npt;
+ cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<"), dlnx="<<dlnx<<endl;
+
+ if(dlxinc<0.) dlxinc=(lnx2-lnx1)/(-dlxinc);
+ if(dlxmax<0.) dlxmax *= -dlxinc; else if(dlxmax==0.) dlxmax=-1.;
+ cout<<"perc="<<perc<<" dlxinc="<<dlxinc<<" dlxmax="<<dlxmax<<" glorder="<<glorder<<endl;
+
+ cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;
+
+ Schechter sch(nstar,mstar,alpha);
+ cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
+
+ InitTim();
+
+ //-------dn/dm
+ cout<<endl;
+ sch.SetOutValue(0);
+ for(int i=0;i<=5;i++) {
+   double num = IntegrateFuncLog(sch,lnx1-i,lnx2+i,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated number = "<<num<<endl;
+ }
+ Histo hdndm(lnx1,lnx2,npt); hdndm.ReCenterBin();
+ FuncToHisto(sch,hdndm,true);
+
+ //-------m*dn/dm
+ cout<<endl;
+ sch.SetOutValue(1);
+ cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
+ for(int i=0;i<=5;i++) {
+   double sum = IntegrateFuncLog(sch,lnx1-i,lnx2+i,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated mass = "<<sum<<" Msol"<<endl;
+ }
+ Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
+ FuncToHisto(sch,hmdndm,true);
+
+ //------- Le flux HI pour M=mstar a 1 Gpc
+ cout<<endl;
+ double d = 1000.;
+ cout<<"FluxHI: d="<<d<<" Mpc:"<<"  M= "<<mstar<<" Msol  "
+     <<"-> Flux= "<<FluxHI(mstar,d)<<" W/m^2"<<endl;
+
+ //------- Random
+ Histo hdna(hdndm); hdna.Zero();
+ Histo hmdna(hmdndm); hmdna.Zero();
+ FunRan tirhdndm = FunRan(hdndm,true);
+ FunRan tirhmdndm = FunRan(hmdndm,true);
+ if(nalea>0) {
+   double summ=0.;
+   PrtTim("--- avant tirage aleatoire");
+   for(int i=0;i<nalea;i++) {
+     double a = tirhdndm.RandomInterp();
+     hdna.Add(a);
+     a = tirhmdndm.RandomInterp();
+     summ += pow(10.,a);
+     hmdna.Add(a);
+   }
+   cout<<"Tirgae aleatoire: <m>="<<summ/(double)nalea<<endl;
+   hdna *= hdndm.Sum()/hdna.Sum();
+   hmdna *= hmdndm.Sum()/hmdna.Sum();
+   PrtTim("--- apres tirage aleatoire");
+ }
+
+ //-------
+ const int n = 3;
+ char *vname[n] = {"m","f","mf"};
+ NTuple nt(n,vname);
+ double xnt[n];
+ for(double lx=lnx1;lx<lnx2+dlnx/2.;lx+=dlnx) {
+   double x = pow(10.,lx);
+   xnt[0]  = x;
+   sch.SetOutValue(0);
+   xnt[1]  = sch(x);
+   sch.SetOutValue(1);
+   xnt[2]  = sch(x);
+   nt.Fill(xnt);
+ }
+
+ cout<<"Ecriture"<<endl;
+ string tag = "cmvtstsch.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+ tag = "hdndm"; pos.PutObject(hdndm,tag);
+ tag = "hmdndm"; pos.PutObject(hmdndm,tag);
+ tag = "hdna"; pos.PutObject(hdna,tag);
+ tag = "hmdna"; pos.PutObject(hmdna,tag);
+ Histo hdum1(tirhdndm);
+ tag = "tirhdndm"; pos.PutObject(hdum1,tag);
+ Histo hdum2(tirhmdndm);
+ tag = "tirhmdndm"; pos.PutObject(hdum2,tag);
+
+ return 0;
+}
+
+/*
+openppf cmvtstsch.ppf
+
+#--- la fonction de masse numerique
+zone 2 2
+n/plot nt.f%m ! ! "nsta crossmarker3"
+n/plot nt.f%log10(m) ! ! "nsta crossmarker3"
+n/plot nt.log10(f)%log10(m) f>0. ! "nsta crossmarker3"
+
+# ce qu'on integre (en log10)
+n/plot nt.m*f%log10(m) ! ! "nsta crossmarker3"
+
+#--- la fonction de masse
+zone 2 2
+n/plot nt.mf%m ! ! "nsta crossmarker3"
+n/plot nt.mf%log10(m) f>0. ! "nsta crossmarker3"
+n/plot nt.log10(mf)%log10(m) mf>0. ! "nsta crossmarker3"
+
+# ce qu'on integre (en log10)
+n/plot nt.m*mf%log10(m) ! ! "nsta crossmarker3"
+
+#--- Check histos
+zone 2 1
+disp hdndm "red"
+n/plot nt.f%log10(m) ! ! "nsta crossmarker3 same"
+
+disp hmdndm "red"
+n/plot nt.m*f%log10(m) ! ! "nsta crossmarker3 same"
+
+#--- Check tirage aleatoire
+zone 2 1
+disp hdndm
+disp hdna "same red"
+
+disp hmdndm
+disp hmdna "same red"
+
+zone
+disp tirhdndm
+disp tirhmdndm "same red"
+
+ */
Index: trunk/Cosmo/SimLSS/cmvtuniv.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3115)
@@ -0,0 +1,248 @@
+#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 "cosmocalc.h"
+#include "geneutils.h"
+
+void usage(void);
+void usage(void) {cout<<"cmvtuniv z1,z2,dz [perc,dzinc,dzmax,glorder]"<<endl;}
+
+void tstprint(CosmoCalc& univ,double z1,double z2,double dz);
+void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz);
+void tstspeed(CosmoCalc& univ,double z1,double z2,double dz);
+void tstntuple(CosmoCalc& univ,double z1,double z2,double dz);
+void tstinterp(CosmoCalc& univ,double z1,double z2,double dz);
+
+const double deg2rad = M_PI/180.;
+
+int main(int narg,char *arg[])
+{
+ unsigned short flat = 0;
+ // -- WMAP
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ // -- ouvert matter only
+ //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ // -- plat matter only
+ //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.; flat = 1;
+ // -- plat lambda only
+ //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.; flat = 2;
+
+ double z1=0., z2=-1., dz=-1.;
+ if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);
+ if(z2<=z1) z2 = z1+1.;
+ if(dz<0.) dz = 10.*(z2-z1); 
+ cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
+
+ double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
+ if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%u",&perc,&dzinc,&dzmax,&glorder);
+ if(dzinc<=0.) dzinc = dz/2.;
+ cout<<" perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
+
+ CosmoCalc univ(flat,true,z2);
+ univ.SetInteg(perc,dzinc,dzmax,glorder);
+ univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print();
+
+ CosmoCalc univ1(univ);
+ univ1.SetInteg(perc,dzinc,dzmax,glorder);
+ univ1.SetDynParam(h100,om0,or0,ol0,w0);
+ univ1.Print();
+
+ //tstprint(univ,z1,z2,dz);
+ //tstprint(univ,univ1,z1,z2,dz);
+ //tstspeed(univ,z1,z2,dz);
+ //tstntuple(univ,z1,z2,dz);
+ tstinterp(univ,z1,z2,dz);
+
+ return 0;
+}
+
+/* ----------------------------------------------------- */
+void tstprint(CosmoCalc& univ,double z1,double z2,double dz)
+{
+ cout<<"\nTSTPRINT()"<<endl;
+ double dzz = dz; if(dzz>z2-z1) dzz = 0.01;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   cout<<"--"<<endl;
+   univ.Print(z);
+   cout<<"Volume comoving in [z,z+"<<dzz<<"] for 4Pi sr: "
+       <<univ.Vol4Pi(z,z+dzz)<<" Mpc^3"<<endl;
+   double dang = univ.Dang(z);
+   double a = deg2rad/3600.;
+   cout<<"1\" -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
+   a = deg2rad/60.;
+   cout<<"1\' -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
+   a = deg2rad;
+   cout<<"1d -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
+   double dloscom = univ.Dhubble() / univ.E(z);
+   cout<<"dz=1 -> "<<dloscom<<" Mpc com"<<endl;
+
+ }
+}
+
+/* ----------------------------------------------------- */
+void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz)
+{
+ cout<<"\nTSTPRINT()"<<endl;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   cout<<endl<<"--spline:"<<endl;
+   univ1.Print(z);
+   cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
+       <<univ1.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
+   cout<<"--integ:"<<endl;
+   univ2.Print(z);
+   cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
+       <<univ2.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
+ }
+}
+
+/* ----------------------------------------------------- */
+void tstspeed(CosmoCalc& univ,double z1,double z2,double dz)
+{
+ cout<<"\nTSTSPEED()"<<endl;
+ double sum = 0.;
+ int_4 n=0;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   sum += univ.Dang(z);
+   n++;
+ }
+ if(n>0) cout<<"... "<<sum/n<<endl;
+}
+
+/* ----------------------------------------------------- */
+void tstntuple(CosmoCalc& univ,double z1,double z2,double dz)
+{
+ cout<<"\nTSTNTUPLE()"<<endl;
+ double norm =  univ.Dhubble();  // 1.
+ double norm3 = pow(norm,3.);
+ const int n = 15;
+ char *vname[n] = {
+   "z","hz","om","or","ol","ok","ot","ob",
+   "dtc","dlc","da","dl",
+   "dvc","vc0","vcdz"
+ };
+ NTuple nt(n,vname);
+ double xnt[n];
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   xnt[0]  = z;
+   xnt[1]  = univ.H(z);
+   xnt[2]  = univ.Omatter(z);
+   xnt[3]  = univ.Orelat(z);
+   xnt[4]  = univ.Olambda(z);
+   xnt[5]  = univ.Ocurv(z);
+   xnt[6]  = univ.Otot(z);
+   xnt[7]  = univ.Obaryon(z);
+   xnt[8]  = univ.Dtrcom(z)/norm;
+   xnt[9]  = univ.Dloscom(z)/norm;
+   xnt[10] = univ.Dang(z)/norm;
+   xnt[11] = univ.Dlum(z)/norm;
+   xnt[12] = univ.dVol(z)/norm3;
+   xnt[13] = univ.Vol4Pi(z)/norm3;
+   xnt[14] = univ.Vol4Pi(z,z+dz)/norm3;
+   nt.Fill(xnt);
+ }
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtuniv.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+}
+
+
+/*
+openppf cmvtuniv.ppf
+
+set cut 1
+set cut z<100.
+set cut z<10.
+
+zone
+n/plot nt.hz%z $cut ! "nsta connectpoints"
+
+zone
+n/plot nt.dl%z $cut ! "nsta connectpoints"
+n/plot nt.da%z $cut ! "nsta connectpoints same red"
+n/plot nt.dtc%z $cut ! "nsta connectpoints same blue"
+n/plot nt.dlc%z $cut ! "nsta connectpoints same green"
+
+zone 2 2
+n/plot nt.dvc%z $cut ! "nsta connectpoints"
+n/plot nt.vc0%z $cut ! "nsta connectpoints"
+n/plot nt.vcdz%z $cut ! "nsta connectpoints"
+
+zone
+n/plot nt.ot%z $cut ! "nsta connectpoints"
+n/plot nt.om%z $cut ! "nsta connectpoints same blue"
+n/plot nt.or%z $cut ! "nsta connectpoints same red"
+n/plot nt.ol%z $cut ! "nsta connectpoints same green"
+n/plot nt.ok%z $cut ! "nsta connectpoints same orange"
+
+ */
+
+/* ----------------------------------------------------- */
+void tstinterp(CosmoCalc& univ,double z1,double z2,double dz)
+{
+ cout<<"\nTSTINTERP()"<<endl;
+
+ // On repmplit les donnes completes
+ vector<double> Z,D;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   Z.push_back(z);
+   D.push_back(univ.Dloscom(z));
+ }
+ // On remplit un sous tableau
+ int ninc = 5;
+ vector<double> Dinterp;
+ double zmin = Z[0], zmax = Z[0];
+ for(int i=0;i<Z.size();i+=ninc) {
+   Dinterp.push_back(D[i]);
+   zmax = Z[i];
+ }
+ InterpFunc interp(zmin,zmax,Dinterp);
+ unsigned short ok;
+
+ const int n = 8;
+ char *vname[n] = {"z","d","di","dl","dp","zi","zl","zp"};
+ NTuple nt(n,vname);
+ double xnt[n];
+   
+ for(int i=0;i<Z.size();i++) {
+   if(Z[i]>zmax) break;
+   xnt[0] = Z[i];
+   xnt[1] = D[i];
+   xnt[2] = interp(Z[i]);
+   xnt[3] = interp.Linear(Z[i],ok);
+   xnt[4] = interp.Parab(Z[i],ok);
+   xnt[5] = xnt[6] = xnt[7] = 0.;
+   nt.Fill(xnt);
+ }
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtuniv.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+}
+
+/*
+openppf cmvtuniv.ppf
+
+n/plot nt.d%z ! ! "nsta connectpoints"
+n/plot nt.di%z ! ! "nsta connectpoints same green"
+n/plot nt.dl%z ! ! "nsta connectpoints same red"
+n/plot nt.dp%z ! ! "nsta connectpoints same blue"
+
+n/plot nt.di-d%z ! ! "nsta connectpoints green"
+n/plot nt.dl-d%z ! ! "nsta connectpoints same red"
+n/plot nt.dp-d%z ! ! "nsta connectpoints same blue"
+
+n/plot nt.(di-d)/d%z d>0 ! "nsta connectpoints green"
+n/plot nt.(dl-d)/d%z d>0 ! "nsta connectpoints same red"
+n/plot nt.(dp-d)/d%z d>0 ! "nsta connectpoints same blue"
+
+*/
Index: trunk/Cosmo/SimLSS/cmvtvarspec.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtvarspec.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cmvtvarspec.cc	(revision 3115)
@@ -0,0 +1,198 @@
+#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 "constcosmo.h"
+#include "pkspectrum.h"
+
+// cmvtvarspec [zref,R]
+void usage(void);
+void usage(void)
+{
+  cout<<"cmvtvarspec -z zref -R R -n as,ns -k kmin,kmax -i eps,npt,fracmax"<<endl;
+}
+
+
+int main(int narg,char *arg[])
+{
+ double ob0 = 0.0444356;
+ // -- WMAP
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ // -- ouvert matter only
+ //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ // -- plat matter only
+ //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;
+
+ double ns = 1., as = 1.;
+
+ double zref = 0.5;
+ double R=8./h100;
+
+ double kmin=1e-5,kmax=1000.;
+ int npt = 10000;
+ double eps=1.e-3, fracmax=1e-4;
+
+ char c;
+ while((c = getopt(narg,arg,"hz:R:n:k:i:")) != -1) {
+  switch (c) {
+  case 'z' :
+    sscanf(optarg,"%lf",&zref);
+    break;
+  case 'R' :
+    sscanf(optarg,"%lf",&R);
+    if(R<0.) R = -R/h100;
+    break;
+  case 'n' :
+    sscanf(optarg,"%lf,%lf",&as,&ns);
+    break;
+  case 'k' :
+    sscanf(optarg,"%lf,%lf",&kmin,&kmax);
+    break;
+  case 'i' :
+    sscanf(optarg,"%lf,%d,%lf",&eps,&npt,&fracmax);
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+
+ double Rg=R/sqrt(5.);
+ double lkmin=log10(kmin), lkmax=log10(kmax), dlk=(lkmax-lkmin)/npt;
+
+ cout<<"zref="<<zref<<endl;
+ cout<<"as="<<as<<" ns="<<ns<<endl;
+ cout<<"R="<<R<<", Rg="<<Rg<<endl;
+ cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<")"
+     <<", dlk="<<dlk<<" npt="<<npt<<endl;
+ cout<<"eps="<<eps<<"  fracmax="<<fracmax<<endl;
+
+ //-----------------------------------------------------------------
+ InitialSpectrum pkini(ns,as);
+
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+ //tf.SetNoOscEnv(2);
+
+ GrowthFactor d1(om0,ol0);
+
+ PkSpectrum0 pk0(pkini,tf);
+
+ PkSpectrumZ pkz(pk0,d1,zref);
+
+ //-----------------------------------------------------------------
+ double k1,k2;
+ int rc = 0;
+ //----
+ cout<<endl<<"Filtrage top_hat"<<endl;
+ VarianceSpectrum varpk_th(pkz,0);
+ double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
+ double pkmax_th = varpk_th(kfind_th);
+ cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
+ k1=kmin, k2=kmax;
+ rc = varpk_th.FindLimits(R,pkmax_th*fracmax,k1,k2,eps);
+ cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
+
+ varpk_th.SetInteg(0.01,dlk,-1.,4);
+ cout<<"varpk_th="<<varpk_th.Variance(R,k1,k2)<<endl;
+
+ //----
+ cout<<endl<<"Filtrage gaussien"<<endl;
+ VarianceSpectrum varpk_ga(pkz,1);
+ double kfind_ga = varpk_ga.FindMaximum(Rg,kmin,kmax,eps);
+ double pkmax_ga = varpk_ga(kfind_ga);
+ cout<<"kfind_ga = "<<kfind_ga<<" -> "<<log10(kfind_ga)<<", integrand="<<pkmax_ga<<endl;
+ k1=kmin, k2=kmax;
+ rc = varpk_ga.FindLimits(Rg,pkmax_ga*fracmax,k1,k2,eps);
+ cout<<"limit_ga: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
+
+ varpk_ga.SetInteg(0.01,dlk,-1.,4);
+ cout<<"varpk_ga="<<varpk_ga.Variance(Rg,k1,k2)<<endl;
+ 
+ //----
+  cout<<endl<<"Filtrage 1 (integrale du sppectre)"<<endl;
+ VarianceSpectrum varpk_int(pkz,2);
+ double kfind_int = varpk_int.FindMaximum(Rg,kmin,kmax,eps);
+ double pkmax_int = varpk_int(kfind_int);
+ cout<<"kfind_int = "<<kfind_int<<" -> "<<log10(kfind_int)<<", integrand="<<pkmax_int<<endl;
+ k1=kmin, k2=kmax;
+ rc = varpk_int.FindLimits(Rg,pkmax_int*fracmax,k1,k2,eps);
+ cout<<"limit_int: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
+
+ varpk_int.SetInteg(0.01,dlk,-1.,4);
+ cout<<"varpk_int="<<varpk_int.Variance(Rg,k1,k2)<<endl;
+ 
+//-----------------------------------------------------------------
+ const int n = 2;
+ char *vname[n] = {"k","pk"};
+ NTuple nt(n,vname);
+ double xnt[n];
+
+ for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
+   double k = pow(10.,lk);
+   xnt[0] = k;
+   xnt[1] = pkz(k);
+   nt.Fill(xnt);
+ }
+
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvtvarspec.ppf"; POutPersist pos(tag); tag = "nt"; pos.PutObject(nt,tag);
+
+ return 0;
+}
+
+/*
+openppf cmvtvarspec.ppf
+
+#### Spectre a z
+n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"
+
+# L'approximation pour j1(x) (top-hat) et la gaussienne equivalente
+set l1 1e-6
+set l2 10
+set j1sx pow(3*(sin(x)-x*cos(x))/(x*x*x),2.)
+set j1sxa3 (1.-x*x/5.)
+set j1sxa5 (1.-x*x/5.*(1.-3.*x*x/35.))
+set j1sxa7 (1.-x*x/5.*(1.-3.*x*x/35.*(1.-4.*x*x/81.)))
+func $j1sx    $l1 $l2 1000
+func $j1sxa3  $l1 $l2 1000 "red same"
+func $j1sxa5  $l1 $l2 1000 "blue same"
+func $j1sxa7  $l1 $l2 1000 "orange same"
+set gauss  pow(exp(-(x/sqrt(5.))*(x/sqrt(5.))/2.),2.)
+func $gauss  $l1 $l2 1000 "green same"
+
+func $j1sx-$j1sxa3  $l1 $l2 1000 "red"
+func $j1sx-$j1sxa5  $l1 $l2 1000 "blue same"
+func $j1sx-$j1sxa7  $l1 $l2 1000 "orange same"
+func $j1sx-$gauss  $l1 $l2 1000 "green same"
+addline $l1 0 $l2 0
+
+# la fonction "Pk" a integrer
+zone
+set D (k*k*pk/(2.*M_PI*M_PI))
+n/plot nt.$D%log10(k) ! ! "nsta crossmarker3"
+
+# La fonction a integrer en dk lineaire (top-hat et gaussien)
+set R 10
+set RG ($R/sqrt(5.))
+set x (k*$R)
+set j1sx pow(3*(sin($x)-$x*cos($x))/($x*$x*$x),2.)
+set y (k*$RG)
+set gauss  pow(exp(-$y*$y/2.),2.)
+
+n/plot nt.$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
+n/plot nt.$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"
+
+# La fonction a integrer en d(log(x)) (top-hat et gaussien)
+# c'est ce que fait Variance()
+n/plot nt.k*$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
+n/plot nt.k*$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"
+n/plot nt.k*$D%log10(k) ! ! "nsta crossmarker3 same blue"
+
+*/
Index: trunk/Cosmo/SimLSS/constcosmo.h
===================================================================
--- trunk/Cosmo/SimLSS/constcosmo.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/constcosmo.h	(revision 3115)
@@ -0,0 +1,24 @@
+#ifndef CONSTCOSMO_SEEN
+#define CONSTCOSMO_SEEN
+
+static const double MpctoMeters_Cst = 3.0856e+22; // Mpc in meters 
+static const double YeartoSec_Cst = 31556925.2; // seconds in a year
+static const double SpeedOfLight_Cst = 299792.458; // speed of light "km/sec"
+static const double LightYear_Cst = 9.46073042e+15; // 1 Light-Year in "m"
+static const double G_Newton_Cst = 6.6742e-11; // G_N in SI units "m^3 / kg /s^2"
+static const double h_Planck_Cst = 6.6260693e-34;   // Planck constant SI units "J s"
+static const double k_Boltzman_Cst = 1.3806503e-23;  // Boltzman constant SI "J / K"
+static const double ProtonMass_Cst = 1.67262171e-24; // Proton mass in "g"
+static const double SolarMass_Cst = 1.98844e+33; // Proton mass in "g"
+//            Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
+static const double StefBoltz_Cst = 5.670400e-8; // constante de Stefan-Boltzman in "W/m^2/K^4"
+static const double Zeta_3_Cst = 1.2020569031595942854; // Zeta(3)
+
+static const double Jansky2Watt_cst = 1.e-26;  // conversion Jansky to "Watt/m^2/Hz == J/m^2"
+static const double Deg2Rad_Cst = (M_PI/180.); // conversion deg to rad
+
+static const double T_CMB_Par = 2.725;  // temperature du CMB en "K"
+static const double T_NU_Par = 1.9;  // temperature des neutrinos en "K"
+static const double Fr_HyperFin_Par = 1.4204;  // frequence de la raie hyperfine HI en "GHz"
+
+#endif
Index: trunk/Cosmo/SimLSS/cosmocalc.cc
===================================================================
--- trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3115)
@@ -0,0 +1,433 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "pexceptions.h"
+#include "cspline.h"
+
+#include "constcosmo.h"
+#include "cosmocalc.h"
+#include "integfunc.h"
+
+
+///////////////////////////////////////////////////////////
+//*********************** CosmoCalc *********************//
+///////////////////////////////////////////////////////////
+
+//----------------------------------------------------------
+// flat = 0 : no compute Otot from other
+//      = 1 : change Omatter to get Otot=1
+//      = 2 : change Olambda to get Otot=1
+//      = 3 : change Orelat to get Otot=1
+// usespline = "true" : compute integrale then use spline to extrapolate
+//           = "false" : compute integrale everytime
+// zmax : maximum redshift(only if usespline=true)
+         
+CosmoCalc::CosmoCalc(unsigned short flat,bool usespline,double zmax)
+  : _flat(flat),
+    _zold(-1.) , _integval(0.),
+    _usespline(usespline), _computespl(true), _zmax(zmax), _spline(NULL),
+    _nspl(0), _xspl(NULL), _yspl(NULL)
+{
+  if(_usespline && _zmax<=0.) {
+    cout<<"CosmoCalc::CosmoCalc: Error bad _zmax value    zmax="<< _zmax<<endl;
+    throw ParmError("CosmoCalc::UseSpline: Error bad _zmax value");
+  }
+  DefaultParam();
+  SetInteg();
+}
+
+CosmoCalc::CosmoCalc(CosmoCalc& univ)
+  : _spline(NULL), _xspl(NULL), _yspl(NULL)
+{
+  Clone(univ);
+  if(_usespline && _zmax<=0.) {
+    cout<<"CosmoCalc::CosmoCalc: Error bad _zmax value    zmax="<< _zmax<<endl;
+    throw ParmError("CosmoCalc::UseSpline: Error bad _zmax value");
+  }
+}
+
+CosmoCalc::~CosmoCalc(void)
+{
+  if(_spline != NULL) delete _spline; _spline = NULL;
+  if(_xspl != NULL) delete [] _xspl; _xspl = NULL;
+  if(_yspl != NULL) delete [] _yspl; _yspl = NULL;
+}
+
+void CosmoCalc::Clone(CosmoCalc& univ)
+{
+ _flat = univ._flat;
+ _h100 = univ._h100;
+ _H0 = univ._H0;
+ _Olambda0 = univ._Olambda0;
+ _W0 = univ._W0;
+ _Omatter0 = univ._Omatter0;
+ _Obaryon0 = univ._Obaryon0;
+ _Orelat0 = univ._Orelat0;
+ _Otot0 = univ._Otot0;
+ _Ocurv0 = univ._Ocurv0;
+ _Dhubble = univ._Dhubble;
+
+ _dperc = univ._dperc;
+ _dzinc = univ._dzinc;
+ _dzmax = univ._dzmax;
+ _glorder = univ._glorder;
+ if(_xgausl.size()>0) {
+   _xgausl.resize(0); _wgausl.resize(0);
+   for(int i=0;i<(int)_xgausl.size();i++) {
+     _xgausl.push_back(univ._xgausl[i]);
+     _wgausl.push_back(univ._wgausl[i]);
+   }
+ }
+
+ _zold = -1.;
+ _integval = -1.;
+
+ _usespline = univ._usespline;
+ _computespl = true;
+ _zmax = univ._zmax;
+ _nspl = 0;
+ if(_spline) delete _spline; _spline = NULL;
+ if(_xspl) delete [] _xspl; _xspl = NULL;
+ if(_yspl) delete [] _yspl; _yspl = NULL;
+}
+
+//----------------------------------------------------------
+// On va couper l'intervalle entre [0,zmax]:
+// On parcourt [0,zmax] par pas de dzinc
+// et on cree un intervalle
+//     - si la fonction varie de plus de "dperc"
+//     - ou si l'increment en z est superieur a "dzmax"
+void CosmoCalc::SetInteg(double dperc,double dzinc,double dzmax,unsigned short order)
+{
+  _dperc = fabs(dperc);
+
+  _dzinc = fabs(dzinc);
+  if(_dzinc>_zmax/5.) _dzinc = _zmax/5.; // Protection against big dzinc
+
+  _dzmax = fabs(dzmax);
+  if(_dzmax<_dzinc) _dzmax = _dzinc;
+
+  _glorder = order;
+  if(_glorder<=1) _glorder = 4;
+  Compute_GaussLeg(_glorder,_xgausl,_wgausl,0.,1.);
+
+  _computespl=true;
+}
+
+void CosmoCalc::SetObaryon0(double v)
+{
+  _Obaryon0 = v;
+  if(_Obaryon0<0.) {
+    cout<<"CosmoCalc::SetObaryon0: Error bad _Obaryon0 value: "<<_Obaryon0<<endl;
+    throw ParmError("CosmoCalc::SetObaryon0: Error bad _Obaryon0 value");
+  }
+}
+
+void CosmoCalc::SetDynParam(double h100,double om0,double or0,double ol0,double w0)
+{
+  _computespl=true;
+
+  _h100 = h100;
+  _H0 = 100. * _h100 ;
+  _Dhubble = SpeedOfLight_Cst / _H0;
+
+  _Omatter0 = om0;
+  _Orelat0 = or0;
+  _Olambda0 = ol0;
+  _W0 = w0;
+  _Otot0 = _Olambda0 + _Omatter0 + _Orelat0;
+  _Ocurv0 = 1. - _Otot0;
+  if( _h100<0. || _Omatter0<0. || _Orelat0<0. || _Olambda0<0. ) {
+    cout<<"CosmoCalc::SetDynParam: Error bad parameter value"<<endl;
+    throw ParmError("CosmoCalc::SetDynParam: Error bad parameter value");
+  }
+  if(_flat==0) return;
+
+  _Otot0 = 1.; _Ocurv0 = 0.;
+  if(_flat==1) {
+    _Omatter0 = _Otot0 - _Olambda0 - _Orelat0;
+    if( _Omatter0<0. ) {
+      cout<<"CosmoCalc::SetDynParam: Error bad _Omatter0 value: "<<_Omatter0<<endl;
+      throw ParmError("CosmoCalc::SetDynParam: Error bad _Omatter0 value");
+    }
+    return;
+  }
+
+  //--
+  if(_flat==2) {
+    _Olambda0 = _Otot0 - _Omatter0 - _Orelat0;
+    if( _Olambda0<0. ) {
+      cout<<"CosmoCalc::SetDynParam: Error bad _Olambda0 value: "<<_Olambda0<<endl;
+      throw ParmError("CosmoCalc::SetDynParam: Error bad _Olambda0 value");
+    }
+    return;
+  }
+
+  //--
+  if(_flat==3) {
+    _Orelat0 = _Otot0 - _Omatter0 - _Olambda0;
+    if( _Orelat0<0. ) {
+      cout<<"CosmoCalc::SetDynParam: Error bad _Orelat0 value: "<<_Orelat0<<endl;
+      throw ParmError("CosmoCalc::SetDynParam: Error bad _Orelat0 value");
+    }
+    return;
+  }
+
+  cout<<"CosmoCalc::SetDynParam: Error bad _flat value: "<<_flat<<endl;
+  throw ParmError("CosmoCalc::SetDynParam: Error bad _flat value");
+
+}
+
+void CosmoCalc::DefaultParam(void)
+{
+  _computespl=true;
+
+  double h100 = 0.71;
+  double ol0 = 0.73, w0 = -1.;
+  double om0 = 0.135/(h100*h100);
+  // Relat = photons (2.725 K) + neutrinos (1.9 K)
+  double or0 = 2.47e-5 * (1. + 21./8.*pow(T_NU_Par/T_CMB_Par,4.)) / (h100*h100);
+  SetDynParam(h100,om0,or0,ol0,w0);
+  _Obaryon0 = 0.0224/(h100*h100);
+
+  return;
+}
+
+//----------------------------------------------------------
+void CosmoCalc::Print(double z)
+{
+  if(z<0.) z = 0.;
+
+  printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g\n",int(_usespline),ZMax(),Flat(),z);
+  printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3\n"
+        ,h100(),H0(),Dhubble(),H(z),Rhoc(z));
+  printf("Olambda=%g  W0=%g\n",Olambda(z),_W0);
+  printf("Omatter=%g  Obaryon=%g\n",Omatter(z),Obaryon(z));
+  printf("Orelat=%g\n",Orelat(z));
+  printf("Otot=%g  Ocurv=%g\n",Otot(z),Ocurv(z));
+  if(z <= 0.) return;
+  printf("Distance comoving: los=%g Mpc, transv=%g Mpc\n",Dloscom(z),Dtrcom(z));
+  printf("Distance: angular=%g Mpc, lum=%g Mpc\n",Dang(z),Dlum(z));
+  printf("Volume comoving element: %g Mpc^3/sr/dz=1\n",dVol(z));
+  printf("Volume comoving in [0,z] for 4Pi sr: %g Mpc^3\n",Vol4Pi(z));
+
+}
+//----------------------------------------------------------
+double CosmoCalc::Olambda(double z)
+{
+  double v = _Olambda0;
+  if(_W0 != -1.) v *= pow((1.+z),3.*(1+_W0));
+  return v / E2(z);
+}
+
+double CosmoCalc::Omatter(double z)
+{
+  return _Omatter0 * (1.+z)*(1.+z)*(1.+z) /E2(z);
+}
+
+double CosmoCalc::Obaryon(double z)
+{
+  return _Obaryon0 * (1.+z)*(1.+z)*(1.+z) /E2(z);
+}
+
+double CosmoCalc::Orelat(double z)
+{
+  double z2 = (1.+z)*(1.+z);
+  return _Orelat0 * z2*z2 /E2(z);
+}
+
+double CosmoCalc::Ocurv(double z)
+{
+  return _Ocurv0 * (1.+z)*(1.+z) /E2(z);
+}
+
+double CosmoCalc::Otot(double z)
+{
+  return Olambda(z)+Omatter(z)+Orelat(z);
+}
+
+double CosmoCalc::Rhoc(double z)
+// Densite critique au redshift "z"  en "g/cm^3"
+{
+  double h2 = H(z) / MpctoMeters_Cst;  h2 *= h2;
+  return 3.* h2 / (8.*M_PI*G_Newton_Cst) * 1000.;
+}
+
+//----------------------------------------------------------
+double CosmoCalc::Dloscom(double z)     /* Mpc */
+{
+  return _Dhubble * NInteg(z);
+}
+
+double CosmoCalc::Dtrcom(double z)     /* Mpc */
+{
+  double v = NInteg(z);  // Zero curvature
+
+  if(_flat) return _Dhubble * v;
+
+  if(_Ocurv0>0.) {  // hyperbolique
+    double s = sqrt(_Ocurv0);
+    v = sinh(s*v)/s;
+  } else if(_Ocurv0<0.) {  // spherique
+    double s = sqrt(-_Ocurv0);
+    v = sin(s*v)/s;
+  }
+
+  return _Dhubble * v;
+
+}
+
+double CosmoCalc::Dang(double z)     /* Mpc */
+{
+  return Dtrcom(z) / (1.+z);
+}
+
+double CosmoCalc::Dlum(double z)     /* Mpc */
+{
+  return (1.+z) * Dtrcom(z);
+}
+
+double CosmoCalc::dVol(double z)     /* Mpc^3/ sr / unite z */
+{
+  double d = Dtrcom(z);
+  return _Dhubble * d*d / E(z);
+}
+
+//----------------------------------------------------------
+double CosmoCalc::Vol4Pi(double z)  /* Mpc^3 pour 4Pi sr entre [0,z] */
+  // --- on pose x = dm/dh = (1+z)*da/dh, s = sqrt(|Ok|)
+  // Ok=0 : V(z)/dh^3 = (4*Pi/3) * x^3
+  // Ok>0 : V(z)/dh^3 = (4*Pi/(2*Ok)) * (x*sqrt(1+Ok*x^2) - 1/s * asinh(s*x))
+  // Ok<0 : V(z)/dh^3 = (4*Pi/(2*Ok)) * (x*sqrt(1+Ok*x^2) - 1/s * asin(s*x) )
+  // --- on pose y = s*x  (y>0)
+  // Ok>0 : V(z)/dh^3 = (4*Pi/(2*Ok*s)) * (y*sqrt(1+y^2) - asinh(y))
+  // Ok<0 : V(z)/dh^3 = (4*Pi/(2*Ok*s)) * (y*sqrt(1-y^2) - asin(y) )
+  // --- a petit "z" on a pour "y -> 0+" en faisant le DL
+  // Ok>0 : V(z)/dh^3 =  4*Pi*y^3 / (3*Ok*s) * (1 - 3*y^2/10) + O(y^7)
+  // Ok<0 : V(z)/dh^3 = -4*Pi*y^3 / (3*Ok*s) * (1 + 3*y^2/10) + O(y^7)
+  //                          (remarque: Ok^(3/2) = s*Ok)
+{
+  double v,x = Dtrcom(z)/_Dhubble;
+
+  if(_flat) {
+    v = 4.*M_PI/3. * x*x*x;
+  } else if(_Ocurv0>0.) {
+    double s = sqrt(_Ocurv0), y = s*x, y2 = y*y;
+    if(y<1e-6) {
+      v = 1. - 3.*y2/10.;
+      v *= 4.*M_PI*y*y2 / (3.*_Ocurv0*s);
+    } else {
+      v = y*sqrt(1.+y2) - asinh(y);
+      v *= 4.*M_PI/(2.*_Ocurv0*s);
+    }
+  } else if (_Ocurv0<0.) {
+    double s = sqrt(-_Ocurv0), y = s*x, y2 = y*y;
+    if(y<1e-6) {
+      v = 1. + 3.*y2/10.;
+      v *= -4.*M_PI*y*y2 / (3.*_Ocurv0*s);
+    } else {
+      v = y*sqrt(1.-y2) - asin(y);
+      v *= 4.*M_PI/(2.*_Ocurv0*s);
+    }
+  } else {
+    v = 4.*M_PI/3. * x*x*x;
+  }
+
+  return v * _Dhubble*_Dhubble*_Dhubble;
+}
+
+double CosmoCalc::Vol4Pi(double z1,double z2)  /* Mpc^3 pour 4Pi sr entre [z1,z2] */
+{
+  return Vol4Pi(z2) - Vol4Pi(z1);
+}
+
+double CosmoCalc::E2(double z) const
+{
+  z += 1.;
+  double oldum = _Olambda0;
+
+  if(oldum>0. && _W0!=-1.) oldum *= pow(z,3.*(1.+_W0));
+
+  return oldum + z*z*(_Ocurv0 + z*(_Omatter0+z*_Orelat0));
+}
+
+//----------------------------------------------------------
+double CosmoCalc::NInteg(double z)
+{
+  if(z<0.) {
+    cout<<"CosmoCalc::NInteg: Error bad z value  z="<<z<<endl;
+    throw ParmError("CosmoCalc::NInteg: Error bad z value");
+  }
+
+  // On utilise le spline
+  if(_usespline) {
+    if( (!_computespl) && (z==_zold) ) return _integval;
+    if(z>_zmax) {_zmax = z; _computespl = true;}
+    if(_computespl) Init_Spline();
+    _integval = _spline->CSplineInt(z);
+    _zold = z;
+    return _integval;
+  }
+
+  // On calcule l'integrale
+  if(z != _zold) {
+    _integval = IntegrateFunc(*this,0.,z,_dperc,_dzinc,_dzmax,_glorder);
+    _zold = z;
+  }
+  return _integval;
+
+}
+
+int_4 CosmoCalc::Init_Spline(void)
+{
+  if(_spline!=NULL)
+    {delete _spline; delete [] _xspl; delete [] _yspl;}
+
+  //-- Look for intervalles and integrate
+  vector<double> x,y;
+  double zbas=0., fbas = Integrand(zbas);
+  x.push_back(zbas); y.push_back(0.);
+  for(double z=_dzinc;z<_zmax+_dzinc && zbas<_zmax;z+=_dzinc) {
+    double f = Integrand(z);
+    bool doit = false;
+    //cout<<fbas<<","<<f<<"  :  "<<fabs(f-fbas)<<" >? "<<_dperc*fabs(fbas)<<endl;
+    if( fabs(f-fbas)>_dperc*fabs(fbas) ) doit = true;
+    if( z-zbas>_dzmax ) doit = true;
+    if( z>_zmax ) {doit = true; z=_zmax;}
+    if( ! doit ) continue;
+    x.push_back(z);
+    double sum = 0., dz = z-zbas;
+    for(unsigned short i=0;i<_glorder;i++) sum += _wgausl[i]*Integrand(zbas+_xgausl[i]*dz);
+    y.push_back(sum*dz);
+    //cout<<"...set... "<<zbas<<","<<z<<"  ,  "<<fbas<<","<<f<<endl;
+    zbas = z;
+    fbas = f;
+  }
+
+  //-- Fill spline
+  _nspl = x.size();
+  _xspl = new double[_nspl];
+  _yspl = new double[_nspl];
+  for(int i=0;i<_nspl;i++) {
+    _xspl[i] = x[i];
+    _yspl[i] = y[i];
+    if(i!=0) _yspl[i] += _yspl[i-1];
+  }
+  double yp1 = Integrand(_xspl[0]);
+  double ypn = Integrand(_xspl[_nspl-1]);
+  _spline = new CSpline(_nspl,_xspl,_yspl,yp1,ypn,CSpline::NoNatural,false);
+  _spline->ComputeCSpline();
+  _computespl = false;
+
+  cout<<"CosmoCalc::Init_Spline called: (zmax="<<_zmax<<") _nspl="<<_nspl<<endl;
+  for(int i=0;i<(10<_nspl?10:_nspl);i++) cout<<_xspl[i]<<" ";
+  cout<<" ... "<<_xspl[_nspl-1]<<endl;
+
+  return _nspl;
+}
Index: trunk/Cosmo/SimLSS/cosmocalc.h
===================================================================
--- trunk/Cosmo/SimLSS/cosmocalc.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/cosmocalc.h	(revision 3115)
@@ -0,0 +1,99 @@
+#ifndef COSMOCALC_SEEN
+#define COSMOCALC_SEEN
+
+#include "machdefs.h"
+#include <vector>
+
+#include "genericfunc.h"
+
+namespace SOPHYA {
+
+  class CSpline;
+
+//-----------------------------------------------------------------------------------
+class CosmoCalc : public GenericFunc {
+public:
+  CosmoCalc(unsigned short flat=0,bool usespline=false,double zmax=0.);
+  CosmoCalc(CosmoCalc& univ);
+  virtual ~CosmoCalc(void);
+
+  unsigned short Flat(void) { return _flat;}
+
+  double ZMax(void) { return _zmax;}
+  void SetInteg(double dperc=0.01,double dzinc=0.05,double dzmax=5.,unsigned short order=4);
+
+  void SetDynParam(double h100,double om0,double or0,double ol0,double w0=-1.);
+  void SetObaryon0(double v);
+  void DefaultParam(void);
+
+  inline double h100(void) {return _h100;}
+  inline double H0(void) {return _H0;}
+  inline double Dhubble(void) {return _Dhubble;}
+  inline double Olambda0(void) {return _Olambda0;}
+  inline double W0(void) {return _W0;}
+  inline double Omatter0(void) {return _Omatter0;}
+  inline double Obaryon0(void) {return _Obaryon0;}
+  inline double Orelat0(void) {return _Orelat0;}
+  inline double Ocurv0(void) {return _Ocurv0;}
+  inline double Otot0(void) {return _Otot0;}
+  double Rhoc(double z);
+
+  inline double H(double z) {return _H0 * E(z);}
+  double Olambda(double z);
+  double Omatter(double z);
+  double Obaryon(double z);
+  double Orelat(double z);
+  double Ocurv(double z);
+  double Otot(double z);
+
+  double Dtrcom(double z);
+  double Dloscom(double z);
+  double Dang(double z);
+  double Dlum(double z);
+
+  double dVol(double z);
+  double Vol4Pi(double z);
+  double Vol4Pi(double z1,double z2);
+
+  void Print(double z=0.);
+
+  double E2(double z) const;
+  inline double E(double z) const {return sqrt(E2(z));}
+
+  double NInteg(double z);
+  inline double Integrand(double z) const {return 1./E(z);}
+
+  virtual double operator() (double z) const {return Integrand(z);}
+protected:
+  int_4 Init_Spline(void);
+  void Clone(CosmoCalc& univ);
+
+  // La cosmology
+  unsigned short _flat;
+  double _h100, _H0;  /* km/s/Mpc */
+  double _Olambda0,_W0;
+  double _Omatter0,_Obaryon0;
+  double _Orelat0;
+  double _Otot0,_Ocurv0;
+  double _Dhubble;  /* Mpc */
+
+  // Parametres d'integration
+  double _dperc,_dzinc,_dzmax;
+  unsigned short _glorder;
+  vector<double> _xgausl, _wgausl;
+
+  // Pour le calcul direct
+  double _zold,_integval;
+
+  // Pour les splines
+  bool _usespline,_computespl;
+  double _zmax;
+  CSpline * _spline;
+  int_4 _nspl;
+  double *_xspl, *_yspl;
+
+};
+
+} // Fin du namespace
+
+#endif
Index: trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3115)
@@ -0,0 +1,684 @@
+#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 "tarray.h"
+#include "pexceptions.h"
+#include "perandom.h"
+#include "srandgen.h"
+
+#include "constcosmo.h"
+#include "integfunc.h"
+#include "geneutils.h"
+
+#include "genefluct3d.h"
+
+//#define FFTW_THREAD
+
+#define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
+
+//-------------------------------------------------------
+GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz)
+  : T_(T) , pkz_(pkz)
+{
+ SetNThread();
+ SetSize(1,-1,1,1.,-1.,1.);
+}
+
+GeneFluct3D::~GeneFluct3D(void)
+{
+ fftw_destroy_plan(pf_);
+ fftw_destroy_plan(pb_);
+#ifdef FFTW_THREAD
+ if(nthread_>0) fftw_cleanup_threads();
+#endif
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::SetSize(int nx,int ny,int nz,double dx,double dy,double dz)
+{
+ if(nx<=0 || dx<=0. ) {
+   cout<<"GeneFluct3D::SetSize_Error: bad value(s)"<<endl;
+   throw ParmError("GeneFluct3D::SetSize_Error: bad value(s)");
+ }
+
+ Tcontent_ = 0;
+
+ // Les taille des tableaux
+ Nx_ = nx;
+ Ny_ = ny;  if(Ny_ <= 0) Ny_ = Nx_;
+ Nz_ = nz;  if(Nz_ <= 0) Nz_ = Nx_;
+ NRtot_ = Nx_*Ny_*Nz_; // nombre de pixels dans le survey
+ NCz_ =  Nz_/2 +1;
+ NTz_ = 2*NCz_;
+ SzK_[2] = Nx_; SzK_[1] = Ny_; SzK_[0] = NCz_; // a l'envers
+
+ // le pas dans l'espace (Mpc)
+ Dx_ = dx;
+ Dy_ = dy; if(Dy_ <= 0.) Dy_ = Dx_;
+ Dz_ = dz; if(Dz_ <= 0.) Dz_ = Dx_;
+ dVol_ = Dx_*Dy_*Dz_;
+ Vol_ = (Nx_*Dx_)*(Ny_*Dy_)*(Nz_*Dz_);
+
+ // Le pas dans l'espace de Fourier (Mpc^-1)
+ Dkx_ = 2.*M_PI/(Nx_*Dx_);
+ Dky_ = 2.*M_PI/(Ny_*Dy_);
+ Dkz_ = 2.*M_PI/(Nz_*Dz_);
+ Dk3_ = Dkx_*Dky_*Dkz_;
+ 
+ // La frequence de Nyquist en k (Mpc^-1)
+ Knyqx_ = M_PI/Dx_;
+ Knyqy_ = M_PI/Dy_;
+ Knyqz_ = M_PI/Dz_;
+
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::Print(void)
+{
+ cout<<"GeneFluct3D(T_typ="<<Tcontent_<<"): z="<<pkz_.GetZ()<<endl;
+ cout<<"Space Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<Nz_<<" ("<<NTz_<<")  size="
+     <<NRtot_<<endl;
+ cout<<"      Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
+     <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
+ cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<endl;
+ cout<<"        Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
+     <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
+ cout<<"          (2Pi/k: "<<2.*M_PI/Dkx_<<" "<<2.*M_PI/Dky_<<" "<<2.*M_PI/Dkz_<<" Mpc)"<<endl;
+ cout<<"      Nyquist: kx="<<Knyqx_<<" ky="<<Knyqy_<<" kz="<<Knyqz_<<" Mpc^-1"
+     <<", Kmax="<<GetKmax()<<" Mpc^-1"<<endl;
+ cout<<"          (2Pi/k: "<<2.*M_PI/Knyqx_<<" "<<2.*M_PI/Knyqy_<<" "<<2.*M_PI/Knyqz_<<" Mpc)"<<endl;
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::ComputeFourier0(void)
+// cf ComputeFourier() mais avec autre methode de realisation du spectre
+//    (attention on fait une fft pour realiser le spectre)
+{
+ int lp=2;
+
+ T_.ReSize(3,SzK_);
+ T_.SetMemoryMapping(BaseArray::CMemoryMapping);
+
+ // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
+ if(lp>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
+ fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
+ double *data = (double *) (&T_(0,0,0));
+#ifdef FFTW_THREAD
+ if(nthread_>0) {
+   cout<<"Computing with "<<nthread_<<" threads"<<endl;
+   fftw_init_threads();
+   fftw_plan_with_nthreads(nthread_);
+ }
+#endif
+ pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data,fdata,FFTW_ESTIMATE);
+ pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata,data,FFTW_ESTIMATE);
+ if(lp>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
+
+ // --- realisation d'un tableau de tirage gaussiens
+ if(lp>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
+ // On tient compte du pb de normalisation de FFTW3
+ double sntot = sqrt((double)NRtot_);
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   data[ip] = NorRand()/sntot;
+ }
+ if(lp>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
+
+ // --- realisation d'un tableau de tirage gaussiens
+ if(lp>1) PrtTim("--- ComputeFourier0: before fft real ---");
+ fftw_execute(pf_);
+ if(lp>1) PrtTim("--- ComputeFourier0: after fft real ---");
+
+ // --- On remplit avec une realisation
+ if(lp>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
+ T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
+ int lmod = Nx_/10; if(lmod<1) lmod=1;
+ for(int i=0;i<Nx_;i++) {
+   int ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;  kx *= kx;
+   if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   for(int j=0;j<Ny_;j++) {
+     int jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     for(int l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;  kz *= kz;
+       if(i==0 && j==0 && l==0) continue; // Suppression du continu
+       double k = sqrt(kx+ky+kz);
+       // cf normalisation: Peacock, Cosmology, formule 16.38 p504
+       double pk = pkz_(k)/Vol_;
+       // ici pas de "/2" a cause de la remarque ci-dessus
+       T_(l,j,i) *= sqrt(pk);
+     }
+   }
+ }
+ if(lp>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
+
+ double p = compute_power_carte();
+ cout<<"Puissance dans la realisation: "<<p<<endl;
+ if(lp>1) PrtTim("--- ComputeFourier0: after Computing power ---");
+
+ Tcontent_ = 1;
+
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::ComputeFourier(void)
+// Calcule une realisation du spectre "pkz_"
+// Attention: dans TArray le premier indice varie le + vite
+// Explication normalisation: see Coles & Lucchin, Cosmology, p264-265
+// FFTW3: on note N=Nx*Ny*Nz
+// f  --(FFT)-->  F = TF(f)  --(FFT^-1)-->  fb = TF^-1(F) = TF^-1(TF(f))
+// sum(f(x_i)^2) = S
+//                sum(F(nu_i)^2) = S*N
+//                                          sum(fb(x_i)^2) = S*N^2
+{
+ int lp=2;
+ 
+ // --- Dimensionnement du tableau
+ // ATTENTION: TArray adresse en memoire a l'envers du C
+ //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
+ T_.ReSize(3,SzK_);
+ T_.SetMemoryMapping(BaseArray::CMemoryMapping);
+
+ // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
+ if(lp>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
+ fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
+ double *data = (double *) (&T_(0,0,0));
+#ifdef FFTW_THREAD
+ if(nthread_>0) {
+   cout<<"Computing with "<<nthread_<<" threads"<<endl;
+   fftw_init_threads();
+   fftw_plan_with_nthreads(nthread_);
+ }
+#endif
+ pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data,fdata,FFTW_ESTIMATE);
+ pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata,data,FFTW_ESTIMATE);
+ //fftw_print_plan(pb_); cout<<endl;
+ if(lp>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
+
+ // --- RaZ du tableau
+ T_ = complex<r_8>(0.);
+
+ // --- On remplit avec une realisation
+ if(lp>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
+ int lmod = Nx_/10; if(lmod<1) lmod=1;
+ for(int i=0;i<Nx_;i++) {
+   int ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;  kx *= kx;
+   if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   for(int j=0;j<Ny_;j++) {
+     int jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     for(int l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;  kz *= kz;
+       if(i==0 && j==0 && l==0) continue; // Suppression du continu
+       double k = sqrt(kx+ky+kz);
+       // cf normalisation: Peacock, Cosmology, formule 16.38 p504
+       double pk = pkz_(k)/Vol_;
+       // Explication de la division par 2: voir perandom.cc
+       // ou egalement Coles & Lucchin, Cosmology formula 13.7.2 p279
+       T_(l,j,i) = ComplexGaussRan(sqrt(pk/2.));
+     }
+   }
+ }
+ if(lp>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
+
+ manage_coefficients();   // gros effet pour les spectres que l'on utilise !
+ if(lp>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
+
+ double p = compute_power_carte();
+ cout<<"Puissance dans la realisation: "<<p<<endl;
+ if(lp>1) PrtTim("--- ComputeFourier: after before Computing power ---");
+
+ Tcontent_ = 1;
+
+}
+
+int GeneFluct3D::manage_coefficients(void)
+// Take into account the real and complexe conjugate coefficients
+// because we want a realization of a real data in real space
+{
+ fftw_complex *fdata = (fftw_complex *) (&T_(0,0,0));
+
+ // 1./ Le Continu et Nyquist sont reels
+ int nreal = 0;
+ for(int kk=0;kk<2;kk++) {
+   int k=0;  // continu
+   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
+   for(int jj=0;jj<2;jj++) {
+     int j=0;
+     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
+     for(int ii=0;ii<2;ii++) {
+       int i=0;
+       if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
+       size_t ip = k+NCz_*(j+Ny_*i);
+       //cout<<"i="<<i<<" j="<<j<<" k="<<k<<" = ("<<fdata[ip][0]<<","<<fdata[ip][1]<<")"<<endl;
+       fdata[ip][1] = 0.; fdata[ip][0] *= M_SQRT2;
+       nreal++;
+     }
+   }
+ }
+ cout<<"Number of forced real number ="<<nreal<<endl;
+
+ // 2./ Les elements complexe conjugues (tous dans le plan k=0,Nyquist)
+
+ // a./ les lignes et colonnes du continu et de nyquist
+ int nconj1 = 0;
+ for(int kk=0;kk<2;kk++) {
+   int k=0;  // continu
+   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
+   for(int jj=0;jj<2;jj++) { // selon j
+     int j=0;
+     if(jj==1) {if( Ny_%2!=0) continue; else j = Ny_/2;}
+     for(int i=1;i<(Nx_+1)/2;i++) {
+       size_t ip = k+NCz_*(j+Ny_*i);
+       size_t ip1 = k+NCz_*(j+Ny_*(Nx_-i));
+       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       nconj1++;
+     }
+   }
+   for(int ii=0;ii<2;ii++) {
+     int i=0;
+     if(ii==1) {if( Nx_%2!=0) continue; else i = Nx_/2;}
+     for(int j=1;j<(Ny_+1)/2;j++) {
+       size_t ip = k+NCz_*(j+Ny_*i);
+       size_t ip1 = k+NCz_*((Ny_-j)+Ny_*i);
+       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       nconj1++;
+     }
+   }
+ }
+ cout<<"Number of forced conjugate on cont+nyq ="<<nconj1<<endl;
+
+ // b./ les lignes et colonnes hors continu et de nyquist
+ int nconj2 = 0;
+ for(int kk=0;kk<2;kk++) {
+   int k=0;  // continu
+   if(kk==1) {if(Nz_%2!=0) continue; else k = Nz_/2;}  // Nyquist
+   for(int j=1;j<(Ny_+1)/2;j++) {
+     if(Ny_%2==0 && j==Ny_/2) continue; // on ne retraite pas nyquist en j
+     for(int i=1;i<Nx_;i++) {
+       if(Nx_%2==0 && i==Nx_/2) continue; // on ne retraite pas nyquist en i
+       size_t ip = k+NCz_*(j+Ny_*i);
+       size_t ip1 = k+NCz_*((Ny_-j)+Ny_*(Nx_-i));
+       fdata[ip1][0] = fdata[ip][0]; fdata[ip1][1] = -fdata[ip][1];
+       nconj2++;
+     }
+   }
+ }
+ cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
+
+ cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
+
+ return nreal+nconj1+nconj2;
+}
+
+double GeneFluct3D::compute_power_carte(void)
+// Calcul la puissance de la realisation du spectre Pk
+{
+ double s2 = 0.;
+ for(int l=0;l<NCz_;l++)
+   for(int j=0;j<Ny_;j++)
+     for(int i=0;i<Nx_;i++) s2 += MODULE2(T_(l,j,i));
+
+ double s20 = 0.;
+ for(int j=0;j<Ny_;j++)
+   for(int i=0;i<Nx_;i++) s20 += MODULE2(T_(0,j,i));
+
+ double s2n = 0.;
+ if(Nz_%2==0)
+   for(int j=0;j<Ny_;j++)
+     for(int i=0;i<Nx_;i++) s2n += MODULE2(T_(NCz_-1,j,i));
+
+ return 2.*s2 -s20 -s2n;
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::FilterByPixel(void)
+// Filtrage par la fonction fenetre du pixel (parallelepipede)
+// TF = 1/(dx*dy*dz)*Int[{0,dx},{0,dy},{0,dz}]
+//                   e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
+// |TF|^2 = 2*(1-cos(k_x*dx)/dx^2/k_x^2 * (idem dy) * (idem dz)
+// et on traite la divergence en K_x = 0:
+// 2*(1-cos(k*d)/d^2/k^2 ~= 1 - (k*d)^2/12 *[1 - (k*d)^2/30]
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- FilterByPixel: before ---");
+
+ for(int i=0;i<Nx_;i++) {
+   int ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_ *Dx_;
+   double pkx2 = pixelfilter(kx); 
+   for(int j=0;j<Ny_;j++) {
+     int jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_ *Dy_;
+     double pky2 =  pixelfilter(ky); 
+     for(int l=0;l<NCz_;l++) {
+       double kz = l*Dkz_ *Dz_;
+       double pkz2 =  pixelfilter(kz); 
+       T_(l,j,i) *= pkx2*pky2*pkz2;
+     }
+   }
+ }
+
+ if(lp>1) PrtTim("--- FilterByPixel: after ---");
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ComputeReal(void)
+// Calcule une realisation dans l'espace reel
+{
+ int lp=2;
+
+ if( Tcontent_==0 ) {
+   cout<<"GeneFluct3D::ComputeReal_Error: empty array"<<endl;
+   throw ParmError("GeneFluct3D::ComputeReal_Error: empty array");
+ }
+
+ // On fait la FFT
+ if(lp>1) PrtTim("--- ComputeReal: before fftw backward ---");
+ fftw_execute(pb_);
+ if(lp>1) PrtTim("--- ComputeReal: after fftw backward ---");
+
+ Tcontent_ = 2;
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ReComputeFourier(void)
+{
+ int lp=2;
+
+ // On fait la FFT
+ if(lp>1) PrtTim("--- ComputeReal: before fftw forward ---");
+ fftw_execute(pf_);
+ // On corrige du pb de la normalisation de FFTW3
+ double v = (double)NRtot_;
+ for(int i=0;i<Nx_;i++)
+   for(int j=0;j<Ny_;j++)
+     for(int l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
+
+ Tcontent_ = 3;
+ if(lp>1) PrtTim("--- ComputeReal: after fftw forward ---");
+}
+
+//-------------------------------------------------------------------
+int GeneFluct3D::ComputeSpectrum(HProf& hp)
+// Compute spectrum from "T" and fill HProf "hp"
+// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
+// cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
+{
+
+ if(hp.NBins()<0) return -1;
+ hp.Zero();
+ hp.SetErrOpt(false);
+
+ // Attention a l'ordre
+ for(int i=0;i<Nx_;i++) {
+   int ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;  kx *= kx;
+   for(int j=0;j<Ny_;j++) {
+     int jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;  ky *= ky;
+     for(int l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;  kz *= kz;
+       double k = sqrt(kx+ky+kz);
+       double pk = MODULE2(T_(l,j,i));
+       hp.Add(k,pk);
+     }
+   }
+ }
+ hp.UpdateHisto();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ hp *= norm;
+
+ return 0;
+}
+
+//-------------------------------------------------------
+size_t GeneFluct3D::VarianceFrReal(double R,double& var)
+// Recompute MASS variance in spherical top-hat (rayon=R)
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- VarianceFrReal: before computation ---");
+
+ double *data = (double *) (&T_(0,0,0));
+ int dnx = int(R/Dx_+0.5); if(dnx<=0) dnx = 1;
+ int dny = int(R/Dy_+0.5); if(dny<=0) dny = 1;
+ int dnz = int(R/Dz_+0.5); if(dnz<=0) dnz = 1;
+ cout<<"dnx="<<dnx<<" dny="<<dny<<" dnz="<<dnz<<endl;
+
+ double sum=0., sum2=0., r2 = R*R; size_t nsum=0;
+
+ for(int i=dnx;i<Nx_-dnx;i+=dnx) {
+   for(int j=dny;j<Ny_-dny;j+=dny) {
+     for(int l=dnz;l<Nz_-dnz;l+=dnz) {
+       double s=0.; size_t n=0;
+       for(int ii=i-dnx;ii<=i+dnx;ii++) {
+         double x = (ii-i)*Dx_; x *= x;
+         for(int jj=j-dny;jj<=j+dny;jj++) {
+           double y = (jj-j)*Dy_; y *= y;
+           for(int ll=l-dnz;ll<=l+dnz;ll++) {
+             double z = (ll-l)*Dz_; z *= z;
+             if(x+y+z>r2) continue;
+             size_t ip = ll+NTz_*(jj+Ny_*ii);
+             s += 1.+data[ip];
+             n++;
+	   }
+	 }
+       }
+       if(n>0) {sum += s; sum2 += s*s; nsum++;}
+       //cout<<i<<","<<j<<","<<l<<" n="<<n<<" s="<<s<<" sum="<<sum<<" sum2="<<sum2<<endl;
+     }
+   }
+ }
+
+ if(nsum<=1) {var=0.; return nsum;}
+
+ sum /= nsum;
+ sum2 = sum2/nsum - sum*sum;
+ if(lp) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+
+ var = sum2/(sum*sum);  // <dM>^2/<M>^2
+ if(lp) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+
+ if(lp>1) PrtTim("--- VarianceFrReal: after computation ---");
+ return nsum;
+}
+
+//-------------------------------------------------------
+size_t GeneFluct3D::NumberOfBad(double vmin,double vmax)
+// number of pixels outside of ]vmin,vmax[ extremites exclues
+//     ->  vmin and vmax are considered as bad
+{
+ double *data = (double *) (&T_(0,0,0));
+
+ size_t nbad = 0;
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   double v = data[ip];
+   if(v<=vmin || v>=vmax) nbad++;
+ }
+
+ return nbad;
+}
+
+size_t GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax)
+// mean,sigma^2 pour pixels avec valeurs ]vmin,vmax[ extremites exclues
+//   -> mean and sigma^2 are NOT computed with pixels values vmin and vmax
+{
+ double *data = (double *) (&T_(0,0,0));
+
+ size_t n = 0;
+ rm = rs2 = 0.;
+
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   double v = data[ip];
+   if(v<=vmin || v>=vmax) continue;
+   rm += v;
+   rs2 += v*v;
+   n++;
+ }
+
+ if(n>1) {
+   rm /= (double)n;
+   rs2 = rs2/(double)n - rm*rm;
+ }
+
+ return n;
+}
+
+size_t GeneFluct3D::SetToVal(double vmin, double vmax,double val0)
+// set to "val0" if out of range ]vmin,vmax[ extremites exclues
+//     ->  vmin and vmax are set to val0
+{
+ double *data = (double *) (&T_(0,0,0));
+
+ size_t nbad = 0;
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   double v = data[ip];
+   if(v<=vmin || v>=vmax) {data[ip] = val0; nbad++;}
+ }
+
+ return nbad;
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::TurnFluct2Mass(void)
+// d_rho/rho -> Mass  (add one!)
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- TurnFluct2Mass: before computation ---");
+ double *data = (double *) (&T_(0,0,0));
+
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   data[ip] += 1.;
+ }
+
+ Tcontent_ = 4;
+ if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
+}
+
+double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
+// do NOT treate negative or nul values
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
+
+ double *data = (double *) (&T_(0,0,0));
+
+ double m,s2;
+ size_t ngood = MeanSigma2(m,s2,0.,1e+200);
+ if(lp) cout<<"TurnMass2MeanNumber: ngood="<<ngood
+            <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
+
+ // On doit mettre n*Vol galaxies dans notre survey
+ // On en met uniquement dans les pixels de masse >0.
+ // On NE met PAS a zero les pixels <0
+ // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
+ //   comme on ne prend que les pixels >0, on doit normaliser
+ //   a la moyenne de <1+d_rho/rho> sur ces pixels
+ //   (rappel sur tout les pixels <1+d_rho/rho>=1)
+ double dn = n_by_mpc3*Vol_/m /(double)ngood;  // nb de gal a mettre ds 1 px
+ if(lp) cout<<"galaxy density move from "
+            <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
+ double sum = 0.;
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   data[ip] *= dn; // par coherence on multiplie aussi les <=0
+   if(data[ip]>0.) sum += data[ip];  // mais on ne les compte pas
+ }
+ if(lp) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
+
+ Tcontent_ = 6;
+ if(lp>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
+ return sum;
+}
+
+double GeneFluct3D::ApplyPoisson(void)
+// do NOT treate negative or nul mass  -> let it as it is
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+
+ double *data = (double *) (&T_(0,0,0));
+
+ double sum = 0.;
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   double v = data[ip];
+   if(v>0.) {
+     unsigned long dn = PoissRandLimit(v,10.);
+     data[ip] = (double)dn;
+     sum += (double)dn;
+   }
+ }
+ if(lp) cout<<sum<<" galaxies put into survey"<<endl;
+ Tcontent_ = 8;
+
+ if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+ return sum;
+}
+
+double GeneFluct3D::TurnNGal2Mass(FunRan& massdist,bool axeslog)
+// do NOT treate negative or nul mass  -> let it as it is
+// INPUT:
+//   massdist : distribution de masse (m*dn/dm)
+//   axeslog = false : retourne la masse
+//           = true  : retourne le log10(mass)
+// RETURN la masse totale
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- TurnNGal2Mass: before computation ---");
+
+  double *data = (double *) (&T_(0,0,0));
+
+ double sum = 0.;
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   double v = data[ip];
+   if(v>0.) {
+     int ngal = int(v+0.1);
+     data[ip] = 0.;
+     for(int i=0;i<ngal;i++) {
+       double m = massdist.RandomInterp();  // massdist.Random();
+       if(axeslog) m = pow(10.,m);
+       data[ip] += m;
+     }
+     sum += data[ip];
+   }
+ }
+ if(lp) cout<<sum<<" MSol HI mass put into survey"<<endl;
+ Tcontent_ = 10;
+
+ if(lp>1) PrtTim("--- TurnNGal2Mass: after computation ---");
+ return sum;
+}
+
+void GeneFluct3D::AddNoise2Real(double snoise)
+// add noise to every pixels (meme les <=0 !)
+{
+ int lp=2;
+ if(lp>1) PrtTim("--- AddNoise2Real: before computation ---");
+
+ double *data = (double *) (&T_(0,0,0));
+
+ for(int i=0;i<Nx_;i++) for(int j=0;j<Ny_;j++) for(int l=0;l<Nz_;l++) {
+   size_t ip = l+NTz_*(j+Ny_*i);
+   data[ip] += snoise*NorRand();
+ }
+ Tcontent_ = 12;
+
+ if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");
+}
Index: trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3115)
@@ -0,0 +1,92 @@
+#ifndef GENEFLUCT3D_SEEN
+#define GENEFLUCT3D_SEEN
+
+#include "machdefs.h"
+#include "genericfunc.h"
+#include "tarray.h"
+#include "hisprof.h"
+#include "perandom.h"
+
+#include <vector>
+
+#include "pkspectrum.h"
+
+#include "FFTW/fftw3.h"
+
+namespace SOPHYA {
+
+//-----------------------------------------------------------------------------------
+class GeneFluct3D {
+public:
+  GeneFluct3D(TArray< complex<r_8 > >& T,PkSpectrumZ& pkz);
+  virtual ~GeneFluct3D(void);
+
+  void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
+  void SetSize(int nx,int ny,int nz,double dx,double dy,double dz);  // Mpc
+
+  inline void SetZ(double z) {pkz_.SetZ(z);}
+  double GetZ(void) {return pkz_.GetZ();}
+  double GetVol(void) {return Vol_;}
+  double GetDVol(void) {return dVol_;}
+  size_t NPix(void) {return NRtot_;}
+  double GetKmax(void) {return sqrt(Knyqx_*Knyqx_+Knyqy_*Knyqy_+Knyqz_*Knyqz_);}
+  vector<r_8> GetKinc(void)
+    {vector<r_8> dk; dk.push_back(Dkx_); dk.push_back(Dky_); dk.push_back(Dkz_); return dk;}
+  vector<r_8> GetKnyq(void)
+    {vector<r_8> kn; kn.push_back(Knyqx_); kn.push_back(Knyqy_); kn.push_back(Knyqz_); return kn;}
+
+  void ComputeFourier0(void);
+  void ComputeFourier(void);
+  void FilterByPixel(void);
+  int  ComputeSpectrum(HProf& hp);
+  void ComputeReal(void);
+  void ReComputeFourier(void);
+
+  size_t VarianceFrReal(double R,double& var);
+  size_t MeanSigma2(double& rm,double& rs2,double vmin=-1.e+150,double vmax=1.e+150);
+  size_t NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
+  size_t SetToVal(double vmin, double vmax,double val0=0.);
+
+  void TurnFluct2Mass(void);
+  double TurnMass2MeanNumber(double n_by_mpc3);
+  double ApplyPoisson(void);
+  double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
+  double TurnMass2Flux(void);
+  void AddNoise2Real(double snoise);
+
+  void Print(void);
+
+protected:
+  int manage_coefficients(void);
+  double compute_power_carte(void);
+  inline double pixelfilter(double kx)
+    {return (kx<0.05) ? 1.-kx*kx/12.*(1.-kx*kx/30.): 2.*(1.-cos(kx))/(kx*kx);}
+
+
+  int Nx_,Ny_,Nz_,SzK_[3];
+  int NCz_,NTz_;
+  size_t NRtot_;
+  double Dx_,Dy_,Dz_;
+  double dVol_, Vol_;
+  double Dkx_,Dky_,Dkz_;
+  double Knyqx_,Knyqy_,Knyqz_;
+  double Dk3_;
+
+  fftw_plan pf_,pb_;
+  unsigned short nthread_;
+
+  //   = 0 : T_ vide
+  //   = 1 : T_ contient une realisation du spectre (espace fourier)
+  //   = 2 : T_ contient une realisation dans l'espace reel
+  //   = 3 : T_ contient la TF d'une realisation dans l'espace reel (d_rho/rho)
+  //   = 4 : T_ contient (1 + d_rho/rho)
+  //   = ... etc ...
+  unsigned short Tcontent_;
+
+  TArray< complex<r_8 > >& T_;
+  PkSpectrumZ& pkz_;
+};
+
+} // Fin du namespace
+
+#endif
Index: trunk/Cosmo/SimLSS/genericfunc.h
===================================================================
--- trunk/Cosmo/SimLSS/genericfunc.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/genericfunc.h	(revision 3115)
@@ -0,0 +1,38 @@
+#ifndef GENERICFUNC_H
+#define GENERICFUNC_H
+
+#include "pexceptions.h"
+
+namespace SOPHYA {
+
+//! Abstract interface definition for functions (double f(double)) viewed as classes
+class GenericFunc {
+public:
+  GenericFunc(void) { } 
+  virtual ~GenericFunc(void) { }
+
+  virtual double operator()(double x) {
+    cout<<"GenericFunc::operator(double) not implemented"<<endl;
+    throw NotAvailableOperation("GenericFunc::operator(double) not implemented");
+  }
+
+  virtual double operator()(double x,double y) {
+    cout<<"GenericFunc::operator(double,double) not implemented"<<endl;
+    throw NotAvailableOperation("GenericFunc::operator(double,double) not implemented");
+  }
+
+  virtual double operator()(double x,double y,double z) {
+    cout<<"GenericFunc::operator(double,double,double) not implemented"<<endl;
+    throw NotAvailableOperation("GenericFunc::operator(double,double,double) not implemented");
+  }
+
+  virtual double operator()(double* x) {
+    cout<<"GenericFunc::operator(double*) not implemented"<<endl;
+    throw NotAvailableOperation("GenericFunc::operator(double*) not implemented");
+  }
+
+};
+
+}
+
+#endif
Index: trunk/Cosmo/SimLSS/geneutils.cc
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/geneutils.cc	(revision 3115)
@@ -0,0 +1,156 @@
+#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 "hisprof.h"
+#include "srandgen.h"
+
+#include "geneutils.h"
+
+//-------------------------------------------------------------------
+// Classe d'interpolation lineaire:
+// Le vecteur y a n elements y_i tels que y_i = f(x_i)
+// Les x_i sont regulierement espaces
+//   et x_0=xmin et x_{n-1}=xmax
+InterpFunc::InterpFunc(double xmin,double xmax,vector<double>& y)
+  : _xmin(xmin), _xmax(xmax), _y(y)
+{
+  if(_xmin>=_xmax || _y.size()<=2) {  // au moins 3 points!
+   cout<<"InterpFunc::InterpFunc : bad arguments values"<<endl;
+   throw ParmError("InterpFunc::InterpFunc : bad arguments values");
+  }
+  _nm1   = _y.size()-1;
+  _xlarg = _xmax-_xmin;
+  _dx    = _xlarg/(double)_nm1;
+}
+
+double InterpFunc::Linear(double x,unsigned short& ok)
+{
+  ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
+  x -= _xmin;
+  long i = long(x/_xlarg*_nm1);  // On prend le "i" juste en dessous
+  if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
+  return _y[i] + (_y[i+1]-_y[i])/_dx*(x-i*_dx);
+}
+
+double InterpFunc::Parab(double x,unsigned short& ok)
+{
+  ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
+  x -= _xmin;
+  long i = long(x/_xlarg*_nm1+0.5);  // On prend le "i" le + proche
+  if(i<1) i=1; else if(i>=_nm1-1) i=_nm1-2;
+  double a = (_y[i+1]-2.*_y[i]+_y[i-1])/(2.*_dx*_dx);
+  double b = (_y[i+1]-_y[i-1])/(2.*_dx);
+  return _y[i] + (x-i*_dx)*(a*(x-i*_dx)+b);
+}
+
+//-------------------------------------------------------------------
+int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex)
+// Remplit l'histo 1D "h" avec la fonction "func"
+// INPUT:
+// logaxex = false : remplissage lineaire
+//        les abscisses "x" des bins sont remplis avec f(x)
+// logaxex = true : remplissage logarithmique (base 10)
+//        les abscisses "x" des bins sont remplis avec f(10^x)
+// RETURN:
+//       0 = OK
+//       1 = error
+{
+  if(h.NBins()<=0) return 1;
+
+  h.Zero();
+
+  for(int_4 i=0;i<h.NBins();i++) {
+    double x = h.BinCenter(i);
+    if(logaxex) x = pow(10.,x);
+    h.SetBin(i,func(x));
+  }
+
+  return 0;
+}
+
+int FuncToVec(GenericFunc& func,TVector<r_8>& v,double xmin,double xmax,bool logaxex)
+// Remplit le TVector avec la fonction "func"
+// INPUT:
+// logaxex = false : remplissage lineaire
+//        les abscisses "x" des bins sont remplis avec f(x)
+// logaxex = true : remplissage logarithmique (base 10)
+//        les abscisses "x" des bins sont remplis avec f(10^x)
+// RETURN:
+//       0 = OK
+//       1 = error
+// Remarque:
+//  v(i) = f(xmin+i*dx) avec dx = (xmax-xmin)/v.NElts()
+{
+  if(v.NElts()<=0 || xmax<=xmin) return 1;
+
+  v = 0.;
+  double dx = (xmax-xmin)/v.NElts();
+   
+  for(int_4 i=0;i<v.NElts();i++) {
+    double x = xmin + i * dx;;
+    if(logaxex) x = pow(10.,x);
+    v(i) = func(x);
+  }
+
+  return 0;
+}
+
+//-------------------------------------------------------------------
+double AngSol(double dtheta,double dphi,double theta0)
+// Retourne l'angle solide d'un "rectangle" et coordonnees spheriques
+// de DEMI-COTE "dtheta" x "dphi" et centre en "theta0"
+// Attention: Le "theta0" de l'equateur est Pi/2  (et non pas zero)
+//            Les unites des angles sont en radians
+//               theta0 in [0,Pi]
+//               dtheta in [0,Pi]
+//               dphi   in [0,2Pi]
+// Return: l'angle solide en steradian
+{
+  double theta1 = theta0-dtheta, theta2 = theta0+dtheta;
+  if(theta1<0.) theta1=0.;
+  if(theta2>M_PI) theta2=M_PI;
+
+  return   2.*dphi * (cos(theta1)-cos(theta2));
+}
+
+double AngSol(double dtheta)
+// Retourne l'angle solide d'une calotte spherique de demi-ouverture "dtheta"
+// Attention: Les unites des angles sont en radians
+//               dtheta in [0,Pi]
+// Return: l'angle solide en steradian
+// Approx pour theta petit: PI theta^2
+{
+  return   2.*M_PI * (1.-cos(dtheta));
+}
+
+//-------------------------------------------------------------------
+unsigned long PoissRandLimit(double mu,double mumax)
+{
+  double pp,ppi,x;
+  unsigned long n;
+
+  if(mu>=mumax) {
+    pp = sqrt(mu);
+    while( (x=pp*NorRand()) < -mu );
+    return (unsigned long)(mu+x+0.5);
+  }
+
+  ppi = pp = exp(-mu);
+  x = drand01();
+  n = 0;
+  while (x > ppi) {
+    n++;
+    pp = mu*pp/(double)n;
+    ppi += pp;
+  }
+  return n;
+}
Index: trunk/Cosmo/SimLSS/geneutils.h
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/geneutils.h	(revision 3115)
@@ -0,0 +1,57 @@
+#ifndef GENEUTILS_SEEN
+#define GENEUTILS_SEEN
+
+#include "machdefs.h"
+#include <math.h>
+#include "genericfunc.h"
+#include "histos.h"
+#include "tvector.h"
+
+#include <vector>
+
+namespace SOPHYA {
+
+class InterpFunc {
+public:
+  InterpFunc(double xmin,double xmax,vector<double>& y);
+  virtual ~InterpFunc(void) { }
+
+  //! Retourne l'element le plus proche de f donnant y=f(x)
+  inline double operator()(double x)
+         {
+         x -= _xmin;
+         long i = long(x/_xlarg*_nm1+0.5);  // On prend le "i" le plus proche
+         if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
+         return _y[i];
+         }
+
+  // idem operator(double) et retourne
+  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
+  inline double operator()(double x,unsigned short& ok)
+    {ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2; return (*this)(x);}
+
+  //! Retourne l'interpolation lineaire de f donnant y=f(x)
+  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
+  double Linear(double x,unsigned short& ok);
+
+  //! Retourne l'interpolation parabolique de f donnant y=f(x)
+  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
+  double Parab(double x,unsigned short& ok);
+
+protected:
+  double _xmin,_xmax,_xlarg,_dx;
+  vector<double>& _y;
+  long _nm1;  // n-1
+};
+
+}
+
+int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex=false);
+int FuncToVec(GenericFunc& func,TVector<r_8>& h,double xmin,double xmax,bool logaxex=false);
+
+double AngSol(double dtheta,double dphi,double theta0=M_PI/2.);
+double AngSol(double dtheta);
+
+unsigned long PoissRandLimit(double mu,double mumax=10.);
+
+#endif
Index: trunk/Cosmo/SimLSS/integfunc.cc
===================================================================
--- trunk/Cosmo/SimLSS/integfunc.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/integfunc.cc	(revision 3115)
@@ -0,0 +1,171 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "pexceptions.h"
+
+#include "integfunc.h"
+
+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/integfunc.h
===================================================================
--- trunk/Cosmo/SimLSS/integfunc.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/integfunc.h	(revision 3115)
@@ -0,0 +1,17 @@
+#ifndef INTEGFUNC_SEEN
+#define INTEGFUNC_SEEN
+
+#include "machdefs.h"
+#include "genericfunc.h"
+#include <vector>
+
+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 3115)
+++ trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3115)
@@ -0,0 +1,565 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "pexceptions.h"
+
+#include "constcosmo.h"
+#include "integfunc.h"
+#include "pkspectrum.h"
+
+
+///////////////////////////////////////////////////////////
+//******************** InitialSpectrum ******************//
+///////////////////////////////////////////////////////////
+
+InitialSpectrum::InitialSpectrum(double n,double a)
+  : n_(n), A_(a)
+{
+}
+
+InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf)
+  : n_(pkinf.n_), A_(pkinf.A_)
+{
+}
+
+InitialSpectrum::~InitialSpectrum(void)
+{
+}
+
+void InitialSpectrum::SetNorm(double a)
+{
+  A_ = a;
+}
+
+void InitialSpectrum::SetSlope(double n)
+{
+  n_ = n;
+}
+
+
+///////////////////////////////////////////////////////////
+//****************** TransfertEisenstein ****************//
+///////////////////////////////////////////////////////////
+
+// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
+TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
+  : Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
+  , nobaryon_(nobaryon) , nooscenv_(0)
+{
+  if(nobaryon_ == false && Ob_>0.) Init_With_Baryons();
+  else Init_Without_Baryon();
+}
+
+TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf)
+  : Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
+  , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_)
+{
+  if(nobaryon_ == false && Ob_>0.) Init_With_Baryons();
+  else Init_Without_Baryon();
+}
+
+void TransfertEisenstein::Init_With_Baryons(void)
+{
+  int lp = 1;
+
+  nobaryon_ = false;
+  O0_ = Oc_ + Ob_;
+    if(lp) cout<<"Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
+
+  double H0 = 100. * h_, h2 = h_*h_;
+
+  th2p7_ = tcmb_/2.7;
+  double th2p7P4 = th2p7_*th2p7_*th2p7_*th2p7_;
+
+  // Formule 2 p 606
+  zeq_ = 2.50e4 * O0_ * h2 / th2p7P4;
+    if(lp) cout<<"zeq = "<<zeq_<<endl;
+
+  // Formule 3 p 607
+  // (attention ici C=1 : H0 -> H0/C si on utilise la premiere formule)
+  //  keq_ = sqrt(2.*O0_*H0*H0*zeq_) / SpeedOfLight_Cst;
+  keq_ = 7.46e-2 * O0_ * h2 / (th2p7_*th2p7_);
+    if(lp) cout<<"keq = "<<keq_<<" Mpc^-1"<<endl;
+
+  // Formule 4 p 607
+  double b1_eq4 = 0.313*pow(O0_*h2,-0.419)*(1. + 0.607*pow(O0_*h2,0.674));
+  double b2_eq4 = 0.238*pow(O0_*h2,0.223);
+  double zd_ = 1291. * pow(O0_*h2,0.251) / (1.+0.659* pow(O0_*h2,0.828))
+             * (1. + b1_eq4*pow(Ob_*h2,b2_eq4));
+    if(lp) cout<<"zd = "<<zd_<<endl;
+
+  // Formule 5 page 607
+  Req_ = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zeq_);
+  Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zd_);
+    if(lp) cout<<"Req = "<<Req_<<" Rd = "<<Rd_<<endl;
+
+  // Formule 6 p 607
+  s_ = 2./(3.*keq_) * sqrt(6./Req_)
+      * log( (sqrt(1.+Rd_) + sqrt(Rd_+Req_)) / (1.+sqrt(Req_)) );
+    if(lp) cout<<"s = "<<s_<<" Mpc"<<endl;
+
+  // Formule 7 page 607
+  ksilk_ = 1.6*pow(Ob_*h2,0.52)*pow(O0_*h2,0.73) * (1. + pow(10.4*O0_*h2,-0.95));
+    if(lp) cout<<"ksilk = "<<ksilk_<<" Mpc^-1"<<endl;
+
+  // Formules 10 page 608
+  double a1 = pow(46.9*O0_*h2,0.670) * (1. + pow(32.1*O0_*h2,-0.532));
+  double a2 = pow(12.0*O0_*h2,0.424) * (1. + pow(45.0*O0_*h2,-0.582));
+  alphac_ = pow(a1,-Ob_/O0_) * pow(a2,-pow(Ob_/O0_,3.));
+  double b1 = 0.944 / (1. + pow(458.*O0_*h2,-0.708));
+  double b2 = pow(0.395*O0_*h2,-0.0266);
+  betac_ = 1 / ( 1. + b1*(pow(Oc_/O0_,b2) - 1.) );
+
+  // Formule 23 page 610
+  bnode_ = 8.41 * pow(O0_*h2,0.435);
+    if(lp) cout<<"bnode = "<<bnode_<<endl;
+
+  // Formule 14 page 608
+  double y = (1.+zeq_)/(1.+zd_);
+  double s1py = sqrt(1.+y);
+  double Gy = y*( -6.*s1py + (2.+3.*y)*log((s1py+1.)/(s1py-1.)) );
+  alphab_ = 2.07*keq_*s_*pow(1.+Rd_,-3./4.)*Gy;
+
+  // Formule 24 page 610
+  betab_ = 0.5 + Ob_/O0_
+         + (3.-2.*Ob_/O0_) * sqrt(pow(17.2*O0_*h2,2.) + 1.);
+
+  // Formule 31 page 612
+  alphag_ = 1.
+          - 0.328*log(431.*O0_*h2)*Ob_/O0_
+          + 0.38*log(22.3*O0_*h2)*pow(Ob_/O0_,2.);
+
+ // --- Pour la positon du premier pic acoustique
+ // Formule 26 page 611
+ double s_peak = 44.5*log(9.83/(O0_*h2)) / sqrt(1.+10.*pow(Ob_*h2,3./4.));  // Mpc
+ // Formule 25 page 611
+ kpeak_ = 5*M_PI/(2.*s_peak) * (1.+0.217*O0_*h2);  // 1/Mpc
+   if(lp) cout<<"for 1er peak: s="<<s_peak<<" kpeak="<<kpeak_<<endl;
+
+  return;
+}
+
+void TransfertEisenstein::Init_Without_Baryon(void)
+{
+  int lp = 1;
+
+  nobaryon_ = true;
+  O0_ = Oc_;
+  if(lp) cout<<"Omatter="<<O0_<<" Ocdm="<<Oc_<<endl;
+  th2p7_ = tcmb_/2.7;
+}
+
+TransfertEisenstein::~TransfertEisenstein(void)
+{
+}
+
+void TransfertEisenstein::SetNoOscEnv(unsigned short nooscenv)
+// To obtain the non-oscillatory part of the spectrum
+// nooscenv = 0 : use the baryon oscillatory part of transfert function
+// nooscenv = 1 : use approx. paragraph 3.3 p610 (middle of right column)
+// nooscenv = 2 : use formulae 29+30+31 page 612
+{
+  if(nooscenv>2) {
+    cout<<"TransfertEisenstein::SetNoOscEnv: Error bad nooscenv value: "<<nooscenv<<endl;
+    throw ParmError("TransfertEisenstein::SetNoOscEnv: Error bad nooscenv value");
+  }
+  nooscenv_ = nooscenv;
+}
+
+double TransfertEisenstein::T0tild(double k,double alphac,double betac)
+{
+  // Formule 10 p 608
+  //double q = k*th2p7_*th2p7_/(O0_*h_*h_);
+  double q = k/(13.41*keq_);
+  // Formule 20 p 610
+  double C = (14.2/alphac) + 386./(1.+69.9*pow(q,1.08));
+  // Formule 19 p 610
+  double x = log(M_E+1.8*betac*q);
+  return x / (x + C*q*q);
+}
+
+double TransfertEisenstein::operator() (double k)
+{
+
+  // --- Pour zero baryon
+  //  OU Pour function lissee sans oscillation baryon
+  if(nobaryon_  || nooscenv_ == 2) {
+    double gamma = O0_*h_;
+    // Formule 30 page 612 (pour fct lissee)
+    if( nobaryon_==false && nooscenv_ == 2 )
+      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*s_,4.)));
+    // Formule 28 page 612
+    double q = k / h_ * th2p7_*th2p7_/gamma;  // Mpc^-1
+    // Formules 29 page 612
+    double l0 = log(2.*M_E + 1.8*q);
+    double c0 = 14.2 + 731./(1.+62.5*q);
+    return l0 / (l0 + c0*q*q);
+  }
+
+  // --- CDM
+  double f = 1. / (1. + pow(k*s_/5.4,4.));
+  double Tc = f*T0tild(k,1.,betac_) + (1.-f)*T0tild(k,alphac_,betac_);
+
+  // --- Baryons
+  // Formule 22 page 610
+  double stk, ksbnode = k*s_/bnode_;
+  if(ksbnode<0.001) stk =s_ * ksbnode;
+     else   stk = s_ / pow(1. + pow(1./ksbnode,3.), 1/.3);
+  // Formule 21 page 610
+  double j0kst = 0.;
+  if(nooscenv_ == 1) j0kst = pow(1.+pow(k*stk,4.) , -1./4.); //lissee sans oscillation baryon
+  else {
+    double x = k*stk;
+    if(x<0.01) j0kst = 1. - x*x/6.*(1.-x*x/20.);
+      else j0kst = sin(x)/x;
+    ////if(k>0.038&&k<0.07) cout<<"k="<<k<<" stk="<<stk<<" x="<<x<<endl;
+  }
+  double Tb = T0tild(k,1.,1.) / (1. + pow(k*s_/5.2,2.));
+  Tb += alphab_/(1.+pow(betab_/(k*s_),3.)) * exp(-pow(k/ksilk_,1/4.));
+  Tb *= j0kst;
+
+  // --- Total
+  double T = (Ob_/O0_)*Tb + (Oc_/O0_)*Tc;
+
+  return T;
+}
+
+double TransfertEisenstein::KPeak(void)
+// Position du premier pic acoustic
+{
+ if(nobaryon_) return -1.;
+ return kpeak_;
+}
+
+
+///////////////////////////////////////////////////////////
+//********************* GrowthFactor ********************//
+///////////////////////////////////////////////////////////
+
+// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
+GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
+  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
+{
+  if(OmegaMatter0==0.) {
+    cout<<"GrowthFactor::GrowthFactor: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
+    throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
+  }
+  norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
+  norm_ = (*this)(0.);
+  cout<<"GrowthFactor::GrowthFactor : norm="<<norm_<<endl;
+}
+
+GrowthFactor::GrowthFactor(GrowthFactor& d1)
+  : O0_(d1.O0_) , Ol_(d1.Ol_) , Ok_(d1.Ok_) , norm_(d1.norm_)
+{
+}
+
+GrowthFactor::~GrowthFactor(void)
+{
+}
+
+double GrowthFactor::operator() (double z)
+// see Formulae A4 + A5 + A6 page 614
+{
+ z += 1.;
+ double z2 = z*z, z3 = z2*z;
+ double den = Ol_ + Ok_*z2 + O0_*z3;
+ double o0z = O0_ *z3 / den;
+ double olz = Ol_ / den;
+
+ // 4./7. = 0.571429
+ double D1z = pow(o0z,0.571429) - olz + (1.+o0z/2.)*(1.+olz/70.);
+ D1z = 2.5*o0z / z / D1z;
+
+ return D1z / norm_;
+}
+
+
+///////////////////////////////////////////////////////////
+//************** PkSpectrum0 et PkSpectrumZ *************//
+///////////////////////////////////////////////////////////
+
+PkSpectrum0::PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf)
+  : pkinf_(pkinf) , tf_(tf)
+{
+}
+
+PkSpectrum0::PkSpectrum0(PkSpectrum0& pk0)
+  : pkinf_(pk0.pkinf_) , tf_(pk0.tf_)
+{
+}
+
+PkSpectrum0::~PkSpectrum0(void)
+{
+}
+
+double PkSpectrum0::operator() (double k)
+{
+  double tf = tf_(k);
+  double pkinf = pkinf_(k);
+  return pkinf *tf*tf;
+}
+
+//------------------------------------
+PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
+  : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(0)
+  , zold_(-1.) , d1old_(1.)
+{
+}
+
+PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
+  : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(0)
+  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
+{
+}
+
+PkSpectrumZ::~PkSpectrumZ(void) 
+{
+}
+
+void PkSpectrumZ::SetTypSpec(unsigned short typspec)
+  // typsec = 0 : compute Pk(k)
+  //        = 1 : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
+{
+  if(typspec>1) {
+    cout<<"PkSpectrumZ::SetTypSpec: Error bad typspec value: "<<typspec<<endl;
+    throw ParmError("PkSpectrumZ::SetTypSpec: Error bad typspec value");
+  }
+  typspec_ = typspec;
+}
+
+double PkSpectrumZ::operator() (double k)
+{
+  return (*this)(k,zref_);
+}
+
+double PkSpectrumZ::operator() (double k,double z)
+{
+  double d1;
+  if(z == zold_) d1 = d1old_;
+    else {d1 = d1old_ = d1_(z); zold_ = z;}
+
+  double v = pk0_(k) * d1*d1;
+  if(typspec_) v *= k*k*k/(2.*M_PI*M_PI);
+
+  return scale_ * v;
+}
+
+
+
+///////////////////////////////////////////////////////////
+//******************* VarianceSpectrum ******************//
+///////////////////////////////////////////////////////////
+
+VarianceSpectrum::VarianceSpectrum(GenericFunc& pk,unsigned short typfilter=0)
+  : pk_(pk) , R_(0.)
+{
+  SetFilter(typfilter);
+}
+
+VarianceSpectrum::VarianceSpectrum(VarianceSpectrum& vpk)
+  : pk_(vpk.pk_) , R_(vpk.R_)
+{
+  SetFilter(vpk.typfilter_);
+}
+
+VarianceSpectrum::~VarianceSpectrum(void)
+{
+}
+
+//------------------------------------
+void VarianceSpectrum::SetFilter(unsigned short typfilter)
+// typfilter = 0 : spherical 3D top-hat
+//           = 1 : spherical 3D gaussian
+//           = 2 : no filter juste integrate spectrum)
+{
+  if(typfilter>2) {
+    cout<<"VarianceSpectrum::SetFilter: Error bad value for type of filter"<<endl;
+    throw ParmError("VarianceSpectrum::SetFilter: Error bad value for type of filter");
+  }
+  typfilter_ = typfilter;
+}
+
+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)
+{
+  dperc_ = dperc;  if(dperc_<=0.) dperc_ = 0.1;
+  dlogkinc_ = dlogkinc;
+  dlogkmax_ = dlogkmax;
+  glorder_ = glorder;
+}
+
+
+//------------------------------------
+double VarianceSpectrum::Filter2(double x)
+// ATTENTION: c'est le filtre au carre qui est renvoye
+{
+  // Just integrate the spectrum without filtering
+  if(typfilter_ == 2) return 1.;
+
+  double x2 = x*x;
+  // Filtre gaussien G(x) = exp(-x^2/2)
+  //        remarque G(x)^2 = exp(-x^2)
+  // on prend le DL de G(x)^2 pour x->0 a l'ordre O(x^6)
+  //             DL(x) = 1-x^2*(1-x^2/2)
+  //             pour x<0.01  |DL(x)-G(X)^2|<2.0e-13
+  if(typfilter_ == 1)
+    if(x<0.01) return 1.-x2*(1.-x2/2.); else return exp(-x2);
+
+  // Filtre top-hat T(x) = 3*(sin(x)-x*cos(x))/x^3
+  // --- Gestion de la pseudo-divergence pour x->0
+  // on prend le DL de T(x)^2 pour x->0 a l'ordre O(x^7)
+  //             DL(x) = 1-x^2/5*(1-3*x^2/35*(1-4*x^2/81))
+  //             pour x<0.1  |DL(x)-T(X)^2|<2.5e-13
+  double f2=0.;
+  if(x<0.1) {
+    f2 = 1.-x2/5.*(1.-3.*x2/35.*(1.-4.*x2/81.));
+  } else {
+    f2 = 3.*(sin(x)-x*cos(x))/(x2*x);
+    f2 *= f2;
+  }
+  return f2;
+  
+}
+
+double VarianceSpectrum::Variance(double R,double kmin,double kmax)
+// Compute variance of spectrum pk_ by integration
+// Input:
+//     R = taille du filter top-hat ou gaussien
+//     kmin,kmax = bornes en k de l'integrale pour calculer la variance
+// Return:
+//     valeur de la variance (sigma^2)
+// Remarque:
+//   la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5))
+//   la variance renvoyee est la variance de la masse
+{
+  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
+    cout<<"VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
+    throw ParmError("VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax");
+  }
+
+  R_ = R;
+  double lkmin = log10(kmin), lkmax = log10(kmax);
+
+  double var = IntegrateFuncLog(*this,lkmin,lkmax,dperc_,dlogkinc_,dlogkmax_,glorder_);
+
+  return var;
+}
+
+//------------------------------------
+double VarianceSpectrum::FindMaximum(double R,double kmin,double kmax,double eps)
+// Retourne le maximum de la fonction a integrer
+// La recherche a lieu entre [kmin,kmax] par pas logarithmiques
+// Input:
+//     R : taille du filter top-hat ou gaussien
+//     kmin,kmax : intervalle de recherche
+//     eps : precision requise sur les valeurs
+// Return:
+//     position (en k) du maximum
+{
+  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
+    cout<<"VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
+    throw ParmError("VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
+  }
+
+  R_ = R;
+
+  int n = 10; // toujours >2
+  double lkmin = log10(kmin), lkmax = log10(kmax), dlk = (lkmax-lkmin)/n;
+
+  double lkfind=lkmin, pkfind=-1.;
+  while(1) {
+    for(int i=0; i<=n; i++) {
+      double lk = lkmin  + i*dlk;
+      double v = (*this)(pow(10.,lk));
+      if(v<pkfind) continue;
+      pkfind = v; lkfind = lk;
+    }
+    //cout<<"VarianceSpectrum::FindMaximum: lkfind="<<lkfind<<" pkfind="<<pkfind
+    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
+    // --- Convergence si l'encadrement de "kfind" est tel que "dk/kfind<eps"
+    // On a dk = 10^(lkfind+dlk) - 10^(lkfind-dlk) = kfind * (10^(dlk) - 10^(-dlk))
+    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
+    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
+    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
+    dlk = (lkmax-lkmin)/n;
+  }
+
+  return pow(10.,lkfind);
+}
+
+int VarianceSpectrum::FindLimits(double R,double high,double &kmin,double &kmax,double eps)
+// Retourne "[kmin,kmax]" tel que la fonction a integrer soit "f(k) <= high"
+//   La recherche a lieu entre [kmin,kmax] par pas logarithmiques
+// Input:
+//     R : taille du filter top-hat ou gaussien
+//     kmin,kmax : intervalle de recherche
+//     eps : precision requise sur les valeurs kmin et kmax
+// Output:
+//     kmin,kmax telles que "f(k) <= high"
+// Return:
+//     rc  = 0 si OK
+//     rc |= 1 "f(kmin) >= high"   (bit0 =1)
+//     rc |= 2 "f(kmax) >= high"   (bit1 =1)
+//     rc |= 4 "f(k) < high pour tout k"   (bit2 =1)
+{
+  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
+    cout<<"VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
+    throw ParmError("VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
+  }
+
+  R_ = R;
+  int n = 10; // toujours >2
+
+  int rc = 0;
+  double lkmin,lkmax,dlk,lkfind;
+
+  // --- Find kmin
+  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
+  while(1) {
+    lkfind = lkmin;
+    for(int i=0;i<=n;i++) {
+      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
+      lkfind = lkmin + i*dlk;
+    }
+    //cout<<"VarianceSpectrum::FindLimits[kmin]: lkfind="<<lkfind
+    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
+    if(fabs(lkfind-lkmax)<dlk/2.) {rc |= 4; return rc;} // protect against f(k)<high for all k
+    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
+    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
+    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
+    dlk = (lkmax-lkmin)/n;
+  }
+  if(lkfind-lkmin<dlk/2.) rc |= 1;  // f(kmin) >= high
+  else kmin = pow(10.,lkmin);
+  //cout<<"rc="<<rc<<" lkmin="<<lkmin<<"  pk="<<(*this)(pow(10.,lkmin))<<endl;
+
+  // --- Find kmax
+  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
+  while(1) {
+    lkfind=lkmax;
+    for(int i=0;i<=n;i++) {
+      if( (*this)(pow(10,lkfind)) >= high ) break;   //cmvbug
+      lkfind -= dlk;
+      lkfind = lkmax - i*dlk;
+    }
+    //cout<<"VarianceSpectrum::FindLimits[kmax]: lkfind="<<lkfind
+    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
+    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
+    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
+    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
+    dlk = (lkmax-lkmin)/n;
+  }
+  if(lkmax-lkfind<dlk/2.) rc |= 2;  // f(kmax) >= high
+  else kmax = pow(10.,lkmax);
+  //cout<<"rc="<<rc<<" lkmax="<<lkmax<<"  pk="<<(*this)(pow(10.,lkmax))<<endl;
+
+  return rc;
+}
Index: trunk/Cosmo/SimLSS/pkspectrum.h
===================================================================
--- trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3115)
@@ -0,0 +1,128 @@
+#ifndef PKSPECTRUM_SEEN
+#define PKSPECTRUM_SEEN
+
+#include "machdefs.h"
+#include "genericfunc.h"
+
+namespace SOPHYA {
+
+//-----------------------------------------------------------------------------------
+class InitialSpectrum : public GenericFunc {
+public:
+  InitialSpectrum(double n,double a=1.);
+  InitialSpectrum(InitialSpectrum& pkinf);
+  virtual ~InitialSpectrum(void);
+  virtual double operator() (double k) {return A_ * pow(k,n_);}
+  void SetNorm(double a);
+  void SetSlope(double n);
+protected:
+  double n_, A_;
+};
+
+//-----------------------------------------------------------------------------------
+class TransfertEisenstein : public GenericFunc {
+public:
+  TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon=false);
+  TransfertEisenstein(TransfertEisenstein& tf);
+  virtual ~TransfertEisenstein(void);
+  virtual double operator() (double k);
+  double KPeak(void);
+  void SetNoOscEnv(unsigned short nooscenv=0);
+protected:
+  double O0_,Oc_,Ob_,h_,tcmb_,th2p7_;
+  double zeq_,keq_,zd_,Req_,Rd_,s_,ksilk_,alphac_,betac_,bnode_,alphab_,betab_;
+  double alphag_;
+  double kpeak_;
+
+  bool nobaryon_;
+  unsigned short nooscenv_;
+  double T0tild(double k,double alphac,double betac);
+  void Init_With_Baryons(void);
+  void Init_Without_Baryon(void);
+};
+
+
+//-----------------------------------------------------------------------------------
+class GrowthFactor : public GenericFunc {
+public:
+  GrowthFactor(double OmegaMatter0,double OmegaLambda0);
+  GrowthFactor(GrowthFactor& d1);
+  virtual ~GrowthFactor(void);
+  virtual double operator() (double z);
+protected:
+  double O0_,Ol_,Ok_;
+  double norm_;
+};
+
+
+//-----------------------------------------------------------------------------------
+class PkSpectrum0 : public GenericFunc {
+public:
+  PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf);
+  PkSpectrum0(PkSpectrum0& pk0);
+  virtual ~PkSpectrum0(void);
+  virtual double operator() (double z);
+  InitialSpectrum& GetPkIni(void) {return pkinf_;}
+  TransfertEisenstein& GetTransfert(void) {return tf_;}
+protected:
+  InitialSpectrum& pkinf_;
+  TransfertEisenstein& tf_;
+};
+
+//-----------------------------------------------------------------------------------
+class PkSpectrumZ : public GenericFunc {
+public:
+  PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref=0.);
+  PkSpectrumZ(PkSpectrumZ& pkz);
+  virtual ~PkSpectrumZ(void);
+  virtual double operator() (double k);
+  virtual double operator() (double k,double z);
+  inline void   SetZ(double z) {zref_ = z;}
+  inline double GetZ(void) {return zref_;}
+  void SetTypSpec(unsigned short typspec=0);
+  inline void SetScale(double scale=1.) {scale_=scale; zold_=-1.;}
+  inline double GetScale(void) {return scale_;}
+  PkSpectrum0& GetPk0(void) {return pk0_;}
+  GrowthFactor& GetGrowthFactor(void) {return d1_;}
+protected:
+  PkSpectrum0& pk0_;
+  GrowthFactor& d1_;
+  double zref_, scale_;
+  unsigned short typspec_;
+  mutable double zold_, d1old_;
+};
+
+//-----------------------------------------------------------------------------------
+class VarianceSpectrum : public GenericFunc {
+public:
+  VarianceSpectrum(GenericFunc& pk,unsigned short typfilter);
+  VarianceSpectrum(VarianceSpectrum& pkinf);
+  virtual ~VarianceSpectrum(void);
+
+  void SetFilter(unsigned short typfilter=0);
+  void SetInteg(double dperc=0.1,double dlogkinc=-1.,double dlogkmax=-1.,unsigned short glorder=4);
+
+  double Variance(double R,double kmin,double kmax);
+
+  // ATTENTION: La fonction a integrer est : f(k)dk = k^3*Pk(k)/(2Pi^2) *filter2(k*R) *dk/k
+  virtual double operator() (double k) {return k*k*pk_(k)*Filter2(k*R_)/(2.*M_PI*M_PI);}
+  double Filter2(double x);
+
+  // Aide a l'integration
+  double FindMaximum(double R,double kmin,double kmax,double eps=1.e-3);
+  int FindLimits(double R,double high,double &kmin,double &kmax,double eps=1.e-3);
+
+protected:
+
+  GenericFunc& pk_;
+  unsigned short typfilter_;
+  double R_;
+
+  double dperc_,dlogkinc_,dlogkmax_;
+  unsigned short glorder_;
+
+};
+
+} // Fin du namespace
+
+#endif
Index: trunk/Cosmo/SimLSS/planckspectra.cc
===================================================================
--- trunk/Cosmo/SimLSS/planckspectra.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/planckspectra.cc	(revision 3115)
@@ -0,0 +1,337 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "pexceptions.h"
+
+#include "constcosmo.h"
+#include "planckspectra.h"
+
+/*
+----------------------------------------------------------------
+- UNITES:
+     Temperature en Kelvin  "K"
+     Frequence en GigaHertz  "Hz"
+     Longueur d'onde en Metre  "m"
+----------------------------------------------------------------
+- Frequence / Longeur d'onde
+ l = c/f ; f = c/l
+ dl = c/f^2 df = l^2/c df
+ df = c/l^2 dl = f^2/c dl
+ P[f] df = P[l] dl = P[f] c/l^2 dl = P[l] c/f^2 df
+----------------------------------------------------------------
+- Lois de Planck Wien et Rayleigh
+     Planck   H*Nu/(K*T) = H*C/(Lambda*K*T) quelconque
+     Wien     H*Nu/(K*T) = H*C/(Lambda*K*T) >> 1
+     Rayleigh H*Nu/(K*T) = H*C/(Lambda*K*T) << 1
+----------------------------------------------------------------
+- Les termes en facteur pour les lois de Planck Wien et Rayleigh
+        (la partie contenant l'exponentielle)
+ ainsi que pour les derivees des ces termes par rapport a la temperature
+ Planck:   1/[exp(H*f/(K*t))-1] = 1/[exp(HC/(l*K*t))-1]
+     der=  H*f/(K*t) /t * exp(HC/(l*K*t))/[exp(HC/(l*K*t))-1]^2
+ Wien:     1/exp(H*f/(K*t))
+     der=  H*f/(K*t) /t /exp(H*f/(K*t))
+ Rayleigh: 1/(H*f/(K*t))
+     der=  K/(H*f)
+ Pour x->0: compute exp(x)=1+x+x^2/2 for 1/(exp(x)-1)
+            soit 1/(exp(x)-1)  -> 1/(x+x^2/2)
+  car exp(x)-1 vaut zero plus vite que son DL a cause de la soustraction
+----------------------------------------------------------------
+ - RADIANCE (Brillance) en fonction de nu ou lambda
+ nombre de watts (ou photons par seconde)
+ radies par unite de surface par steradian et par hertz (ou metre)
+    2*f^2/c^2   * 1/(exp()-1)    ph/s/m^2/Sr/Hz
+    2*h*f^3/c^2 * 1/(exp()-1)    W/m^2/Sr/Hz
+    2*c/l^4     * 1/(exp()-1)    ph/s/m^2/Sr/m
+    2*h*c^2/l^5 * 1/(exp()-1)    W/m^2/Sr/m
+----------------------------------------------------------------
+ - EMITTANCE en fonction de nu ou lambda
+ nombre de watts (ou photons par seconde)
+ radies dans un hemisphere (2*PI sr)
+ par unite de surface et par hertz (ou metre)
+ ..quand on integre sur 2Pi sr il faut
+ ..multiplier par cos(theta)  -> facteur Pi uniquement
+    EMITTANCE = Pi * RADIANCE
+    2*Pi*f^2/c^2   * 1/(exp()-1)    ph/s/m^2/Hz
+    2*Pi*h*f^3/c^2 * 1/(exp()-1)    W/m^2/Hz
+    2*Pi*c/l^4     * 1/(exp()-1)    ph/s/m^2/m
+    2*Pi*h*c^2/l^5 * 1/(exp()-1)    W/m^2/m
+----------------------------------------------------------------
+ - DENSITE en fonction de nu ou lambda
+ nombre de joules (ou photons) par m^3  et par hertz (ou metre)
+   DENSITE = 4*Pi/c * RADIANCE
+   8*Pi*f^2/c^3    * 1/(exp()-1)    ph/m^3/Hz
+   8*Pi*h*f^3/c^3  * 1/(exp()-1)    J/m^3/Hz
+   8*Pi/l^4        * 1/(exp()-1)    ph/m^3/m
+   8*Pi*h*c/l^5    * 1/(exp()-1)    J/m^3/m
+----------------------------------------------------------------
+ - DENSITE d'energie totale  (loi de Stefan-Boltzman)
+   4*sigma*T^4/C   J/m^3
+   Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
+----------------------------------------------------------------
+*/
+
+/*
+POUR TESTER LES LIMITES DE L'EXPONENTIELLE:
+1-/ copier/coller les lignes suivantes dans toto.icc
+double x=10.;
+for(int i=0;i<35;i++) {
+  double ex = exp(x);
+  double exm1a = x*(1.+x/2.*(1.+x/3.)),;
+  printf("x=%g exp(x)=%e exp(x)-1=%e dl=%e (d=%e)\n",x,ex,ex-1.,exm1a,ex-1.-exm1a);
+  x /= 10.;
+}
+2-/ executer sous spiapp: > c++execfrf toto.icc
+*/
+
+
+#define MAXARGEXP 100.    /* maxi = 709.7827128933840 */
+#define MINARGEXP 1.e-10  /* Replace exp(x) by DL for 1/(exp(x)-1) when x->0 */
+
+//************************************************************************
+// vitesse de la lumiere en m/s
+double CinM_ = SpeedOfLight_Cst*1e+3;
+// h/k pour une frequence en Hz
+double HsK_ = h_Planck_Cst/k_Boltzman_Cst;
+// hc/k pour une longeur d'onde en metre
+double HCsK_ = h_Planck_Cst*CinM_/k_Boltzman_Cst;
+
+//************************************************************************
+PlanckSpectra::PlanckSpectra(double T)
+  : T_(T) , approx_(0) , deriv_(0) , typvar_(0) , typspec_(0) , unitout_(0)
+{
+}
+
+PlanckSpectra::PlanckSpectra(PlanckSpectra& s)
+  : T_(s.T_) , approx_(s.approx_) , deriv_(s.deriv_) , typvar_(s.typvar_) , typspec_(s.typspec_)
+  , unitout_(s.unitout_)
+{
+}
+
+PlanckSpectra::~PlanckSpectra(void)
+{
+}
+
+//************************************************************************
+void PlanckSpectra::SetApprox(unsigned short approx)
+// approx = 0 : Planck spectra
+//        = 1 : Rayleigh spectra
+//        = 2 : Wien spectra
+{
+  if(approx>2) {
+    cout<<"PlanckSpectra::SetApprox : unknown approximation"<<endl;
+    throw ParmError("PlanckSpectra::SetApprox : unknown approximation");
+  }
+  approx_ = approx;
+}
+
+void PlanckSpectra::SetDeriv(unsigned short deriv)
+// deriv = 0 : spectra
+//       = 1 : derivative of spectra relative to temperature
+{
+  if(deriv>1) {
+    cout<<"PlanckSpectra::SetDeriv : parameter out of range"<<endl;
+    throw ParmError("PlanckSpectra::SetDeriv : parameter out of range");
+  }
+  deriv_ = deriv;
+}
+
+void PlanckSpectra::SetVar(unsigned short typvar)
+// typvar = 0 : variable is frequency (Hz)
+//        = 1 : variable is wavelength (m)
+{
+  if(typvar>1) {
+    cout<<"PlanckSpectra::SetVar : parameter out of range"<<endl;
+    throw ParmError("PlanckSpectra::SetVar : parameter out of range");
+  }
+  typvar_ = typvar;
+}
+
+void PlanckSpectra::SetUnitOut(unsigned short unitout)
+// unitout = 0 : in W/...
+//         = 1 : in photon/s/...
+{
+  if(unitout>1) {
+    cout<<"PlanckSpectra::SetUnitOut : parameter out of range"<<endl;
+    throw ParmError("PlanckSpectra::SetUnitOut : parameter out of range");
+  }
+  unitout_ = unitout;
+}
+
+void PlanckSpectra::SetTypSpectra(unsigned short typspec)
+// typspec = 0 : radiance   W/m^2/Sr/Hz  or ...
+//         = 1 : emittance  W/m^2/Hz     or ... 
+//         = 2 : density    J/m^3/Hz     or ...
+//         = 3 : just the function exp() or equivalent for Rayleigh
+{
+  if(typspec>3) {
+    cout<<"PlanckSpectra::SetTypSpectra : parameter out of range"<<endl;
+    throw ParmError("PlanckSpectra::SetTypSpectra : parameter out of range");
+  }
+  typspec_ = typspec;
+}
+
+
+//************************************************************************
+double PlanckSpectra::F2L(double nu_or_lambda)
+// Convert frequency (GHz) to/from Wavelength (m)
+// Lambda(m) * Nu(GHz) = C(km/s) * 10^-6
+{
+  return CinM_ / nu_or_lambda;
+}
+
+//************************************************************************
+double PlanckSpectra::PlanckExp(double fl)
+{
+  double x = (typvar_) ? HCsK_/(fl*T_) : HsK_*fl/T_;
+
+  if(approx_ == 1) return 1./x;  // Rayleigh
+
+  if(x>MAXARGEXP) return 0.;
+
+  if(approx_ == 2) return exp(-x);  // Wien
+
+  if(x<MINARGEXP) return 1./(x*(1.+x/2.));
+  return 1./(exp(x)-1.);
+}
+
+double PlanckSpectra::DPlanckExp_DT(double fl)
+{
+  if(approx_ == 1) return (typvar_) ? fl/HCsK_: 1. / (HsK_*fl);  // Rayleigh
+
+  double x = (typvar_) ? HCsK_/(fl*T_) : HsK_*fl/T_;
+
+  if(x>MAXARGEXP) return 0.;
+
+  if(approx_ == 2) return x/T_ * exp(-x);  // Wien
+
+  double v;
+  if(x<MINARGEXP) { // 1/(exp(x)-1)  -derive-> -exp(x)/(exp(x)-1)^2
+    v = x*(1.+x/2.);
+    return x/T_ * (1.+v)/(v*v);
+  }
+  v = exp(x);
+  return x/T_ * v/((v-1.)*(v-1.));  
+}
+
+//------------------------------------------------------------
+double PlanckSpectra::operator() (double x)
+{
+  // --- type de spectre et d'approximation selon choix variable
+  double fex = (deriv_) ? DPlanckExp_DT(x): PlanckExp(x);
+
+  // --- le coefficient multiplicateur devant l'expo
+  if(typspec_ == 3) return fex;  // expo seule
+
+  // radiance en ph/...       longeur d'onde         frequency
+  double coeff = (typvar_) ? 2.*CinM_/(x*x*x*x): 2.*x*x/(CinM_*CinM_);
+
+  if(typspec_ == 1)      coeff *= M_PI;  // emittance en ph/...
+  else if(typspec_ == 2) coeff *= 4.*M_PI/CinM_;  // densite en ph/...
+
+  // unite en W/...                       longeur d'onde         frequency
+  if(unitout_ == 0) coeff *= (typvar_) ? h_Planck_Cst*CinM_/x: h_Planck_Cst*x;
+
+  return coeff * fex;
+}
+
+//------------------------------------------------------------
+double PlanckSpectra::PlanckEnergie(void)
+//  densite d'energie: energie par m^3  en "J / m^3"
+// cste de Stefan-Boltzmann en W/m^2/K^4
+{
+  return 4.*StefBoltz_Cst*T_*T_*T_*T_/(SpeedOfLight_Cst*1000.);
+}
+
+//------------------------------------------------------------
+double PlanckSpectra::PlanckPhoton(void)
+//  densite de photons: nombre de photons par m^3  en "ph / m^3"
+{
+  double x = k_Boltzman_Cst*T_/(h_Planck_Cst*SpeedOfLight_Cst*1000.);
+  return 16.*M_PI*Zeta_3_Cst* x*x*x;
+}
+
+double PlanckSpectra::WienLaw(void)
+// Retourne la frequence (ou longeur d'onde) du maximum
+// du spectre de radiance (ou emittance ou densite).
+// La valeur du maximum depend de l'unite du spectre (W/... ou Ph/s/...)
+/*
+--- Les fonctions:
+En frequence: n = 3  en W/...  ,  n = 2  en Ph/...
+En longueur d'onde: n = 5  en W/...  ,  n = 4  en Ph/...
+f[x_] := x^n*Exp[-x]; Solve[f'[x]==0,x]   // Wien
+f[x_] := x^n/(Exp[x]-1); NSolve[f'[x]==0,x,17]  // Planck
+La derivee par rapport a T des fonctions en frequence:
+df/dT[x_] := x^(n+1)*Exp[-x]; Solve[f'[x]==0,x]  // Wien
+df/dT[x_] := x^(n+1)*Exp[x]/(Exp[x]-1)^2; NSolve[f'[x]==0,x,17]  // Planck
+*/
+{
+  double x=-1.;
+  if(typvar_ == 0) { // Frequence
+    if(deriv_ == 0) {  // Le spectre
+      if(unitout_ == 0) //  W/...
+        {if(approx_==0) x=2.82143938074; else if(approx_==2) x=3.;}
+      else if(unitout_ == 1) //  Ph/s/...
+        {if(approx_==0) x=1.59362426; else if(approx_==2) x=2.;}
+    } else if(deriv_ == 1) {  // La derivee du spectre
+      if(unitout_ == 0) //  W/... 
+        {if(approx_==0) x=3.830016095; else if(approx_==2) x=4.;}
+      else if(unitout_ == 1) //  Ph/s/...
+        {if(approx_==0) x=2.575678888; else if(approx_==2) x=3.;}
+    }
+  } else { // Longeur d'onde
+    if(deriv_ == 0) {  // Le spectre
+      if(unitout_ == 0) //  W/...
+        {if(approx_==0) x=4.96511420; else if(approx_==2) x=5.;}
+      else if(unitout_ == 1) //  Ph/s/...
+        {if(approx_==0) x=3.9206904099; else if(approx_==2) x=4.;}
+    } else if(deriv_ == 1) {  // La derivee du spectre
+      if(unitout_ == 0) //  W/...
+        {if(approx_==0) x=5.969409190; else if(approx_==2) x=6.;}
+     else if(unitout_ == 1) //  Ph/s/...
+        {if(approx_==0) x=4.928119345; else if(approx_==2) x=5.;}
+    }
+  }
+
+  if(x<0.) return x;  // Rayleigh n'a pas de maximum
+  if(typvar_ == 0) return x / HsK_ *T_; // Frequence
+  return HCsK_ / x /T_; // Longeur d'onde
+}
+
+double PlanckSpectra::FindMaximum(double eps)
+// Retourne le maximum des spectres (a la precicion eps)
+{
+  if(approx_ == 1) {  // Pas de maximum pour les spectres de Raleigh
+    cout<<"PlanckSpectra::FindMaximum : no maximum for Raleigh spectra"<<endl;
+    throw ParmError("PlanckSpectra::FindMaximum : no maximum for Raleigh spectra");
+  }
+  if(eps<=0.) eps = 0.001;
+
+  double xmax = WienLaw(), wmax=(*this)(xmax);
+  double xlim[2]={xmax,xmax};
+
+  // Encadrement du maximum par xlim[0] et xlim[1]
+  while( (*this)(xlim[1])>=0.01*wmax ) xlim[1] *= 2.;
+  while( (*this)(xlim[0])>=0.01*wmax ) xlim[0] /= 2.;
+  //cout<<"FindMaximum: wmax="<<wmax<<" pour x=["<<xlim[0]<<","<<xlim[1]<<"]"<<endl;
+
+  while( (xlim[1]-xlim[0])/xmax > eps ) {
+   wmax=-1.;
+   double dx = (xlim[1]-xlim[0])/10.;
+   for(double x=xlim[0]; x<=xlim[1]; x+=dx) {
+     double f=(*this)(x);
+     if(f>wmax) {wmax=f; xmax=x;}
+   }
+   //cout<<"...wmax="<<wmax<<" pour xmax="<<xmax<<" dx="<<dx<<endl;
+   xlim[0] = xmax-dx*1.1; xlim[1] = xmax+dx*1.1;
+ }
+
+
+  //cout<<"wmax="<<wmax<<" pour x="<<xmax<<"  dx="<<xlim[1]-xlim[0]<<endl;
+  return xmax;
+}
Index: trunk/Cosmo/SimLSS/planckspectra.h
===================================================================
--- trunk/Cosmo/SimLSS/planckspectra.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/planckspectra.h	(revision 3115)
@@ -0,0 +1,45 @@
+#ifndef PLANCKSPECTRA_SEEN
+#define PLANCKSPECTRA_SEEN
+
+#include "machdefs.h"
+#include "genericfunc.h"
+
+namespace SOPHYA {
+
+//-----------------------------------------------------------------------------------
+class PlanckSpectra : public GenericFunc {
+public:
+  PlanckSpectra(double T);
+  PlanckSpectra(PlanckSpectra& s);
+  virtual ~PlanckSpectra(void);
+
+  inline void SetTemperature(double T) {T_ = T;}
+  inline double GetTemperature(void) {return T_;}
+
+  void SetApprox(unsigned short approx=0);
+  void SetDeriv(unsigned short deriv=0);
+  void SetVar(unsigned short typvar=0);
+  void SetTypSpectra(unsigned short typspec=0);
+  void SetUnitOut(unsigned short unitout=0);
+
+
+  virtual double operator() (double x);
+
+  double F2L(double nu_or_lambda);
+  double PlanckEnergie(void);
+  double PlanckPhoton(void);
+  double WienLaw(void);
+  double FindMaximum(double eps=0.001);
+
+protected:
+  double T_;
+  unsigned short approx_,deriv_,typvar_,typspec_,unitout_;
+
+  double PlanckExp(double fl);
+  double DPlanckExp_DT(double fl);
+
+};
+
+} // Fin du namespace
+
+#endif
Index: trunk/Cosmo/SimLSS/schechter.cc
===================================================================
--- trunk/Cosmo/SimLSS/schechter.cc	(revision 3115)
+++ trunk/Cosmo/SimLSS/schechter.cc	(revision 3115)
@@ -0,0 +1,74 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "pexceptions.h"
+
+#include "constcosmo.h"
+#include "schechter.h"
+
+///////////////////////////////////////////////////////////
+//***************** Schechter Functions *****************//
+///////////////////////////////////////////////////////////
+
+// HI mass function:
+//    see M.Zwaan astroph-0502257 (MNRAS Letters, Volume 359, Issue 1, pp. L30-L34.)
+//    see M.Zwaan astroph-9707109 (ApJ, Volume 509, Issue 2, pp. 687-702.)
+
+Schechter::Schechter(double nstar,double mstar,double alpha)
+  : nstar_(nstar) , mstar_(mstar) , alpha_(alpha) , outvalue_(0)
+{
+}
+
+Schechter::Schechter(Schechter& f)
+  : nstar_(f.nstar_) , mstar_(f.mstar_) , alpha_(f.alpha_) , outvalue_(f.outvalue_)
+{
+}
+
+Schechter::~Schechter(void)
+{
+}
+
+void Schechter::SetOutValue(unsigned short outvalue)
+// outvalue = 0 : give dn/dm
+//          = 1 : give m*dn/dm
+{
+  if(outvalue>1) {
+    cout<<"Schechter::SetOutValue: Error bad outvalue: "<<outvalue<<endl;
+    throw ParmError("Schechter::SetOutValue: Error bad outvalue");
+  }
+  outvalue_ = outvalue;
+}
+
+double Schechter::operator() (double m)
+// Return : "dn/dm = f(m)" or "m*dn/dm = f(m)"
+{
+  double x = m/mstar_;
+  double dndm = nstar_/mstar_ * pow(x,alpha_) * exp(-x);
+  if(outvalue_) dndm *= m;  // on veut m*dn/dm
+  return dndm;
+}
+
+///////////////////////////////////////////////////////////
+//******************* Le Flux a 21cm ********************//
+///////////////////////////////////////////////////////////
+
+double FluxHI(double m,double d)
+// Input:
+//    m : masse de HI en "Msol"
+//    d : distance en "Mpc"  (si cosmo d=DLum(z))
+// Return:
+//    le flux total emis a 21 cm en W/m^2
+// Ref: Binney & Merrifield, Galactic Astronomy p474 (ed:1998)
+//      J.Peterson & K.Bandura, astroph-0606104  (eq 7)
+//      F.Abdalla & Rawlings, astroph-0411342 (eq 7 may pb de d'unites?)
+{
+  return  m / (2.35e+5 * d*d )
+          * Jansky2Watt_cst
+          * (Fr_HyperFin_Par*1e+9)/SpeedOfLight_Cst;
+}
+
Index: trunk/Cosmo/SimLSS/schechter.h
===================================================================
--- trunk/Cosmo/SimLSS/schechter.h	(revision 3115)
+++ trunk/Cosmo/SimLSS/schechter.h	(revision 3115)
@@ -0,0 +1,28 @@
+#ifndef SCHECHTER_SEEN
+#define SCHECHTER_SEEN
+
+#include "machdefs.h"
+#include "genericfunc.h"
+
+
+namespace SOPHYA {
+
+//-----------------------------------------------------------------------------------
+class Schechter : public GenericFunc {
+public:
+  Schechter(double nstar,double mstar,double alpha);
+  Schechter(Schechter& f);
+  virtual ~Schechter(void);
+  void SetOutValue(unsigned short outvalue=0);
+  virtual double operator() (double m);
+protected:
+  double nstar_,mstar_,alpha_;
+  unsigned short outvalue_;
+};
+
+} // Fin du namespace
+
+//-----------------------------------------------------------------------------------
+double FluxHI(double m,double d);
+
+#endif
