Index: /trunk/Cosmo/SimLSS/Makefile
===================================================================
--- /trunk/Cosmo/SimLSS/Makefile	(revision 3789)
+++ /trunk/Cosmo/SimLSS/Makefile	(revision 3790)
@@ -18,18 +18,20 @@
 PROGS = \
        $(EXE)cmvobserv3d $(EXE)cmvobserv3df \
-       $(EXE)cmvginit3d $(EXE)cmvginit3df $(EXE)cmvrvloscorf \
+       $(EXE)cmvginit3d $(EXE)cmvginit3df $(EXE)cmvrvloscor $(EXE)cmvrvloscorf \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
        $(EXE)cmvtintfun $(EXE)cmvconcherr $(EXE)cmvtinterp $(EXE)cmvtstagn \
-       $(EXE)cmvschdist $(EXE)cmvhshsns $(EXE)cmvhshsew
+       $(EXE)cmvschdist $(EXE)cmvhshsns $(EXE)cmvhshsew \
+       $(EXE)cmvgsm2sph $(EXE)cmvgsm2cube
  
 PROGSOBJ = \
           $(OBJ)cmvobserv3d.o $(OBJ)cmvobserv3df.o \
-          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o $(OBJ)cmvrvloscorf.o \
+          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o $(OBJ)cmvrvloscor.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 \
           $(OBJ)cmvtintfun.o $(OBJ)cmvtinterp.o \
           $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o $(OBJ)cmvtstagn.o $(OBJ)cmvschdist.o \
-          $(OBJ)cmvhshsns.o $(OBJ)cmvhshsew.o
+          $(OBJ)cmvhshsns.o $(OBJ)cmvhshsew.o \
+          $(OBJ)cmvgsm2sph.o $(OBJ)cmvgsm2cube.o
  
 LIBROBJ = \
@@ -57,8 +59,8 @@
 #---- Les programmes utilisant des librairies non standard SOPHYA
 SPROGS = \
-       $(EXE)cmvfitpk
+       $(EXE)cmvfitpk $(EXE)cmvgsmcfit
 
 SPROGSOBJ = \
-       $(OBJ)cmvfitpk.o
+       $(OBJ)cmvfitpk.o $(OBJ)cmvgsmcfit.o
 
 #----
@@ -206,4 +208,11 @@
 
 ##############################################################################
+cmvrvloscor: $(EXE)cmvrvloscor
+	echo $@ " done"
+$(EXE)cmvrvloscor: $(OBJ)cmvrvloscor.o $(LIBR) $(LIBG)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvrvloscorf.o $(MYLIB) 
+$(OBJ)cmvrvloscor.o: cmvrvloscor.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvrvloscor.cc
+
 cmvrvloscorf: $(EXE)cmvrvloscorf
 	echo $@ " done"
@@ -298,4 +307,26 @@
 
 ##############################################################################
+cmvgsm2sph: $(EXE)cmvgsm2sph
+	echo $@ " done"
+$(EXE)cmvgsm2sph: $(OBJ)cmvgsm2sph.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvgsm2sph.o $(MYLIB) 
+$(OBJ)cmvgsm2sph.o: cmvgsm2sph.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvgsm2sph.cc
+
+cmvgsm2cube: $(EXE)cmvgsm2cube
+	echo $@ " done"
+$(EXE)cmvgsm2cube: $(OBJ)cmvgsm2cube.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvgsm2cube.o $(MYLIB) 
+$(OBJ)cmvgsm2cube.o: cmvgsm2cube.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvgsm2cube.cc
+
+cmvgsmcfit: $(EXE)cmvgsmcfit
+	echo $@ " done"
+$(EXE)cmvgsmcfit: $(OBJ)cmvgsmcfit.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvgsmcfit.o $(MYLIB) -lMinuit2
+$(OBJ)cmvgsmcfit.o: cmvgsmcfit.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvgsmcfit.cc
+
+##############################################################################
 cmvhshsns: $(EXE)cmvhshsns
 	echo $@ " done"
