Index: trunk/Cosmo/SimLSS/Makefile
===================================================================
--- trunk/Cosmo/SimLSS/Makefile	(revision 3119)
+++ trunk/Cosmo/SimLSS/Makefile	(revision 3120)
@@ -19,5 +19,6 @@
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
        $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
-       $(EXE)cmvtluc $(EXE)cmvtpoisson $(EXE)cmvconchprof
+       $(EXE)cmvtpoisson $(EXE)cmvconchprof
+#$(EXE)cmvtluc
  
 PROGSOBJ = \
@@ -25,5 +26,5 @@
           $(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
+          $(OBJ)cmvtpoisson.o $(OBJ)cmvconchprof.o $(OBJ)cmvtluc.o
  
 LIBROBJ = \
Index: trunk/Cosmo/SimLSS/cmvconchprof.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvconchprof.cc	(revision 3119)
+++ trunk/Cosmo/SimLSS/cmvconchprof.cc	(revision 3120)
@@ -88,4 +88,8 @@
 
 disp hpconc
+
 n/plot hpconc.val%log10(x) x>0 ! "connectpoints"
+n/plot hpconc.err%log10(x) x>0 ! "connectpoints same red"
+
+n/plot hpconc.err/val%log10(x) x>0&&val>0 ! "connectpoints"
 */
Index: trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3119)
+++ trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3120)
@@ -12,4 +12,5 @@
 #include "srandgen.h"
 #include "perandom.h"
+#include "fabtwriter.h"
 
 #include "arrctcast.h"
@@ -31,5 +32,5 @@
      <<" -z nz,dz : taille en z (npix,Mpc)"<<endl
      <<" -Z zref : redshift median"<<endl
-     <<" -s scale : sigma du bruit en Msol"<<endl
+     <<" -s snoise : sigma du bruit en Msol"<<endl
      <<endl;
 }
@@ -163,5 +164,5 @@
  pkz.SetZ(zref);
 
