Index: /trunk/Cosmo/SimLSS/Makefile
===================================================================
--- /trunk/Cosmo/SimLSS/Makefile	(revision 3767)
+++ /trunk/Cosmo/SimLSS/Makefile	(revision 3768)
@@ -11,6 +11,6 @@
 
 MYEXTINC = ${EXTLIBDIR}/${MACH}/include
-MYLIB  = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lcmvgenfluc -lfftw3_threads -lfftw3 -lm
-MYLIB4 = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lcmvgenfluc4 -lfftw3f_threads -lfftw3f -lm
+MYLIB  = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvgenfluc -lcmvsimbao -lfftw3_threads -lfftw3 -lm
+MYLIB4 = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvgenfluc4 -lcmvsimbao -lfftw3f_threads -lfftw3f -lm
 
 #--------------------------------------------------------------------------
@@ -18,4 +18,5 @@
 PROGS = \
        $(EXE)cmvobserv3d $(EXE)cmvobserv3df \
+       $(EXE)cmvginit3d $(EXE)cmvginit3df \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
@@ -25,4 +26,5 @@
 PROGSOBJ = \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvobserv3df.o \
+          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o \
           $(OBJ)cmvtuniv.o $(OBJ)cmvtransf.o $(OBJ)cmvtgrowth.o $(OBJ)cmvtstpk.o \
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
@@ -189,4 +191,19 @@
 
 ##############################################################################
+cmvginit3d: $(EXE)cmvginit3d
+	echo $@ " done"
+$(EXE)cmvginit3d: $(OBJ)cmvginit3d.o $(LIBR) $(LIBG)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvginit3d.o $(MYLIB) 
+$(OBJ)cmvginit3d.o: cmvginit3d.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvginit3d.cc
+
+cmvginit3df: $(EXE)cmvginit3df
+	echo $@ " done"
+$(EXE)cmvginit3df: $(OBJ)cmvginit3df.o $(LIBR) $(LIBG4)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvginit3df.o $(MYLIB4) 
+$(OBJ)cmvginit3df.o: cmvginit3d.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -DGEN3D_FLOAT -o $@ cmvginit3d.cc
+
+##############################################################################
 cmvtintfun: $(EXE)cmvtintfun
 	echo $@ " done"