Index: /trunk/Cosmo/SimLSS/cmvgsm2cube.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvgsm2cube.cc	(revision 3790)
+++ /trunk/Cosmo/SimLSS/cmvgsm2cube.cc	(revision 3790)
@@ -0,0 +1,107 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream>
+#include <string.h>
+#include <unistd.h>
+#include <vector>
+
+#include "fioarr.h"
+#include "skymap.h"
+
+// set fname = ( `cat gsm820948.list | awk '{print $3}' | sed 's/data/ppf/'` )
+// cmvgsm2cube -o gsm820948_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
+
+void usage(void)
+{
+cout
+<<"cmvgsm2cube [options] liste_des_spheres_healpix"<<endl
+<<" -o outname.ppf : ppf cube output file name (def=cmvgsm2cube.ppf)"<<endl
+<<" -p nphi,dphi,phi0 : phi0 -> phi0+nphi*dphi (degres)"<<endl
+<<" -t ntheta,dtheta,theta0 : theta0 -> theta0+ntheta*dtheta (degres)"<<endl
+<<endl;
+}
+
+
+int main(int narg,char *arg[])
+{
+ string ppfrn = "cmvgsm2cube.ppf";
+ int nphi=0, ntheta=0;
+ double dphi=-1., phi0=0.;
+ double dtheta=-1., theta0=0.;
+
+ if(narg<=1) {usage(); return -1;}
+
+// Read arguments
+ char c;
+ while ((c = getopt(narg, arg, "ho:p:t:")) != -1) {
+  switch (c) {
+  case 'o' :
+    ppfrn = optarg;
+    break;
+  case 'p' :
+    sscanf(optarg,"%d,%lf,%lf",&nphi,&dphi,&phi0);
+    break;
+  case 't' :
+    sscanf(optarg,"%d,%lf,%lf",&ntheta,&dtheta,&theta0);
+    break;
+  case 'h' :
+  default:
+    usage();
+    return -1;
+  }
+ }
+
+ if(nphi<=0 || dphi<=0. || ntheta<=0 || dtheta<=0.) {
+   cout<<"bad value for transverse dim."<<endl;
+   return -2;
+ }
+ cout<<"nphi="<<nphi<<" dphi="<<dphi<<" phi0="<<phi0<<endl;
+ cout<<"ntheta="<<ntheta<<" dtheta="<<dtheta<<" theta0="<<theta0<<" deg"<<endl;
+
+ vector<string> sphname;
+ for(int i=optind;i<narg;i++) sphname.push_back(arg[i]);
+ if(sphname.size()==0) {usage(); return -3;}
+ cout<<"Number of spheres "<<sphname.size()<<endl;
+
+ TArray<r_4> Cube(sphname.size(),nphi,ntheta);
+ Cube = 0.;
+ for(unsigned int is=0;is<sphname.size();is++) {
+
+   cout<<"... "<<sphname[is]<<endl;
+   PInPersist pis(sphname[is]);
+   SphereHEALPix<r_4> sph;
+   pis>>PPFNameTag("sph")>>sph;
+
+   for(int it=0;it<ntheta;it++) {
+     double t = theta0 + it*dtheta; t *= M_PI/180.;
+     if(t<0. || t>M_PI) continue;
+   for(int ip=0;ip<nphi;ip++) {
+     double p = phi0 + ip*dphi; p *= M_PI/180.;
+     while(p<0.) p += 2.*M_PI;
+     while(p>=2.*M_PI) p -= 2.*M_PI;
+     int_4 ksph = sph.PixIndexSph(t,p);
+     Cube(is,ip,it) = sph(ksph);
+   }
+   }
+
+ }
+ 
+ POutPersist pos(ppfrn);
+ pos<<PPFNameTag("cube")<<Cube;
+
+ return 0;
+}
+
+/*
+openppf cmvgsm2cube.ppf
+
+print cube
+
+objaoper cube sliceyz 0
+objaoper cube sliceyz 128
+objaoper cube sliceyz 255
+objaoper cube slicexy 90
+
+ */
Index: /trunk/Cosmo/SimLSS/cmvgsm2sph.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvgsm2sph.cc	(revision 3790)
+++ /trunk/Cosmo/SimLSS/cmvgsm2sph.cc	(revision 3790)
@@ -0,0 +1,73 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <fstream>
+#include <iostream>
+#include <string.h>
+#include <unistd.h>
+
+#include "skymap.h"
+
+
+int main(int narg,char *arg[])
+{
+ int nside = 512;
+ string ppfrn = "cmvgsm2sph";
+
+ if(narg<=1) {
+   cout<<"cmvgsm2sph liste_gsm.list [ppfroot_name]"<<endl;
+   return -1;
+ }
+
+ ifstream flist(arg[1], ofstream::in);
+ if(!flist)
+   {cout<<"ERROR: open failed for "<<arg[1]<<endl; return -2;}
+ if(narg>2) ppfrn = arg[2];
+
+ int nfilrd = 0, nfrderr = 0;;
+ while(1) {
+   string fn; double f=-1.; int is = -1;
+   flist>>is>>f>>fn;
+   if(flist.eof()) break;
+   cout<<nfilrd<<" : is="<<is<<" f="<<f<<" fn="<<fn<<endl;
+
+   SphereHEALPix<r_4> sph(nside);
+   int npix = sph.NbPixels();
+   if(nfilrd==0) cout<<"Healpix sph nside="<<nside<<" npix="<<npix<<endl;
+
+   ifstream fin(fn.c_str(), ofstream::in);
+   if(!fin) {
+     cout<<"ERROR: open failed for "<<fn<<endl;
+     nfrderr++;
+     continue;
+   }
+
+   int nread = 0;
+   while(1) {
+     double v = 0;
+     fin>>v;
+     if(fin.eof()) break;
+     if(nread<npix) sph(nread) = v;
+     nread++;
+   }
+
+   if(nread!=npix) {
+     cout<<"ERROR: mismatch nread="<<nread<<" / "<<npix<<endl;
+     nfrderr++;
+   } else {
+     char str[1024];
+     sprintf(str,"%s_%d.ppf",ppfrn.c_str(),is);
+     POutPersist pos(str);
+     sph.Info()["ISeq"] = (int_4)is;
+     sph.Info()["FMHz"] = (double)f;
+     sph.Info()["FName"] = fn;
+     pos<<PPFNameTag("sph")<<sph;
+     nfilrd++;
+   }
+ }
+ cout<<"\nread "<<nfilrd<<" files, with "<<nfrderr<<" echecs"<<endl;
+
+ return 0;
+}
+
Index: /trunk/Cosmo/SimLSS/cmvgsmcfit.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvgsmcfit.cc	(revision 3790)
+++ /trunk/Cosmo/SimLSS/cmvgsmcfit.cc	(revision 3790)
@@ -0,0 +1,183 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "fioarr.h"
+#include "ntuple.h"
+
+#include "Minuit2/FunctionMinimum.h"
+#include "Minuit2/MnMigrad.h"
+#include "Minuit2/MnUserParameters.h"
+#include "Minuit2/MnPrint.h"
+#include "Minuit2/MnMinos.h"
+#include "Minuit2/MinosError.h"
+#include "Minuit2/FCNBase.h"
+using namespace ROOT::Minuit2;
+
+//--------------------------------------------------
+namespace ROOT { namespace Minuit2 {
+class MyFCN : public FCNBase {
+public:
+  MyFCN(double nu0,double dnu,TVector<r_8>& V,TVector<r_8>& EV);
+  void InitFit(MnUserParameters& upar);
+  double Func(int ix,const vector<double>& par) const;
+  double operator()(const vector<double>& par) const;
+  double Up(void) const {return 1.;}
+  //protected:
+  mutable int _ndata;
+  double _nu0, _dnu;
+  TVector<r_8> _V;
+  TVector<r_8> _EV;
+};
+} }  // namespace ROOT + Minuit2
+
+//--------------------------------------------------
+int main(int narg,char *arg[])
+{
+  double nu0 = 820.0, dnu = 0.5; // MHz
+  if(narg<=1) {
+    cout<<"cmvgsmcfit cubefile.ppf"<<endl;
+    return -1;
+  }
+
+  POutPersist pos("cmvgsmcfit.ppf");
+  PInPersist pis(arg[1]);
+  TArray<r_4> Cube;
+  pis>>PPFNameTag("cube")>>Cube;
+  TVector<r_8> V(Cube.SizeX()), EV(Cube.SizeX());
+
+  const int nvar = 6;
+  const char *vname[nvar] = {"x2r","ok","ph","th","A","N"};
+  NTuple nt(nvar,vname);
+  double xnt[nvar]; for(int i=0;i<nvar;i++) xnt[i] = 0.;
+
+  int nfit = 0;
+  for(int j=0;j<Cube.SizeY();j++) {
+    cout<<"..."<<j<<endl;
+  for(int k=0;k<Cube.SizeZ();k++) {
+    for(int i=0;i<Cube.SizeX();i++) {
+      V(i) = Cube(i,j,k);
+      EV(i) = V(i) * 0.05;
+    }
+
+    MyFCN myfcn(nu0,dnu,V,EV);
+    MnUserParameters upar;
+    myfcn.InitFit(upar);
+    if(nfit==0) cout<<upar;
+
+    unsigned int nfcall = 10000;
+    double tolerance = 0.0001;
+    MnMigrad migrad(myfcn,upar,MnStrategy(1));
+    FunctionMinimum fmini = migrad(nfcall,tolerance);
+
+    MnUserParameterState ustate = fmini.UserState();
+    bool okfit = fmini.IsValid();
+    double xi2r = ustate.Fval()/(V.Size()-ustate.VariableParameters());
+    vector<double> par = ustate.Params();
+    if(nfit==0) {
+      cout<<ustate;
+      // cout<<ustate.Parameters(); cout<<ustate.Covariance(); cout<<ustate.GlobalCC();
+      cout<<"nFCN = "<<ustate.NFcn()<<", FVal = "<<ustate.Fval()<<", Edm = "<<ustate.Edm()<<endl;
+      cout<<"xi2r = "<<xi2r<<" ndata="<<myfcn._ndata<<" nvarp="<<ustate.VariableParameters()<<endl;
+    }
+
+    xnt[0] = xi2r; xnt[1] = (okfit)? 1.: 0.;
+    xnt[2] = j; xnt[3] = k;
+    for(unsigned int i=0;i<par.size();i++) xnt[4+i] = par[i];
+    nt.Fill(xnt);
+    //vector<double> err = ustate.Errors();
+    for(int i=0;i<Cube.SizeX();i++) Cube(i,j,k) -= myfcn.Func(i,par);
+
+    if(nfit<10) {
+      TVector<r_8> F(Cube.SizeX()), R(Cube.SizeX());
+      for(int i=0;i<Cube.SizeX();i++) {
+        F(i) = myfcn.Func(i,par);
+        R(i) = V(i) - F(i);
+      }
+      char str[64];
+      sprintf(str,"v_%d_%d",j,k);
+      pos<<PPFNameTag(str)<<V;
+      sprintf(str,"f_%d_%d",j,k);
+      pos<<PPFNameTag(str)<<F;
+      sprintf(str,"r_%d_%d",j,k);
+      pos<<PPFNameTag(str)<<R;
+    }
+    nfit++;
+    if(nfit>10) break;
+  }
+  }
+
+  pos<<PPFNameTag("nt")<<nt;
+  pos<<PPFNameTag("cublis")<<Cube;
+
+}
+
+//--------------------------------------------------
+MyFCN::MyFCN(double nu0, double dnu,TVector<r_8>& V,TVector<r_8>& EV)
+  : _V(V,false), _EV(EV,false)
+{
+  _ndata =_V.Size();
+  _nu0 = nu0;
+  _dnu = dnu;
+}
+
+void MyFCN::InitFit(MnUserParameters& upar)
+{
+  int n = _V.Size()-1;
+  double a = _V(0);
+  double ni = log(_V(n)/a) / log(1.+n*_dnu/_nu0);
+
+  upar.Add("A",a,a/50.);
+  upar.Add("N",ni,fabs(ni/50.));
+}
+
+double MyFCN::Func(int ix,const vector<double>& par) const
+// f = A*(nu/nu0)^N   avec  nu = nu0 + i*dnu
+{
+ double x = (_nu0 + ix*_dnu)/_nu0;
+ double f = par[0]*pow(x,par[1]);
+ return f;
+}
+
+double MyFCN::operator()(const vector<double>& par) const
+{
+ double xi2 = 0.;
+ _ndata = 0;
+ for(int i=0;i<_V.Size();i++) {
+   double e2 = _EV(i)*_EV(i);
+   if(e2<=0.) continue;
+   double v = _V(i);
+   double f = Func(i,par);
+   xi2 += (v-f)*(v-f)/e2;
+   _ndata++;
+ }
+ //cout<<"par: "<<par[0]<<" "<<par[1]<<endl;
+ return xi2;
+}
+
+/*
+openppf cmvgsmcfit.ppf
+
+set j 0
+set k 0
+n/plot v_${j}_${k}.val%n ! ! "nsta cpts"
+n/plot f_${j}_${k}.val%n ! ! "nsta cpts same red"
+n/plot r_${j}_${k}.val%n ! ! "nsta cpts"
+
+n/plot nt.ok%_nl ! ! "cpts"
+n/plot nt.x2r%_nl ! ! "cpts"
+
+n/plot nt.A%_nl ! ! "cpts"
+n/plot nt.N%_nl ! ! "cpts"
+
+set n 0
+objaoper cublis slicexz $n
+disp slicexz_$n
+
+n/plot cublis.val%x fabs(y-5)<0.1&&fabs(z-5)<0.1
+*/
Index: /trunk/Cosmo/SimLSS/cmvrvloscor.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvrvloscor.cc	(revision 3789)
+++ /trunk/Cosmo/SimLSS/cmvrvloscor.cc	(revision 3790)
@@ -20,5 +20,5 @@
 // set simul = 6_0p0_780
 // cmvrvloscorf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
 
 void usage(void);
