Index: /trunk/Cosmo/SimLSS/Makefile
===================================================================
--- /trunk/Cosmo/SimLSS/Makefile	(revision 3769)
+++ /trunk/Cosmo/SimLSS/Makefile	(revision 3770)
@@ -18,5 +18,5 @@
 PROGS = \
        $(EXE)cmvobserv3d $(EXE)cmvobserv3df \
-       $(EXE)cmvginit3d $(EXE)cmvginit3df \
+       $(EXE)cmvginit3d $(EXE)cmvginit3df $(EXE)cmvrvloscorf \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
@@ -26,5 +26,5 @@
 PROGSOBJ = \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvobserv3df.o \
-          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o \
+          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o $(OBJ)cmvrvloscorf.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 \
@@ -206,4 +206,12 @@
 
 ##############################################################################
+cmvrvloscorf: $(EXE)cmvrvloscorf
+	echo $@ " done"
+$(EXE)cmvrvloscorf: $(OBJ)cmvrvloscorf.o $(LIBR) $(LIBG4)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvrvloscorf.o $(MYLIB4) 
+$(OBJ)cmvrvloscorf.o: cmvrvloscor.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -DGEN3D_FLOAT -o $@ cmvrvloscor.cc
+
+##############################################################################
 cmvtintfun: $(EXE)cmvtintfun
 	echo $@ " done"
Index: /trunk/Cosmo/SimLSS/cmvginit3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3769)
+++ /trunk/Cosmo/SimLSS/cmvginit3d.cc	(revision 3770)
@@ -509,5 +509,5 @@
    //-----------------------------------------------------------------
    cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
-   fluct3d.ToVelTrans();
+   fluct3d.ToVelLoS();
    fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
    PrtTim(">>>> End Modifying cube for Transverse velocity");
@@ -601,13 +601,13 @@
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
  } catch (PException& exc) {
-   cerr<<"cmvobserv3d.cc catched PException"<<exc.Msg()<<endl;
+   cerr<<"cmvginit3d.cc catched PException"<<exc.Msg()<<endl;
    return 77;
  } catch (std::exception& sex) {
-   cerr << "cmvobserv3d.cc std::exception :" 
+   cerr << "cmvginit3d.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    return 78;
  } catch (...) {
-   cerr << "cmvobserv3d.cc catched unknown (...) exception  " << endl; 
+   cerr << "cmvginit3d.cc catched unknown (...) exception  " << endl; 
    return 79;
  }