Index: /trunk/Cosmo/SimLSS/cmvginit3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3768)
+++ /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3768)
@@ -0,0 +1,618 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "sophyainit.h"
+#include "timing.h"
+#include "dvlist.h"
+#include "ntuple.h"
+#include "matharr.h"
+#include "randfmt.h"
+//#include "randr48.h"
+#include "srandgen.h"
+
+#include "constcosmo.h"
+#include "geneutils.h"
+#include "genefluct3d.h"
+
+void usage(void);
+void usage(void)
+{
+ cout<<"cmvobserv3d [...options...]"<<endl
+     <<" -a : auto init random seed (needed for multiple simul)"<<endl
+     <<" -G typevol: compute Pk(z=0) and apply growth factor in real space"<<endl
+     <<"       typevol=1 evolved with distance / observateur (def)"<<endl
+     <<"       typevol=2 evolved with distance to middle of Z planes"<<endl
+     <<"       else : no evol, spectrum Pk(z=z_median) for all cube (def)"<<endl
+     <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
+     <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
+     <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
+     <<"            if ny or dy <=0 take same value as for x"<<endl
+     <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
+     <<" -Z zref : redshift for the center of the simulation cube"<<endl
+     <<" -1 : compute 1D spectrum (default: no)"<<endl
+     <<" -2 : compute 2D spectrum (default: no)"<<endl
+     <<" -C : go back to cube in FT space and compute P(k) for check (def: do nothing)"<<endl
+     <<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl
+     <<"               (default sigmaR=1, R=8/h100 Mpc)"<<endl
+     <<" -v temp_file.ppf: generate los velocity cube"<<endl
+     <<"                   temporary save cube in temp_file.ppf"<<endl
+     <<" -n nent : fill control ntuple with nent entries (def: do not fill)"<<endl
+     <<" -O a,b : tell what you want to write (with the -W and -P options)"<<endl
+     <<"      a : write generated fourier cube (_k0)"<<endl
+     <<"      b : write real space cube dRho/Rho at z (_r0)"<<endl
+     <<"      a,b= 0 no write, 1 ppf write, 2 fits write, 3 ppf+fits write"<<endl
+     <<" -S : write cube slices in PPF format"<<endl
+     <<" -o root_name_out : root string for output file name (def: cmvobserv3d)"<<endl
+     <<" -V : compute variance from real space (for check, default: no)"<<endl
+     <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
+     <<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ SophyaInit();
+ InitTim();
+
+ //-----------------------------------------------------------------
+ // *** Survey definition
+ long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
+ //long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
+ //long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
+
+ // *** Cosmography definition   (WMAP)
+ unsigned short flat = 0;
+ double ob0 = 0.0444356;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double zref = 0.5;
+ double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
+
+ // *** 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;
+
+ // *** type de generation
+ int use_growth_factor = 0;
+ unsigned short nthread=0;
+ int filter_by_pixel = 1;
+
+ // *** What to do
+ bool comptransveloc = false;
+ string temporay_file = "cmvginit3d_tmp.ppf";
+ // compute: 1=1D spectra,  2=2D spectra,
+ //          4=recompute spectra after real space generation
+ unsigned short compdspec = 0;
+ long ntnent = 0;  // fill control NTuple
+ bool wslice = false;
+ bool compvarreal = false;
+ unsigned short whattowrt[2] = {1,1};
+ string rootnameout = "cmvobserv3d";
+
+ // --- Decodage arguments
+ if(narg>0) {
+   cout<<"\n--- Arguments: "<<endl;
+   for(int i=0;i<narg;i++) cout<<arg[i]<<" ";
+   cout<<endl;
+ }
+ system("date -u");
+
+ // --- Choix du generateur aleatoire (a faire ICI imperativement avant AutoInitRand)
+ FMTRandGen *RandGen = new FMTRandGen;
+ RandGen->SelectGaussianAlgo(C_Gaussian_RandLibSNorm);
+ RandomGeneratorInterface::SetGlobalRandGenP(RandGen);
+ 
+ // --- Decodage des arguments
+ char c;
+ while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:")) != -1) {
+  int nth = 0;
+  switch (c) {
+  case 'a' :
+    AutoInitRand(5);
+    break;
+  case 'G' :
+    use_growth_factor = atoi(optarg);
+    break;
+  case 'F' :
+    filter_by_pixel = atoi(optarg);
+    break;
+  case 'x' :
+    sscanf(optarg,"%ld,%lf",&nx,&dx);
+    break;
+  case 'y' :
+    sscanf(optarg,"%ld,%lf",&ny,&dy);
+    break;
+  case 'z' :
+    sscanf(optarg,"%ld,%lf",&nz,&dz);
+    break;
+  case 'Z' :
+    zref = atof(optarg);
+    break;
+  case 'v' :
+    comptransveloc = true;
+    temporay_file = optarg;
+    break;
+  case '1' :
+    compdspec |= 1;
+    break;
+  case '2' :
+    compdspec |= 2;
+    break;
+  case 'C' :
+    compdspec |= 4;
+    break;
+  case '8' :
+    sscanf(optarg,"%lf,%lf",&sigmaR,&R);
+    break;
+  case 'V' :
+    compvarreal = true;
+    break;
+  case 'n':
+    ntnent = atol(optarg);
+    if(ntnent<0) ntnent = 0;
+    break;
+  case 'O' :
+    sscanf(optarg,"%hu,%hu",&whattowrt[0],&whattowrt[1]);
+    break;
+  case 'S' :
+    wslice = true;
+    break;
+  case 'o' :
+    rootnameout = optarg;
+    break;
+  case 'T' :
+    nth = atoi(optarg);
+    nthread = (nth<1)? 0: nth;
+    break;
+  case 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+ try {
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Initialisation and parameters"<<endl;
+ 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<<"Filter by pixel = "<<filter_by_pixel<<endl;
+ if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl;
+ cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
+ cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
+ cout<<"wslice="<<wslice<<" what?="<<whattowrt[0]<<","<<whattowrt[1]<<endl;
+ cout<<"rootnameout="<<rootnameout<<endl;
+ ShowRandom();
+ cout<<"   First random is: "<<drand01()<<endl;
+
+ string tagobs = rootnameout + ".ppf";
+ cout<<"Writing result in "<<tagobs<<endl;
+ POutPersist posobs(tagobs);
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Cosmology"<<endl;
+
+ CosmoCalc univ(flat,true,zref+1.);
+ univ.SetInteg(perc,dzinc,dzmax,glorder);
+ univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.PrtInteg();
+ univ.Print();
+ double loscomref = univ.Dloscom(zref);
+ cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
+ univ.Print(zref);
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum"<<endl;
+
+ InitialSpectrum pkini(ns,as);
+
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+ //tf.SetNoOscEnv(2);
+
+ GrowthFactor growth(om0,ol0);
+ // GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
+ double growth_at_z = growth(zref);
+ cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
+
+ PkSpectrum0 pk0(pkini,tf);
+
+ PkSpectrumZ pkz(pk0,growth,zref);
+
+ //-----------------------------------------------------------------
+ pkz.SetZ(0.);
+ cout<<endl<<"\n--- Compute variance for top-hat R="<<R<<" at z="<<pkz.GetZ()<<endl;
+ VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
+ double kfind_th = varpk_th.FindMaximum(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(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(k1,k2);
+ cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
+
+ double normpkz = sigmaR*sigmaR/sr2;
+ pkz.SetScale(normpkz);
+ cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
+ {
+ Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
+ FuncToHisto(pkz,hpkz0,true);
+ posobs.PutObject(hpkz0,"hpkz0");
+ }
+
+ pkz.SetZ(zref);
+ {
+ Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
+ FuncToHisto(pkz,hpkz,true);
+ posobs.PutObject(hpkz,"hpkz");
+ }
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
+ VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
+ double kfind_int = varpk_int.FindMaximum(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(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(k1int,k2int);
+ cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
+
+ //-----------------------------------------------------------------
+
+ PrtTim(">>>> End of definition");
+
+ //-----------------------------------------------------------------
+ // Le cube et sa cosmlogie
+ // 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--- Initialisation de GeneFluct3D"<<endl;
+
+ GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
+ fluct3d.SetCosmology(univ);
+ fluct3d.SetGrowthFactor(growth);
+ fluct3d.SetObservator(zref,-nz/2.);
+ fluct3d.LosComRedshift(0.001,-1);
+ //TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();
+ //TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
+ cout<<endl; fluct3d.Print();
+
+ //-----------------------------------------------------------------
+ // Les bins et ranges en k pour les histos 1D et 2D
+ double dkmin = fluct3d.GetKincMin();
+ double knyqmax = fluct3d.GetKmax();
+ long nherr = long(knyqmax/dkmin+0.5);
+ cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
+
+ double dktmin = fluct3d.GetKTincMin();
+ double ktnyqmax = fluct3d.GetKTmax();
+ long nherrt = long(ktnyqmax/dktmin+0.5);
+ double dkzmin = fluct3d.GetKinc()[2];
+ double kznyqmax = fluct3d.GetKnyq()[2];
+ long nherrz = long(kznyqmax/dkzmin+0.5);
+ cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
+     <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
+
+ //-----------------------------------------------------------------
+ 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
+ double knyqmax_mod = 0.877*knyqmax;
+ ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
+ varpk_int.SetInteg(0.01,ldlkint,-1.,4);
+ double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
+ cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
+
+ PrtTim(">>>> End Initialisation of GeneFluct3D");
+
+ //-----------------------------------------------------------------
+ cout<<"\n--- Computing a realization in Fourier space"<<endl;
+ if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
+ cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
+ fluct3d.ComputeFourier(pkz);
+ fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
+ PrtTim(">>>> End Computing a realization in Fourier space");
+
+ HistoErr hpkgen(0.,knyqmax,nherr);
+ if(compdspec&1) {
+   cout<<"\n--- Checking realization spectra"<<endl;
+   hpkgen.ReCenterBin(); hpkgen.Zero(); hpkgen.Show();
+   fluct3d.ComputeSpectrum(hpkgen);
+   posobs.PutObject(hpkgen,"hpkgen");
+   PrtTim(">>>> End Checking realization 1D spectra");
+ }
+ Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+ if(compdspec&2) {
+   cout<<"\n--- Checking realization 2D spectra"<<endl;
+   hpkgen2.ReCenterBin(); hpkgen2.Zero(); hpkgen2.Show();
+   fluct3d.ComputeSpectrum2D(hpkgen2);
+   posobs.PutObject(hpkgen2,"hpkgen2");
+   PrtTim(">>>> End Checking realization 2D spectra");
+ }
+
+ //-----------------------------------------------------------------
+ if(filter_by_pixel!=0) {
+   cout<<"\n--- Computing convolution by pixel shape"<<endl;
+   fluct3d.FilterByPixel();
+   fluct3d.NTupleCheck(posobs,string("ntpkgenf"),ntnent);
+   PrtTim(">>>> End Computing convolution by pixel shape");
+ }
+
+ if(compdspec&1) {
+   HistoErr hpkgenf(hpkgen);
+   if(filter_by_pixel!=0) {
+     cout<<"\n--- Checking realization spectra"<<endl;
+     hpkgenf.Zero(); hpkgenf.Show();
+     fluct3d.ComputeSpectrum(hpkgenf);
+     PrtTim(">>>> End Checking realization 1D spectra");
+   }
+   posobs.PutObject(hpkgenf,"hpkgenf");
+ }
+ if(compdspec&2) {
+   Histo2DErr hpkgenf2(hpkgen2);
+   if(filter_by_pixel!=0) {
+     cout<<"\n--- Checking realization 2D spectra"<<endl;
+     hpkgenf2.Zero(); hpkgenf2.Show();
+     fluct3d.ComputeSpectrum2D(hpkgenf2);
+     PrtTim(">>>> End Checking realization 2D spectra");
+   }
+   posobs.PutObject(hpkgenf2,"hpkgenf2");
+ }
+
+ //-----------------------------------------------------------------
+ if(whattowrt[0]&1) {
+   tagobs = rootnameout + "_k.ppf";
+   fluct3d.WritePPF(tagobs,false);
+   PrtTim(">>>> End WritePPF");
+   temporay_file = tagobs;
+ }
+ if(whattowrt[0]&2) {
+   tagobs = "!" + rootnameout + "_k.fits";
+   fluct3d.WriteFits(tagobs);
+   PrtTim(">>>> End WriteFits");
+ }
+ if(comptransveloc && !(whattowrt[0]&1)) {
+   cout<<">>> Writing FT cube (for transv veloc.) in "<<temporay_file<<endl;
+   fluct3d.WritePPF(temporay_file,false);
+   PrtTim(">>>> End Writing FT cube");
+ }
+
+ //-----------------------------------------------------------------
+ cout<<"\n--- Computing a realization in real space"<<endl;
+ fluct3d.ComputeReal();
+ double rmin,rmax; fluct3d.MinMax(rmin,rmax);
+ cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+ fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
+ PrtTim(">>>> End Computing a realization in real space");
+
+ if(use_growth_factor>0) {
+   cout<<"\n--- Apply Growth factor"<<endl;
+   cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
+   fluct3d.ApplyGrowthFactor(use_growth_factor);
+   fluct3d.MinMax(rmin,rmax);
+   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
+   PrtTim(">>>> End Applying growth factor");
+ }
+
+ int_8 nm;
+ double rmref,rs2ref;
+ cout<<"\n--- Computing reference variance in real space"<<endl;
+ nm = fluct3d.MeanSigma2(rmref,rs2ref);
+ cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
+ PrtTim(">>>> End Computing reference variance in real space");
+
+ if(whattowrt[1]&1) {
+   tagobs = rootnameout + "_r.ppf";
+   fluct3d.WritePPF(tagobs,true);
+   PrtTim(">>>> End WritePPF");
+ }
+ if(whattowrt[1]&2) {
+     tagobs = "!" + rootnameout + "_r.fits";
+     fluct3d.WriteFits(tagobs);
+     PrtTim(">>>> End WriteFits");
+ }
+ if(wslice) {
+   tagobs = rootnameout + "_s_r.ppf";
+   fluct3d.WriteSlicePPF(tagobs);
+   PrtTim(">>>> End WriteSlicePPF");
+ }
+
+ //-----------------------------------------------------------------
+ if(compvarreal) {
+   cout<<"\n--- Check variance sigmaR in real space"<<endl;
+   double varr;
+   fluct3d.VarianceFrReal(R,varr);
+   cout<<"...Computed variance = "<<varr
+       <<" , Theorical variance at (z=0) = "<<pow(sigmaR,2.)
+       <<" , at (z="<<zref<<") = "<<pow(sigmaR*growth_at_z,2.)<<endl;
+   PrtTim(">>>> End Check variance sigmaR in real space");
+ }
+
+ //-----------------------------------------------------------------
+ if(compdspec&4) { // ATTENTION:  d_rho/rho ecrase
+   cout<<endl<<"\n--- ReComputing spectrum from real space"
+       <<"\n    Warning: REAL SPACE CUBE OVERWRITTEN"<<endl;
+   fluct3d.ReComputeFourier();
+   fluct3d.NTupleCheck(posobs,string("ntpkrec"),ntnent);
+   PrtTim(">>>> End ReComputing spectrum");
+
+   cout<<endl<<"\n--- Computing final spectrum"<<endl;
+   HistoErr hpkrecb(0.,knyqmax,nherr); hpkrecb.Zero();
+   hpkrecb.ReCenterBin(); hpkrecb.Show();
+   fluct3d.ComputeSpectrum(hpkrecb);
+   posobs.PutObject(hpkrecb,"hpkrecb");
+   PrtTim(">>>> End Computing final spectrum");
+
+   HistoErr hpkrec(hpkrecb);
+   if(filter_by_pixel!=0) {
+     cout<<endl<<"\n--- Computing final spectrum with pixel deconv."<<endl;
+     hpkrec.Zero();
+     fluct3d.ComputeSpectrum(hpkrec,0.,filter_by_pixel);
+     PrtTim(">>>> End Computing final spectrum with pixel deconv.");
+   }
+   posobs.PutObject(hpkrec,"hpkrec");
+
+   cout<<"\n--- Computing final 2D spectrum"<<endl;
+   Histo2DErr hpkrecb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkrecb2.ReCenterBin(); hpkrecb2.Zero(); hpkrecb2.Show();
+   fluct3d.ComputeSpectrum2D(hpkrecb2);
+   posobs.PutObject(hpkrecb2,"hpkrecb2");
+   PrtTim(">>>> End Computing final 2D spectrum");
+
+   Histo2DErr hpkrec2(hpkrecb2);
+   if(filter_by_pixel!=0) {
+     cout<<"\n--- Computing final 2D spectrum with pixel deconv."<<endl;
+     hpkrec2.Zero();
+     fluct3d.ComputeSpectrum2D(hpkrec2,0.,filter_by_pixel);
+     PrtTim(">>>> End Computing final 2D spectrum with pixel deconv.");
+   }
+   posobs.PutObject(hpkrec2,"hpkrec2");
+ }
+
+ //-----------------------------------------------------------------
+ //-----------------------------------------------------------------
+ //-----------------------------------------------------------------
+
+ if(comptransveloc) {
+   cout<<"\n\n\n"<<"---------------------------------------------\n"
+                 <<"--- Transverse velocity space computation ---\n"
+                 <<"---------------------------------------------\n"<<endl;
+
+   //-----------------------------------------------------------------
+   cout<<"\n--- Reading back the Pk(vec(k)) cube from "<<temporay_file<<endl;
+   fluct3d.ReadPPF(temporay_file);
+   PrtTim(">>>> End Reading the Pk(vec(k)) cube");
+
+   //-----------------------------------------------------------------
+   cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
+   fluct3d.ToVelTrans();
+   fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
+   PrtTim(">>>> End Modifying cube for Transverse velocity");
+
+   if(compdspec&1) {
+     cout<<"\n--- Checking realization spectra"<<endl;
+     HistoErr hpvgen(0.,knyqmax,nherr);
+     hpvgen.ReCenterBin(); hpvgen.Zero(); hpvgen.Show();
+     fluct3d.ComputeSpectrum(hpvgen);
+     posobs.PutObject(hpvgen,"hpvgen");
+     PrtTim(">>>> End Checking realization 1D spectra");
+   }
+   if(compdspec&2) {
+     cout<<"\n--- Checking realization 2D spectra"<<endl;
+     Histo2DErr hpvgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+     hpvgen2.ReCenterBin(); hpvgen2.Zero(); hpvgen2.Show();
+     fluct3d.ComputeSpectrum2D(hpvgen2);
+     posobs.PutObject(hpvgen2,"hpvgen2");
+     PrtTim(">>>> End Checking realization 2D spectra");
+   }
+
+    if(whattowrt[0]&1) {
+      tagobs = rootnameout + "_kv.ppf";
+      fluct3d.WritePPF(tagobs,false);
+      PrtTim(">>>> End WritePPF");
+    }
+    if(whattowrt[0]&2) {
+        tagobs = "!" + rootnameout + "_kv.fits";
+        fluct3d.WriteFits(tagobs);
+        PrtTim(">>>> End WriteFits");
+    }
+
+   //-----------------------------------------------------------------
+   cout<<"\n--- Computing a realization in real space for Transverse velocity"<<endl;
+   fluct3d.ComputeReal();
+   double rmin,rmax; fluct3d.MinMax(rmin,rmax);
+   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.NTupleCheck(posobs,string("nvreal"),ntnent);
+   PrtTim(">>>> End Computing a realization in real space");
+
+    if(use_growth_factor>0) {
+      cout<<"\n--- Apply Growth factor"<<endl;
+      cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
+      fluct3d.ApplyDerGrowthFactor(use_growth_factor);
+      fluct3d.MinMax(rmin,rmax);
+      cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+      fluct3d.NTupleCheck(posobs,string("nvgrow"),ntnent);
+      PrtTim(">>>> End Applying growth factor");
+    }
+
+    int_8 nmv;
+    double vmref,vs2ref;
+    cout<<"\n--- Computing reference variance in real space"<<endl;
+    nmv = fluct3d.MeanSigma2(vmref,vs2ref);
+    cout<<" vs2ref= "<<vs2ref<<" , vmref="<<vmref<<" ("<<nmv<<")"<<endl;
+    PrtTim(">>>> End Computing reference variance in real space");
+
+    if(whattowrt[1]&1) {
+      tagobs = rootnameout + "_rv.ppf";
+      fluct3d.WritePPF(tagobs,true);
+      PrtTim(">>>> End WritePPF");
+    }
+    if(whattowrt[1]&2) {
+        tagobs = "!" + rootnameout + "_rv.fits";
+        fluct3d.WriteFits(tagobs);
+        PrtTim(">>>> End WriteFits");
+    }
+    if(wslice) {
+      tagobs = rootnameout + "_s_rv.ppf";
+      fluct3d.WriteSlicePPF(tagobs);
+      PrtTim(">>>> End WriteSlicePPF");
+    }
+
+ }
+
+ //-----------------------------------------------------------------
+ {
+ DVList dvl;
+ dvl("Nx") = (int_4)fluct3d.Nx(); dvl("Ny") = (int_4)fluct3d.Ny(); dvl("Nz") = (int_4)fluct3d.Nz();
+ dvl("Dx") = fluct3d.Dx(); dvl("Dy") = fluct3d.Dy(); dvl("Dz") = fluct3d.Dz();
+ dvl("Z") = fluct3d.Zref(); dvl("Zpk") = fluct3d.ZrefPk();
+ dvl("D") = fluct3d.Dref(); dvl("Dpk") = fluct3d.DrefPk();
+ dvl("s8") = sigmaR;
+ dvl("Dtype") = (int_4)use_growth_factor; dvl("Pxfilt") = (int_4)filter_by_pixel; 
+ posobs.PutObject(dvl,"siminfo");
+ }
+
+ delete RandGen;
+ PrtTim(">>>> End Of Job");
+
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+ } catch (PException& exc) {
+   cerr<<"cmvobserv3d.cc catched PException"<<exc.Msg()<<endl;
+   return 77;
+ } catch (std::exception& sex) {
+   cerr << "cmvobserv3d.cc std::exception :" 
+        << (string)typeid(sex).name() << "\n msg= " 
+        << sex.what() << endl;
+   return 78;
+ } catch (...) {
+   cerr << "cmvobserv3d.cc catched unknown (...) exception  " << endl; 
+   return 79;
+ }
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+
+ return 0;
+}
+
Index: /trunk/Cosmo/SimLSS/cmvhshsns.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvhshsns.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/cmvhshsns.cc	(revision 3768)
@@ -3,4 +3,6 @@
 // > cmvhshsns -n -d 0.105 -g 2,0.105,0. -i 2,0.41,0. -t 50,180,0 -f 1420 1410 1430
 // > cmvhshsns -n -D -1 -d 10 -g 1 -i 10,10,0. -t 50,180,0 -f 1420 1410 1430
+// - Pittsburgh Sun 24 Nov
+//   > cmvhshsns -n -d 0.147 -D 0 -g 1 -i 8,0.147,61.06 -t 360,90,0 -f 1420 
 #include "sopnamsp.h"
 #include "machdefs.h"
@@ -28,4 +30,5 @@
   <<"  -D : +1 lobe du dipole approximation de Hertz"<<endl
   <<"       -1 le lobe du dipole est remplace par le sinc() du cylindre (etude E-W)"<<endl
+  <<"        0 lobe du dipole a 2 brins de longueur \"L\""<<endl
   <<"  -g N_g,D_g,Theta_g : regroupement de dipoles"<<endl
   <<"                       N_g<=1 pas de regroupement"<<endl
Index: /trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3768)
@@ -17,5 +17,4 @@
 #include "constcosmo.h"
 #include "cosmocalc.h"
-#include "schechter.h"
 #include "geneutils.h"
 #include "genefluct3d.h"
@@ -141,4 +140,5 @@
  string rootnameout = "cmvobserv3d";
 
+ // *** Number of entries in control NTuple
  unsigned long ntnent = 10000;  // 0 = do not fill NTuple
 
@@ -415,7 +415,7 @@
 
  GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
- fluct3d.SetObservator(zref,-nz/2.);
  fluct3d.SetCosmology(univ);
  fluct3d.SetGrowthFactor(growth);
+ fluct3d.SetObservator(zref,-nz/2.);
  fluct3d.LosComRedshift(0.001,-1);
  //TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();
Index: /trunk/Cosmo/SimLSS/cmvtgrowth.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/cmvtgrowth.cc	(revision 3768)
@@ -14,16 +14,15 @@
 
 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);
+void usage(void) {cout<<"cmvtgrowth z1,z2,nz [Omatter0,Lambda0]"<<endl;}
 
 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 z1=0., z2=10.;
+ int nz = 100;;
+ if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&z1,&z2,&nz);
+ if(nz<=0) nz = 100;
+ double dz = (z2-z1)/nz;
+ cout<<"z1="<<z1<<"  z2="<<z2<<"  nz="<<nz<<"  dz="<<dz<<endl;
 
  double om0=0.267804, ol0=0.73;
@@ -34,6 +33,6 @@
  cout<<"D1(z=0) = "<<growth(0.)<<endl;
 
- const int n = 2;
- const char *vname[n] = {"z","d1"};
+ const int n = 3;
+ const char *vname[n] = {"z","d1","d1dz"};
  NTuple nt(n,vname);
  double xnt[n];
@@ -41,4 +40,5 @@
    xnt[0]  = z;
    xnt[1]  = growth(z);
+   xnt[2]  = growth.DsDz(z,dz/2.);
    nt.Fill(xnt);
  }
@@ -59,5 +59,23 @@
 set cut z<5
 
+# --- growth
 n/plot nt.d1%z $cut ! "nsta connectpoints"
 n/plot nt.1./(1.+z)%z $cut ! "nsta connectpoints red same"
+
+# --- growth'/growth
+n/plot nt.d1dz/d1%z $cut ! "nsta connectpoints"
+n/plot nt.-1./(1.+z)%z $cut ! "nsta connectpoints red same"
+
+# --- d(growth)/dz
+zmin = 0.
+zmax = 10.
+set npt ${dv.size}
+dd = ($zmax-$zmin)/(${npt}-1.)
+exptovec dv nt d1
+c++exec for(int i=1;i<dv.Size(); i++) dv(i-1) = dv(i)-dv(i-1);
+
+n/plot nt.-1./pow(1.+z,2.)%z $cut ! "nsta connectpoints red"
+n/plot nt.d1dz%z $cut ! "nsta connectpoints same"
+n/plot dv.val/${dd}%${zmin}+(n+1)*$dd ! ! "nsta plusmarker5 green same"
+
 */
Index: /trunk/Cosmo/SimLSS/cosmocalc.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/cosmocalc.cc	(revision 3768)
@@ -21,5 +21,5 @@
 
 //----------------------------------------------------------
-// flat = 0 : no compute Otot from other
+// flat = 0 : do not set Otot=1, compute Otot from others
 //      = 1 : change Omatter to get Otot=1
 //      = 2 : change Olambda to get Otot=1
@@ -96,4 +96,15 @@
  if(_xspl) delete [] _xspl; _xspl = NULL;
  if(_yspl) delete [] _yspl; _yspl = NULL;
+}
+
+//----------------------------------------------------------
+// On change le zmax de la cosmologie
+void CosmoCalc::ChangeZmax(double zmax,double dzinc,double dzmax)
+{
+  _zmax = zmax;
+  if(dzinc<=0. ) dzinc = _dzinc;
+  if(dzmax<=0.)  dzmax = _dzmax;
+  cout<<"CosmoCalc::ChangeZmax: zmax="<<_zmax<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<endl;
+  SetInteg(_dperc,dzinc,_dzmax,_glorder);
 }
 