@@ -33,4 +33,5 @@
 <<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
 <<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
+<<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
 <<endl;
 }
@@ -39,9 +40,9 @@
 {
  int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
- bool docube=true, dopk = false;
+ bool docube=true, dopk = false, do2d = false;
  
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"hn:K:SN")) != -1) {
+ while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
   switch (c) {
   case 'n' :
@@ -58,4 +59,7 @@
   case 'N' :
     docube = false;
+    break;
+  case '2' :
+    do2d = true;
     break;
   case 'h' :
@@ -93,5 +97,5 @@
  double nmax = max(Nx,max(Ny,Nz));
  cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
- Histo hmpc(-dmin*nmax,dmin*nmax,4.*nmax);
+ Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);
 
  POutPersist pos("cmvrvloscor.ppf");
@@ -107,4 +111,10 @@
    rgen = &(fluct3d->GetRealArray());
    *rgen = 0.;
+   cout<<"rgen:  size [1]="<<rgen->SizeX()
+                           <<" [2]="<<rgen->SizeY()
+                           <<" [3]="<<rgen->SizeZ()<<endl;
+   cout<<"pkgen: size [1]="<<fluct3d->GetComplexArray().SizeX()
+                           <<" [2]="<<fluct3d->GetComplexArray().SizeY()
+                           <<" [3]="<<fluct3d->GetComplexArray().SizeZ()<<endl;
  }
 