Index: /trunk/Cosmo/SimLSS/cmvrvloscor.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvrvloscor.cc	(revision 3770)
+++ /trunk/Cosmo/SimLSS/cmvrvloscor.cc	(revision 3770)
@@ -0,0 +1,152 @@
+#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 "fabtcolread.h"
+
+#include "constcosmo.h"
+#include "geneutils.h"
+#include "genefluct3d.h"
+
+void usage(void);
+void usage(void)
+{
+  cout<<"cmvrvloscor rho.fits vlos.fits"<<endl;
+}
+
+int main(int narg,char *arg[])
+{
+ int nthread = 1;
+ int_8 ntfill = 100000;
+ SophyaInit();
+ InitTim();
+
+ if(narg<=2) {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<<"> read rho: "<<arg[1]<<endl;
+ FitsImg3DRead f3dr(arg[1],0,5);
+ cout<<"> read vlos: "<<arg[2]<<endl;
+ FitsImg3DRead f3dv(arg[2],0,5);
+ long naxis1 = f3dr.ReadKeyL("NAXIS1");
+ long naxis2 = f3dr.ReadKeyL("NAXIS2");
+ long naxis3 = f3dr.ReadKeyL("NAXIS3");
+ cout<<"Naxis: 1="<<naxis1<<" 2="<<naxis2<<" 3="<<naxis3<<endl;
+ long nx = f3dr.ReadKeyL("Nx");
+ long ny = f3dr.ReadKeyL("Ny");;
+ long nz = f3dr.ReadKeyL("Nz");;
+ cout<<"N: x="<<nx<<" y="<<ny<<" z="<<nz<<endl;
+ double dx = f3dr.ReadKey("Dx");
+ double dy = f3dr.ReadKey("Dy");
+ double dz = f3dr.ReadKey("Dz");
+ cout<<"D: x="<<dx<<" y="<<dy<<" z="<<dz<<endl;
+ double zref = f3dr.ReadKey("ZREF");
+ double href = f3dr.ReadKey("HREF");
+ cout<<"Zref="<<zref<<" Href="<<href<<endl;
+
+ POutPersist pos("cmvrvloscor.ppf");
+
+ GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
+ fluct3d.Print();
+ TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
+ rgen = 0.;
+ TVector<r_8> R(nz), Rdis(nz);
+
+ const int nvar = 3;
+ const char *vname[nvar] = {"l","ll","d"};
+ NTuple nt(nvar,vname);
+ r_4 xnt[nvar];
+ int_8 ntmod = fluct3d.NPix()/ntfill; if(ntmod==0) ntmod = 1;
+
+ cout<<"> filling redshift distorted cube"<<endl;
+ int_8 nf = 0;
+ for(int i=0;i<nx;i++) {
+   if(i%(nx/10)==0) cout<<"i="<<i<<endl;
+ for(int j=0;j<ny;j++) {
+   for(int l=0;l<nz;l++) R(l) = f3dr.Read(l,j,i);
+   Rdis = 0.;
+   for(int l=0;l<nz;l++) {
+     double d = (1.+zref) / href * f3dv.Read(l,j,i);
+     int ll = int((double)l + d/dz + 0.5);
+     if(ll<0 || ll>=nz) continue;
+     Rdis(ll) += R(l);
+     if(ntfill>0 && nf%ntmod==0)
+       {xnt[0]=l; xnt[1]=ll; xnt[2]=d; nt.Fill(xnt);}
+     nf++;
+   }
+   for(int l=0;l<nz;l++) rgen(l,j,i) += Rdis(l);
+ }
+ }
+ PrtTim(">>>> End filling redshift distorted cube");
+ pos.PutObject(nt,"nt");
+
+ fluct3d.ReComputeFourier();
+ PrtTim(">>>> End ReComputing spectrum");
+
+ cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
+ 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;
+ HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
+ hpkrec.ReCenterBin(); hpkrec.Show();
+ fluct3d.ComputeSpectrum(hpkrec);
+ pos.PutObject(hpkrec,"hpkrec");
+ PrtTim(">>>> End Computing final spectrum");
+
+ cout<<"\n--- Computing final 2D spectrum"<<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;
+ Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+ hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
+ fluct3d.ComputeSpectrum2D(hpkrec2);
+ pos.PutObject(hpkrec2,"hpkrec2");
+ PrtTim(">>>> End Computing final 2D spectrum");
+
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+ } catch (PException& exc) {
+   cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
+   return 77;
+ } catch (std::exception& sex) {
+   cerr << "cmvrvloscor.cc std::exception :" 
+        << (string)typeid(sex).name() << "\n msg= " 
+        << sex.what() << endl;
+   return 78;
+ } catch (...) {
+   cerr << "cmvrvloscor.cc catched unknown (...) exception  " << endl; 
+   return 79;
+ }
+ //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
+
+ return 0;
+}
+
+/*
+openppf cmvrvloscor.ppf
+
+n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
+
+imag hpkrec2
+
+# proj selon kT (black), selon kZ (red)
+n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta cpts"
+
+ */
Index: /trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3769)
+++ /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3770)
@@ -92,4 +92,5 @@
  kredsh_ref_ = 0.;
  dred_ref_ = -999.;
+ h_ref_ = -999.;
  loscom_ref_ = -999.;
  dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
@@ -277,5 +278,5 @@
  dlum_ref_ = cosmo_->Dlum(redsh_ref_);
  dang_ref_ = cosmo_->Dang(redsh_ref_);
- h_ref = cosmo_->H(redsh_ref_);
+ h_ref_ = cosmo_->H(redsh_ref_);
  growth_ref_ = (*growth_)(redsh_ref_);
  dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
@@ -287,5 +288,5 @@
        <<", nuref="<<nu_ref_ <<" GHz"
        <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
-       <<"   H="<<h_ref<<" km/s/Mpc"
+       <<"   H="<<h_ref_<<" km/s/Mpc"
        <<", D="<<growth_ref_
        <<", dD/dz="<<dsdz_growth_ref_<<endl
@@ -445,4 +446,5 @@
    fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
    fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
+   fwrt.WriteKey("HREF",h_ref_," ref H(z)");
    fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
    fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
@@ -503,4 +505,5 @@
    R_.Info()["ZREF"] = (r_8)redsh_ref_;
    R_.Info()["KZREF"] = (r_8)kredsh_ref_;
+   R_.Info()["HREF"] = (r_8)h_ref_;
    R_.Info()["Growth"] = (r_8)Dref();
    R_.Info()["GrowthPk"] = (r_8)DrefPk();
@@ -691,10 +694,10 @@
  long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   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;
+   double kx = AbsKx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);  kz *= kz;
+       double kz = AbsKz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -734,10 +737,10 @@
  long lmod = Nx_/20; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