- Histo hpkz (lkmin,lkmax,npt); hpkz.ReCenterBin();
+ Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
  FuncToHisto(pkz,hpkz,true);
  {
@@ -330,4 +331,25 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
+ if(1) {
+   cout<<"\n---Writing Cube to Fits file"<<endl;
+   FitsImg3DWriter fwrt("!cmvobserv3dr.fits",FLOAT_IMG,5);
+   fwrt.WriteKey("ZREF",zref," redshift");
+   vector<long> n = fluct3d.GetNpix();
+   fwrt.WriteKey("NX",n[0]," axe transverse 1");
+   fwrt.WriteKey("NY",n[1]," axe transverse 2");
+   fwrt.WriteKey("NZ",n[2]," axe longitudinal (redshift)");
+   vector<r_8> d;
+   d = fluct3d.GetDinc();
+   fwrt.WriteKey("DX",d[0]," Mpc");
+   fwrt.WriteKey("DY",d[1]," Mpc");
+   fwrt.WriteKey("DZ",d[2]," Mpc");
+   d = fluct3d.GetKinc();
+   fwrt.WriteKey("DKX",d[0]," Mpc^-1");
+   fwrt.WriteKey("DKY",d[1]," Mpc^-1");
+   fwrt.WriteKey("DKZ",d[2]," Mpc^-1");
+   fwrt.WriteKey("SNOISE",snoise," Msol");
+   fwrt.Write(rgen);
+ }
+
  cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
  fluct3d.AddNoise2Real(snoise);
@@ -364,8 +386,8 @@
 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"
+n/plot hpkz.val%x ! ! "nsta connectpoints"
+n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
+n/plot hpkgenf.val%log10(x) x>0 ! "nsta same orange connectpoints"
+n/plot hpkrec.val%log10(x) x>0 ! "nsta same blue connectpoints"
 
 set k pow(10.,x)
@@ -379,3 +401,7 @@
 rgensl 0
 
+disp hmdndm
+disp tirhmdndm
+addline 0 1 20 1 "red"
+
  */
Index: trunk/Cosmo/SimLSS/cmvtstpk.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3119)
+++ trunk/Cosmo/SimLSS/cmvtstpk.cc	(revision 3120)
@@ -19,5 +19,6 @@
 void usage(void);
 void usage(void) {
- cout<<"cmvtstpk -k npt,lkmin,lkmax -s scale zval"<<endl
+ cout<<"cmvtstpk [options] z_redshift"<<endl
+     <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
      <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
      <<" -s scale : on multiplie pkz par scale"<<endl;
@@ -26,21 +27,34 @@
 int main(int narg,char *arg[])
 {
- double ob0 = 0.0444356;
+ double Ob0 = 0.0444356;
  // -- WMAP
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double h100=0.71, Om0=0.267804, Ol0=0.73,w0=-1.;
  // -- ouvert matter only
- //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ //double h100=0.71, Om0=0.3, Ol0=0.,w0=-1.;
  // -- plat matter only
- //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;
+ //double h100=0.71, Om0=1., Ol0=0.,w0=-1.;
+
 
  double ns = 1., as = 1.;
 
- int npt = 1000;
- double lkmin = -4., lkmax=2.;
+ int npt = 10000;
+ double lkmin = -3., lkmax=2.;
  double scale = 1.;
 
  char c;
-  while((c = getopt(narg,arg,"hk:s:")) != -1) {
+  while((c = getopt(narg,arg,"hk:s:H:M:B:L:")) != -1) {
   switch (c) {
+  case 'H' :
+    sscanf(optarg,"%lf",&h100);
+    break;
+  case 'M' :
+    sscanf(optarg,"%lf",&Om0);
+    break;
+  case 'B' :
+    sscanf(optarg,"%lf",&Ob0);
+    break;
+  case 'L' :
+    sscanf(optarg,"%lf,%lf",&Ol0,&w0);
+    break;
   case 'k' :
     sscanf(optarg,"%d,%lf,%lf",&npt,&lkmin,&lkmax);
@@ -62,4 +76,5 @@
  if(optind<narg) zval = atof(arg[optind]);
 
+ cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Ol0="<<Ol0<<" w0="<<w0<<endl;
  cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
  cout<<"scale="<<scale<<endl;
@@ -69,11 +84,11 @@
  InitialSpectrum pkini(ns,as);
 
- TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+ 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);
+ TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
+
+ GrowthFactor d1(Om0,Ol0);
  cout<<"GrowthFactor: "<<d1(zval)<<endl;
 
@@ -89,5 +104,5 @@
 
  //--------------------------
- Histo hd1(0.,1000.,10000); hd1.ReCenterBin();
+ Histo hd1(0.,20.,10000); hd1.ReCenterBin();
  FuncToHisto(d1,hd1,false);
 
@@ -99,5 +114,5 @@
  FunRan talea(hpkz,true);
  Histo halea(hpkz); halea.Zero();
- int nalea = 1000000;
+ int nalea = 100000;
  for(int i=0;i<nalea;i++) halea.Add(talea.Random());
  halea *= hpkz.Sum()/halea.Sum();
@@ -169,51 +184,55 @@
 #### growth-factor
 zone
-n/plot hd1.log10(val)%log10(x) val>0.&&x>0. ! "nsta connectpoints"
+n/plot hd1.val%x ! ! "nsta connectpoints"
 
 #### Spectre initial
 zone
-n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta crossmarker3"
+n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta connectpoints"
+
+#### Gestion des abscisses
+set xk log10(k)
+set xk k
 
 #### 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"
+n/plot nt.tf%$xk ! ! "nsta connectpoints"
+n/plot nt.tfnosc2%$xk ! ! "nsta connectpoints same red"
+n/plot nt.tfnosc1%$xk ! ! "nsta connectpoints same blue"
+n/plot nt.tfnob%$xk ! ! "nsta connectpoints same green"
+
+n/plot nt.tf/tfnosc2%$xk ! ! "nsta connectpoints red"
+n/plot nt.tf/tfnosc1%$xk ! ! "nsta connectpoints same blue"
+n/plot nt.tf/tfnob%$xk ! ! "nsta connectpoints same green"
 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"
+n/plot nt.pk0%$xk ! ! "nsta connectpoints"
+n/plot nt.pk0nosc2%$xk ! ! "nsta connectpoints same red"
+n/plot nt.pk0nosc1%$xk ! ! "nsta connectpoints same blue"
+n/plot nt.pk0nob%$xk ! ! "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"
+n/plot nt.pk0/pkini-tf*tf%$xk pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%$xk pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%$xk pkini>0 ! "nsta crossmarker3"
+n/plot nt.pk0nob/pkini-tfnob*tfnob%$xk 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"
+n/plot nt.pk%$xk ! ! "nsta connectpoints"
+n/plot nt.pknosc2%$xk ! ! "nsta connectpoints same red"
+n/plot nt.pknosc1%$xk ! ! "nsta connectpoints same blue"
+n/plot nt.pknob%$xk ! ! "nsta connectpoints same green"
+
+n/plot nt.pk/pknosc2%$xk ! ! "nsta connectpoints red"
+n/plot nt.pk/pknosc1%$xk ! ! "nsta connectpoints same blue"
+n/plot nt.pk/pknob%$xk ! ! "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"
+n/plot nt.$D2%$xk  !  ! "nsta crossmarker3 connectpoints"
 
 #### Test des transferts dans Histo et TVector
Index: trunk/Cosmo/SimLSS/cmvtstsch.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3119)
+++ trunk/Cosmo/SimLSS/cmvtstsch.cc	(revision 3120)
@@ -94,5 +94,5 @@
  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;
+   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
  }
  Histo hdndm(lnx1,lnx2,npt); hdndm.ReCenterBin();
@@ -105,5 +105,5 @@
  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;
+   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
  }
  Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
Index: trunk/Cosmo/SimLSS/constcosmo.h
===================================================================
--- trunk/Cosmo/SimLSS/constcosmo.h	(revision 3119)
+++ trunk/Cosmo/SimLSS/constcosmo.h	(revision 3120)
@@ -10,5 +10,5 @@
 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"
+static const double SolarMass_Cst = 1.98844e+33; // Solar 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"
Index: trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3119)
+++ trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3120)
@@ -344,9 +344,9 @@
 void GeneFluct3D::FilterByPixel(void)
 // Filtrage par la fonction fenetre du pixel (parallelepipede)