@@ -147,4 +157,9 @@
  for(int i=0;i<Nx;i++) {
    if(i%(Nx/10)==0) cout<<"   i="<<i<<endl;
+   TMatrix<r_4> M2d, M2dc;
+   if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
+     M2d.ReSize(Ny,Nz); M2d = 0.;
+     M2dc.ReSize(Ny,Nz); M2dc = 0.;
+   }
  for(int j=0;j<Ny;j+=nplany) {
    bool fhis = false;
@@ -161,10 +176,13 @@
      double lpd = (double)l + d/Dz; // valeur du deplacee
      // on repartit proportionellement au recouvrement sur 2 pixels
-     long l1 = long(lpd); // pixel de droite
-     long l2 = l1 + 1;  // pixel de gauche
+     long l1 = long(lpd); // pixel de gauche
+     long l2 = l1 + 1;  // pixel de droite
      lpd -= (double)l1;  // recouvrement du pixel du dessus
      if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
      if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
    }
+   // On remplit eventuellement les matrices 2D
+   if(do2d && M2d.Size()>0)
+     for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dc(j,l) = Rdis(l);}
    // On remplit le cube avec le champ R redshift distordu
    if(fluct3d) for(int l=0;l<Nz;l++) (*rgen)(l,j,i) += Rdis(l);