-   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;
+   double kx = AbsKx(i);  kx *= kx;
+   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<Kx(i)<<"\t pk="<<pk_at_z(kx)<<endl;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);  kz *= kz;
+       double kz = AbsKz(l);  kz *= kz;
        if(i==0 && j==0 && l==0) continue; // Suppression du continu
        double k = sqrt(kx+ky+kz);
@@ -877,11 +880,11 @@
 
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i) *Dx_/2;
+   double kx = AbsKx(i) *Dx_/2;
    double pk_x = pixelfilter(kx);
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j) *Dy_/2;
+     double ky = AbsKy(j) *Dy_/2;
      double pk_y = pixelfilter(ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l) *Dz_/2;
+       double kz = AbsKz(l) *Dz_/2;
        double pk_z =  pixelfilter(kz);
        T_(l,j,i) *= pk_x*pk_y*pk_z;
@@ -893,11 +896,9 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToVelTrans(void)
-{
-  cout<<"-------------------- TO BE VERIFIED ToVelTrans"<<endl;
-
+void GeneFluct3D::ToVelLoS(void)
+{
  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
+ if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", D'/D="<<dpsd
                <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
  check_array_alloc();
@@ -963,6 +964,4 @@
 // 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();
@@ -1044,9 +1043,9 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
+   double kx = AbsKx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i));
@@ -1074,10 +1073,10 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);  kx *= kx;
+   double kx = AbsKx(i);  kx *= kx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);  ky *= ky;
+     double ky = AbsKy(j);  ky *= ky;
      double kt = sqrt(kx+ky);
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double pk = MODULE2(T_(l,j,i));
        herr.Add(kt,kz,pk);
@@ -1118,13 +1117,13 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);
+   double kx = AbsKx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);
+     double ky = AbsKy(j);
      double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
      ky *= ky; fy *= fy;
      for(long l=0;l<NCz_;l++) {
-       double kz = Kz(l);
+       double kz = AbsKz(l);
        double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i)) - sigma2;
@@ -1163,14 +1162,14 @@
  // Attention a l'ordre
  for(long i=0;i<Nx_;i++) {
-   double kx = Kx(i);
+   double kx = AbsKx(i);
    double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
    kx *= kx; fx *= fx;
    for(long j=0;j<Ny_;j++) {
-     double ky = Ky(j);
+     double ky = AbsKy(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 = Kz(l);
+       double kz = AbsKz(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 3769)
+++ /trunk/Cosmo/SimLSS/genefluct3d.h	(revision 3770)
@@ -22,5 +22,5 @@
 #define WITH_FFTW_THREAD
 
-//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -DGEN3D_FLOAT ...")
+//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -c -DGEN3D_FLOAT ...")
 #if defined(GEN3D_FLOAT)
 #define GEN3D_TYPE r_4
@@ -85,6 +85,9 @@
 
   // Return |K_i| module relative to pixel indices
-  inline r_8 Kx(long i) {long ii=(i>Nx_/2)? Nx_-i :i; return ii*Dkx_;}
-  inline r_8 Ky(long j) {long jj=(j>Ny_/2)? Ny_-j :j; return jj*Dky_;}
+  inline r_8 AbsKx(long i) {long ii=(i>Nx_/2)? Nx_-i :i; return ii*Dkx_;}
+  inline r_8 AbsKy(long j) {long jj=(j>Ny_/2)? Ny_-j :j; return jj*Dky_;}
+  inline r_8 AbsKz(long l) {return l*Dkz_;}
+  inline r_8 Kx(long i) {long ii=(i>Nx_/2)? i-Nx_ :i; return ii*Dkx_;}
+  inline r_8 Ky(long j) {long jj=(j>Ny_/2)? j-Ny_ :j; return jj*Dky_;}
   inline r_8 Kz(long l) {return l*Dkz_;}
 
@@ -109,9 +112,10 @@
   double DrefPk(void)
     {return (growth_!=NULL && compute_pk_redsh_ref_>=0.) ? (*growth_)(compute_pk_redsh_ref_): -999.;}
+  double Href(void) {return h_ref_;}
 
   void ComputeFourier0(PkSpectrumZ& pk_at_z);
   void ComputeFourier(PkSpectrumZ& pk_at_z);
   void FilterByPixel(void);
-  void ToVelTrans(void);
+  void ToVelLoS(void);
 
   void ComputeReal(void);
@@ -202,5 +206,5 @@
   double redsh_ref_,kredsh_ref_,dred_ref_;
   double compute_pk_redsh_ref_;
-  double h_ref, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
+  double h_ref_, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
   double growth_ref_, dsdz_growth_ref_;
   double nu_ref_, dnu_ref_ ;