Index: /trunk/Cosmo/SimLSS/cosmocalc.h
===================================================================
--- /trunk/Cosmo/SimLSS/cosmocalc.h	(revision 3767)
+++ /trunk/Cosmo/SimLSS/cosmocalc.h	(revision 3768)
@@ -21,4 +21,5 @@
 
   double ZMax(void) { return _zmax;}
+  void ChangeZmax(double zmax,double dzinc=-1.,double dzmax=-1.);
   void SetInteg(double dperc=-1.,double dzinc=-1.,double dzmax=-1.,unsigned short order=4);
   void PrtInteg(void);
Index: /trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3768)
@@ -1,3 +1,2 @@
-#include "sopnamsp.h"
 #include "machdefs.h"
 #include <iostream>
@@ -88,4 +87,6 @@
  cosmo_ = NULL;
  growth_ = NULL;
+ good_dzinc_ = 0.01;
+ compute_pk_redsh_ref_ = -999.;
  redsh_ref_ = -999.;
  kredsh_ref_ = 0.;
@@ -262,8 +263,9 @@
 // 
 {
+ if(zinc<=0.) zinc = good_dzinc_;
  if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
 
- if(cosmo_ == NULL || redsh_ref_<0.) {
-   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
+ if(cosmo_ == NULL || growth_==NULL || redsh_ref_<0.) {
+   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator, Cosmology and Growth first";
    cout<<bla<<endl; throw ParmError(bla);
  }
@@ -275,4 +277,7 @@
  dlum_ref_ = cosmo_->Dlum(redsh_ref_);
  dang_ref_ = cosmo_->Dang(redsh_ref_);
+ h_ref = cosmo_->H(redsh_ref_);
+ growth_ref_ = (*growth_)(redsh_ref_);
+ dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
  nu_ref_   = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
  dnu_ref_  = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
@@ -282,4 +287,7 @@
        <<", nuref="<<nu_ref_ <<" GHz"
        <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
+       <<"   H="<<h_ref<<" km/s/Mpc"
+       <<", D="<<growth_ref_
+       <<", dD/dz="<<dsdz_growth_ref_<<endl
        <<"   dlosc="<<loscom_ref_<<" Mpc com"
        <<", dtrc="<<dtrc_ref_<<" Mpc com"
@@ -335,6 +343,5 @@
 
  // Fill the corresponding vectors for loscom and zred
- // Be shure to have one dlc <loscom_min and one >loscom_max
- if(zinc<=0.) zinc = 0.01;
+ // Be shure to have one dlc<loscom_min and one dlc>loscom_max
  double zmin = 0., dlcmin=0.;
  while(1) {
@@ -397,4 +404,24 @@
  }
 
+
+ // Compute the table for D'(z)/D(z) where D'(z)=dD/dEta (Eta is conformal time)
+ zredmin_dpsd_ = zred_[0];
+ zredmax_dpsd_ = zred_[zred_.size()-1];
+ long nptd = long(sqrt(Nx_*Nx_ + Ny_*Ny_ +Nz_*Nz_));
+ if(nptd<10) nptd = 10;
+ good_dzinc_ = (zredmax_dpsd_ - zredmin_dpsd_)/nptd;
+ if(lp_>0) cout<<"...good_dzinc changed to "<<good_dzinc_<<endl;
+ dpsdfrzred_.resize(0);
+ double dz = (zredmax_dpsd_ - zredmin_dpsd_)/(nptd-1);
+ int nmod = nptd/5; if(nmod==0) nmod = 1;
+ if(lp_>0) cout<<"...Compute table D'/D on "<<nptd<<" pts for z["
+               <<zredmin_dpsd_<<","<<zredmax_dpsd_<<"]"<<endl;
+ for(int i=0;i<nptd;i++) {
+   double z = zredmin_dpsd_ + i*dz;
+   double v = -cosmo_->H(z) * growth_->DsDz(z,good_dzinc_) / (*growth_)(z);
+   dpsdfrzred_.push_back(v);
+   if(lp_ && (i%nmod==0 || i==nptd-1)) cout<<"    z="<<z<<"  D'/D="<<v<<" km/s/Mpc"<<endl;
+ }
+
  return zred_.size();
 }
@@ -415,6 +442,9 @@
    fwrt.WriteKey("DKY",Dky_," Mpc^-1");
    fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
+   fwrt.WriteKey("ZREFPK",compute_pk_redsh_ref_," Pk computed redshift");
    fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
    fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
+   fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
+   fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
    fwrt.Write(R_);
  } catch (PThrowable & exc) {
@@ -440,4 +470,5 @@
    double dy = fimg.ReadKey("DY");
    double dz = fimg.ReadKey("DZ");
+   double pkzref = fimg.ReadKey("ZREFPK");
    double zref = fimg.ReadKey("ZREF");
    double kzref = fimg.ReadKey("KZREF");
@@ -447,4 +478,5 @@
    SetObservator(zref,kzref);
    array_allocated_ = true;
+   compute_pk_redsh_ref_ = pkzref;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -468,6 +500,9 @@
    R_.Info()["DY"] = (r_8)Dy_;
    R_.Info()["DZ"] = (r_8)Dz_;
+   R_.Info()["ZREFPK"] = (r_8)compute_pk_redsh_ref_;
    R_.Info()["ZREF"] = (r_8)redsh_ref_;
    R_.Info()["KZREF"] = (r_8)kredsh_ref_;
+   R_.Info()["Growth"] = (r_8)Dref();
+   R_.Info()["GrowthPk"] = (r_8)DrefPk();
    POutPersist pos(cfname.c_str());
    if(write_real) pos << PPFNameTag("rgen")  << R_;
@@ -506,4 +541,5 @@
    r_8 dy = R_.Info()["DY"];
    r_8 dz = R_.Info()["DZ"];
+   r_8 pkzref = R_.Info()["ZREFPK"];
    r_8 zref = R_.Info()["ZREF"];
    r_8 kzref = R_.Info()["KZREF"];
@@ -512,4 +548,5 @@
    SetObservator(zref,kzref);
    array_allocated_ = true;
+   compute_pk_redsh_ref_ = pkzref;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -630,11 +667,12 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
+void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
 // cf ComputeFourier() mais avec autre methode de realisation du spectre
 //    (attention on fait une fft pour realiser le spectre)
 {
-
+  compute_pk_redsh_ref_ = pk_at_z.GetZ();
  // --- realisation d'un tableau de tirage gaussiens
- if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl;
+  if(lp_>0) cout<<"--- ComputeFourier0 at z="<<compute_pk_redsh_ref_
+                <<": before gaussian filling ---"<<endl;
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
@@ -651,14 +689,12 @@
  if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
  T_(0,0,0) = complex<GEN3D_TYPE>(0.);  // on coupe le continue et on l'initialise
- long lmod = Nx_/10; if(lmod<1) lmod=1;
+ long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   double kx = Kx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;  kz *= kz;
+       double kz = Kz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -680,5 +716,5 @@
 
 //-------------------------------------------------------
-void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
+void GeneFluct3D::ComputeFourier(PkSpectrumZ& pk_at_z)
 // Calcule une realisation du spectre "pk_at_z"
 // Attention: dans TArray le premier indice varie le + vite
@@ -692,17 +728,16 @@
  // --- RaZ du tableau
  T_ = complex<GEN3D_TYPE>(0.);
+ compute_pk_redsh_ref_ = pk_at_z.GetZ();
 
  // --- On remplit avec une realisation
- if(lp_>0) cout<<"--- ComputeFourier ---"<<endl;
- long lmod = Nx_/10; if(lmod<1) lmod=1;
+ if(lp_>0) cout<<"--- ComputeFourier at z="<<compute_pk_redsh_ref_<<" ---"<<endl;
+ long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
-   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   double kx = Kx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;  kz *= kz;
+       double kz = Kz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -728,4 +763,6 @@
 // Take into account the real and complexe conjugate coefficients
 // because we want a realization of a real data in real space
+// On ecrit que: P(k_x,k_y,k_z) = P(-k_x,-k_y,-k_z)
+//               avec k_x = i, -k_x = N_x - i  etc...
 {
  if(lp_>1) cout<<"...managing coefficients"<<endl;
@@ -800,5 +837,5 @@
  if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
 
- if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
+ if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-nreal<<endl;
 
  return nreal+nconj1+nconj2;
@@ -840,13 +877,11 @@
 
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_ *Dx_/2;
+   double kx = Kx(i) *Dx_/2;
    double pk_x = pixelfilter(kx);
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_ *Dy_/2;
+     double ky = Ky(j) *Dy_/2;
      double pk_y = pixelfilter(ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_ *Dz_/2;
+       double kz = Kz(l) *Dz_/2;
        double pk_z =  pixelfilter(kz);
        T_(l,j,i) *= pk_x*pk_y*pk_z;
@@ -855,4 +890,30 @@
  }
 
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ToVelTrans(void)
+{
+  cout<<"-------------------- TO BE VERIFIED ToVelTrans"<<endl;
+
+ double zpk = compute_pk_redsh_ref_;
+ double dpsd = -cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+ if(lp_>0) cout<<"--- ToVelTrans --- at z="<<zpk<<", D'/D="<<dpsd
+               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
+ check_array_alloc();
+
+ for(long i=0;i<Nx_;i++) {
+   double kx = Kx(i);
+   for(long j=0;j<Ny_;j++) {
+     double ky = Ky(j);
+     double kt2 = kx*kx + ky*ky;
+     for(long l=0;l<NCz_;l++) {
+       double kz = Kz(l);
+       double k2 = kt2 + kz*kz;
+       if(k2<=0.) continue;
+       T_(l,j,i) *= complex<double>(0.,dpsd*kz/k2);
+     }
+   }
+ }
 }
 
@@ -879,7 +940,5 @@
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
- //unsigned short ok;
-
- //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
+
  for(long i=0;i<Nx_;i++) {
    double dx2 = DXcom(i); dx2 *= dx2;
@@ -890,9 +949,6 @@
        if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
          else dz = fabs(dz); // tous les plans Z au meme redshift
-       double z = interpinv(dz);
-       //CHECK: hgr.Add(z);
-       double dzgr = (*growth_)(z);   // interpolation par morceau
-       //double dzgr = growth_->Linear(z,ok);  // interpolation lineaire
-       //double dzgr = growth_->Parab(z,ok);  // interpolation parabolique
+       double z = interpinv(dz);   // interpolation par morceau
+       double dzgr = (*growth_)(z);
        int_8 ip = IndexR(i,j,l);
        data_[ip] *= dzgr;
@@ -901,5 +957,45 @@
  }
 
- //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
+}
+
+//-------------------------------------------------------------------
+void GeneFluct3D::ApplyDerGrowthFactor(int type_evol)
+// Apply Conformal derivative of Growth to real space for transverse velocity cube
+{
+  cout<<"-------------------- TO BE VERIFIED ApplyDerGrowthFactor"<<endl;
+
+ if(lp_>0) cout<<"--- ApplyDerGrowthFactor: evol="<<type_evol<<endl;
+ check_array_alloc();
+
+ if(growth_ == NULL) {
+   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+ if(type_evol<1 || type_evol>2) {
+   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ double zpk = compute_pk_redsh_ref_;
+ double dpsd_orig = - cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+
+ InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+ InterpFunc interdpd(zredmin_dpsd_,zredmax_dpsd_,dpsdfrzred_);
+
+ for(long i=0;i<Nx_;i++) {
+   double dx2 = DXcom(i); dx2 *= dx2;
+   for(long j=0;j<Ny_;j++) {
+     double dy2 = DYcom(j); dy2 *= dy2;
+     for(long l=0;l<Nz_;l++) {
+       double dz = DZcom(l);
+       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
+         else dz = fabs(dz); // tous les plans Z au meme redshift
+       double z = interpinv(dz);   // interpolation par morceau
+       double dpsd = interdpd(z);
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] *= dpsd / dpsd_orig;
+     }
+   }
+ }
 
 }
@@ -948,11 +1044,9 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
+   double kx = Kx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i));
@@ -980,12 +1074,10 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;  kx *= kx;
+   double kx = Kx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;  ky *= ky;
+     double ky = Ky(j);  ky *= ky;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double pk = MODULE2(T_(l,j,i));
        herr.Add(kt,kz,pk);
@@ -1026,15 +1118,13 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;
+   double kx = Kx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;
+     double ky = Ky(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
      ky *= ky; fy *= fy;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
@@ -1073,16 +1163,14 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   long ii = (i>Nx_/2) ? Nx_-i : i;
-   double kx = ii*Dkx_;
+   double kx = Kx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     long jj = (j>Ny_/2) ? Ny_-j : j;
-     double ky = jj*Dky_;
+     double ky = Ky(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
      ky *= ky; fy *= fy;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;
+       double kz = Kz(l);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
        double fz = (pixcor) ? vfz(l): 1.;
Index: /trunk/Cosmo/SimLSS/genefluct3d.h
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3767)
+++ /trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3768)
@@ -18,9 +18,9 @@
 #include "cosmocalc.h"
 #include "pkspectrum.h"
+#include "schechter.h"
 
 #define WITH_FFTW_THREAD
 
-//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE #define GEN3D_FLOAT
-
+//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -DGEN3D_FLOAT ...")
 #if defined(GEN3D_FLOAT)
 #define GEN3D_TYPE r_4
@@ -77,7 +77,10 @@
   vector<long> GetNpix(void) {return N_;}
   int_8 NPix(void) {return NRtot_;}
-  long GetNx(void) {return Nx_;}
-  long GetNy(void) {return Ny_;}
-  long GetNz(void) {return Nz_;}
+  long Nx(void) {return Nx_;}
+  long Ny(void) {return Ny_;}
+  long Nz(void) {return Nz_;}
+  double Dx(void) {return Dx_;}
+  double Dy(void) {return Dy_;}
+  double Dz(void) {return Dz_;}
 
   // Return |K_i| module relative to pixel indices
@@ -101,10 +104,18 @@
   double GetKTincMax(void) {return max(Dk_[0],Dk_[1]);}
 
-  void ComputeFourier0(GenericFunc& pk_at_z);
-  void ComputeFourier(GenericFunc& pk_at_z);
+  double Zref(void) {return redsh_ref_;}
+  double ZrefPk(void) {return compute_pk_redsh_ref_;}
+  double Dref(void) {return growth_ref_;}
+  double DrefPk(void)
+    {return (growth_!=NULL && compute_pk_redsh_ref_>=0.) ? (*growth_)(compute_pk_redsh_ref_): -999.;}
+
+  void ComputeFourier0(PkSpectrumZ& pk_at_z);
+  void ComputeFourier(PkSpectrumZ& pk_at_z);
   void FilterByPixel(void);
+  void ToVelTrans(void);
 
   void ComputeReal(void);
   void ApplyGrowthFactor(int type_evol=1);
+  void ApplyDerGrowthFactor(int type_evol=1);
 
   void ReComputeFourier(void);
@@ -188,8 +199,12 @@
   CosmoCalc *cosmo_;
   GrowthFactor *growth_;
+  double good_dzinc_;
   double redsh_ref_,kredsh_ref_,dred_ref_;
-  double loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
+  double compute_pk_redsh_ref_;
+  double h_ref, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
+  double growth_ref_, dsdz_growth_ref_;
   double nu_ref_, dnu_ref_ ;
   double xobs_[3];
+
   double loscom_min_, loscom_max_;
   vector<double> zred_, loscom_;
@@ -197,4 +212,7 @@
   vector<double> loscom2zred_;
 
+  double zredmin_dpsd_, zredmax_dpsd_;
+  vector<double> dpsdfrzred_;
+
 };
 