@@ -203,4 +221,11 @@
    nlosread++;
  }
+   if(do2d && M2d.Size()>0) {
+     char str[64];
+     sprintf(str,"mx_%d",i);
+     pos.PutObject(M2d,str);
+     sprintf(str,"mxc_%d",i);
+     pos.PutObject(M2dc,str);
+   }
  }
 
@@ -324,3 +349,9 @@
 # proj selon kT (black), selon kZ (red)
 n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
+
+# les matrices 2D
+set n 0
+disp mx_$n
+newwin
+disp mxc_$n
  */
Index: /trunk/Cosmo/SimLSS/cmvtuniv.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3789)
+++ /trunk/Cosmo/SimLSS/cmvtuniv.cc	(revision 3790)
@@ -10,4 +10,5 @@
 #include "ntuple.h"
 
+#include "constcosmo.h"
 #include "cosmocalc.h"
 #include "geneutils.h"
@@ -124,7 +125,12 @@
 {
  if(dz>z2-z1) dz = z2-z1;
+ double fz1 = 1000.*Fr_HyperFin_Par/(1.+z1);
+ double fz1pdz = 1000.*Fr_HyperFin_Par/(1.+z1+dz);
+ double fz2 = 1000.*Fr_HyperFin_Par/(1.+z2);
+ double fz2pdz = 1000.*Fr_HyperFin_Par/(1.+z2+dz);
  cout<<"\nTSTPRINT(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
- for(double z=z1;z<z2+dz/2.;z+=dz) {
-   cout<<"-- z = "<<z<<endl;
+ cout<<"  "<<fz1<<" -> "<<fz2<<" MHz,    df: "<<fz1-fz1pdz<<" , "<<fz2-fz2pdz<<" MHz"<<endl;
+ for(double z=z1;z<z2+dz/2.;z+=dz) {
+   cout<<"\n-- z = "<<z<<"      ("<<1000.*Fr_HyperFin_Par/(1.+z)<<" MHz)"<<endl;
    univ.Print(z);
    cout<<"Volume comoving in [z,z+"<<dz<<"] for 4Pi sr: "
@@ -138,5 +144,5 @@
    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;
+   cout<<"dLOScom/dz -> "<<dloscom<<" Mpc com"<<endl;
  }
 
Index: /trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3789)
+++ /trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3790)
@@ -435,7 +435,7 @@
  try {
    FitsImg3DWriter fwrt(cfname.c_str(),bitpix,5);
-   fwrt.WriteKey("NX",Nx_," axe transverse 1");
-   fwrt.WriteKey("NY",Ny_," axe transverse 2");
-   fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");
+   fwrt.WriteKey("NX",Nx_," axe transverse 1");  // NAXIS3
+   fwrt.WriteKey("NY",Ny_," axe transverse 2");  // NAXIS2
+   fwrt.WriteKey("NZ",Nz_," axe longitudinal (redshift)");  // NAXIS1
    fwrt.WriteKey("DX",Dx_," Mpc");
    fwrt.WriteKey("DY",Dy_," Mpc");
@@ -467,7 +467,7 @@
    FitsImg3DRead fimg(cfname.c_str(),0,5);
    fimg.Read(R_);
-   long nx = fimg.ReadKeyL("NX");
-   long ny = fimg.ReadKeyL("NY");
-   long nz = fimg.ReadKeyL("NZ");
+   long nx = fimg.ReadKeyL("NX");  // NAXIS3
+   long ny = fimg.ReadKeyL("NY");  // NAXIS2
+   long nz = fimg.ReadKeyL("NZ");  // NAXIS1
    double dx = fimg.ReadKey("DX");
    double dy = fimg.ReadKey("DY");