-// TF = 1/(dx*dy*dz)*Int[{0,dx},{0,dy},{0,dz}]
+// TF = 1/(dx*dy*dz)*Int[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
 //                   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]
+//    = 2/(k_x*dx) * sin(k_x*dx/2)  * (idem y) * (idem z)
+// Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
+//                          avec y = k_x*dx/2
 {
  int lp=2;
@@ -355,14 +355,14 @@
  for(int i=0;i<Nx_;i++) {
    int ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_ *Dx_;
-   double pkx2 = pixelfilter(kx); 
+   double kx = ii*Dkx_ *Dx_/2;
+   double pkx = 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); 
+     double ky = jj*Dky_ *Dy_/2;
+     double pky =  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;
+       double kz = l*Dkz_ *Dz_/2;
+       double pkz =  pixelfilter(kz);
+       T_(l,j,i) *= pkx*pky*pkz;
      }
    }
Index: trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3119)
+++ trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3120)
@@ -35,4 +35,8 @@
   vector<r_8> GetKnyq(void)
     {vector<r_8> kn; kn.push_back(Knyqx_); kn.push_back(Knyqy_); kn.push_back(Knyqz_); return kn;}
+  vector<r_8> GetDinc(void)
+    {vector<r_8> d; d.push_back(Dx_); d.push_back(Dy_); d.push_back(Dz_); return d;}
+  vector<long> GetNpix(void)
+    {vector<long> d; d.push_back(Nx_); d.push_back(Ny_); d.push_back(Nz_); return d;}
 
   void ComputeFourier0(void);
@@ -60,6 +64,6 @@
   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);}
+  inline double pixelfilter(double x)
+    {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
 
 
Index: trunk/Cosmo/SimLSS/geneutils.h
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.h	(revision 3119)
+++ trunk/Cosmo/SimLSS/geneutils.h	(revision 3120)
@@ -16,4 +16,7 @@
   InterpFunc(double xmin,double xmax,vector<double>& y);
   virtual ~InterpFunc(void) { }
+
+  double XMin(void) {return _xmin;}
+  double XMax(void) {return _xmax;}
 
   //! Retourne l'element le plus proche de f donnant y=f(x)