Index: /trunk/Cosmo/SimLSS/pkspectrum.cc
===================================================================
--- /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3767)
+++ /trunk/Cosmo/SimLSS/pkspectrum.cc	(revision 3768)
@@ -399,4 +399,36 @@
 }
 
+double GrowthFactor::DsDz(double z,double dzinc)
+// y-y0 = a*(x-x0)^2 + b*(x-x0)
+// dy = a*dx^2 + b*dx   avec  dx = x-x0
+// dy'(dx) = 2*a*dx + b -> pour x=x0 on a dy'(0) = b
+{
+ if(z<0. || dzinc<=0.) {
+   cout<<"GrowthFactor::DsDz: z<0 or dzinc<=0. !"<<endl;
+   throw ParmError("GrowthFactor::DsDz: z<0 or dzinc<=0. !");
+ }
+
+ double z1, z2;
+ if(z>dzinc/2.) {
+   // cas ou z est suffisamment loin de zero
+   // resolution avec 2 points encadrant x0=z
+   z1 = z - dzinc; if(z1<0.) z1 = 0.;
+   z2 = z + dzinc;
+ } else {
+   // cas ou z est proche de zero
+   // resolution avec 2 points au dessus, x0=z1
+   z1 = z + dzinc;
+   z2 = z + 2.*dzinc;
+ }
+
+ double gz  = (*this)(z);
+ double dgz1 = (*this)(z1) - gz;
+ double dgz2 = (*this)(z2) - gz;
+
+ z1 -= z;
+ z2 -= z;
+ return (dgz2*z1*z1 - dgz1*z2*z2)/(z1*z2*(z1-z2));
+}
+
 void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
 {
Index: /trunk/Cosmo/SimLSS/pkspectrum.h
===================================================================
--- /trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3767)
+++ /trunk/Cosmo/SimLSS/pkspectrum.h	(revision 3768)
@@ -77,4 +77,5 @@
   virtual ~GrowthFactor(void);
   virtual double operator() (double z);
+  virtual double DsDz(double z,double dzinc=0.01);
   void SetParTo(double OmegaMatter0,double OmegaLambda0);
   bool SetParTo(double OmegaMatter0);
