source: Sophya/trunk/Cosmo/SimLSS/cmvrvloscor.cc@ 3796

#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 "histos.h"
#include "fabtcolread.h"
#include "fftwserver.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "genefluct3d.h"
// set simul = 6_0p0_780
// cmvrvloscorf -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
// cmvrvloscorf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits

void usage(void);
void usage(void)
{
cout
<<"cmvrvloscor [options] rho.fits vlos.fits"<<endl
<<"-n nplany,nhfill: process one Y plane every \"nplany\" (def:1(all))"<<endl
<<"                  fill histos with \"nhfill\" los (def:25)"<<endl
<<"-K npix: compute correlation R*V at +/- npix pixels (def: no)"<<endl
<<"         (very time consuming!!!)"<<endl
<<"-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;
}

int main(int narg,char *arg[])
{
 int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
 bool docube=true, dopk = false, do2d = false;
 
 // --- Decodage des arguments
 char c;
 while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
  switch (c) {
  case 'n' :
    sscanf(optarg,"%d,%d",&nplany,&nhfilllos);
    if(nplany<=0) nplany = 1;
    if(nhfilllos<=0) nhfilllos = 0;
    break;
  case 'K' :
    npixcor = atoi(optarg);
    break;
  case 'S' :
    dopk = true;
    break;
  case 'N' :
    docube = false;
    break;
  case '2' :
    do2d = true;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }
 if(optind>=narg-1) {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

 SophyaInit();
 InitTim();

 // --- open FITS files (dRho/Rho and Vlos)
 cout<<"> read rho: "<<arg[optind]<<endl;
 FitsImg3DRead f3dr(arg[optind],0,5);
 cout<<"> read vlos: "<<arg[optind+1]<<endl;
 FitsImg3DRead f3dv(arg[optind+1],0,5);
 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;

 double dmin = min(Dx,min(Dy,Dz));
 long nmax = max(Nx,max(Ny,Nz));
 cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
 Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);

 POutPersist pos("cmvrvloscor.ppf");
 DVList dvlcor;

 // --- Create a Cube for analysis
 GeneFluct3D *fluct3d = NULL;
 TArray<GEN3D_TYPE>* rgen = NULL;
 if(docube) {
   cout<<"...Create and fill 3D cube"<<endl;
   fluct3d = new GeneFluct3D(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
   fluct3d->Print();
   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;
 }

 // --- Vector for real-space correlation computation
 int imil = Nz-1;
 dvlcor("imil") = (int_4)imil;
 TVector<int_4> nKsi;
 TVector<r_8> Ksirv, Ksirvc, Ksirr, Ksirrc;
 if(npixcor>0) {
   Ksirv.ReSize(2*Nz-1);  Ksirv  = 0.;
   Ksirvc.ReSize(2*Nz-1); Ksirvc = 0.;
   Ksirr.ReSize(2*Nz-1);  Ksirr  = 0.;
   Ksirrc.ReSize(2*Nz-1); Ksirrc = 0.;
   nKsi.ReSize(2*Nz-1);   nKsi   = 0;
   cout<<"...Compute R*V correlation on +/-"<<npixcor<<" px"<<endl;
 }

 // --- Vector for PK cross-correlation computation
 int npk = 0;
 TVector< complex<r_4> > FR, FV;
 TVector< complex<r_8> > pkvr, FRdis;
 TVector<r_8> pkr, pkrc;
 FFTWServer fftserv;
 if(dopk) cout<<"...compute V*conj(R) cross-power spectrum"<<endl;

 // --- Read and process data
 TVector<r_4> R(Nz), V(Nz);
 TVector<r_8> Rdis(Nz);
 if(nplany>Ny) nplany = Ny;
 cout<<"...Will read one Y plane every "<<nplany<<endl;
 if(nhfilllos) {
   cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
   nhfilllos = int((double)Nx*Ny/nplany/nhfilllos + 0.5);
   if(nhfilllos<=0) nhfilllos = 1;
   cout<<"   -> fill one los every "<<nhfilllos<<endl;
 }

 cout<<">>> filling redshift distorted cube"<<endl;
 int nlosread = 0;
 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;
   if(nhfilllos) if(nlosread%nhfilllos==0) fhis = true;
   //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
   //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
   f3dr.Read(j,i,R);
   f3dv.Read(j,i,V);
   Rdis = 0.;
   // Calcul du champ R redshift distordu
   for(int l=0;l<Nz;l++) {
     double d = (1.+Zref) / Href * V(l);
     if(fhis) hmpc.Add(d);
     double lpd = (double)l + d/Dz; // valeur du deplacee
     // on repartit proportionellement au recouvrement sur 2 pixels
     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);
   // Calcul eventuel de la fonction de correlation R*V
   if(npixcor>0) {
     for(long l1=0;l1<Nz;l1++) {
       for(long l2=max(0L,l1-npixcor);l2<min(Nz,l1+npixcor);l2++) {
         int lc = imil+(l2-l1);
         Ksirr(lc)  += R(l1)*R(l2);
         Ksirrc(lc) += Rdis(l1)*R(l2);
         Ksirv(lc)  += R(l1)*V(l2);
         Ksirvc(lc) += Rdis(l1)*V(l2);
         nKsi(lc)++;
       }
     }
   }
   // Cross-power spectrum computation
   if(dopk) {
     fftserv.FFTForward(V,FV);
     int nf = FV.Size();
     if(pkvr.Size()<=0) {
       cout<<"...Create vector for cross-power spectrum computation"<<endl;
       pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
       pkr.ReSize(nf);  pkr = 0.;
       pkrc.ReSize(nf); pkrc = 0.;
     }
     fftserv.FFTForward(R,FR);
     for(int l=0;l<nf;l++) {
       pkvr(l) += FV(l)*conj(FR(l));
       pkr(l)  += norm(FR(l));
     }
     fftserv.FFTForward(Rdis,FRdis);
     for(int l=0;l<nf;l++) pkrc(l) += norm(FRdis(l));
     npk++;
   }
   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);
   }
 }

 cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
 dvlcor("nlosread") = (int_4)nlosread;
 if(hmpc.NEntries()>0) {
   hmpc.Show();
   pos.PutObject(hmpc,"hmpc");
 }
 if(Ksirr.Size()>0) {
   for(int l=0;l<Ksirr.Size();l++) if(nKsi(l)>0) {
       Ksirr(l)  /= nKsi(l);
       Ksirrc(l) /= nKsi(l);
       Ksirv(l)  /= nKsi(l);
       Ksirvc(l) /= nKsi(l);
     }
   pos.PutObject(Ksirr,"ksirr");
   pos.PutObject(Ksirrc,"ksirrc");
   pos.PutObject(Ksirv,"ksirv");
   pos.PutObject(Ksirvc,"ksirvc");
   pos.PutObject(nKsi,"nksi");
 }
 if(npk>0) {
   pkvr /= (double)npk;
   pkr  /= (double)npk;
   pkrc /= (double)npk;
   pos.PutObject(pkvr,"pkvr");
   pos.PutObject(pkr,"pkr");
   pos.PutObject(pkrc,"pkrc");
 }
 PrtTim(">>>> End filling redshift distorted cube");

 // --- Fourier transform 3D cube and compute 1D and 2D power spectra
 if(fluct3d) {
   cout<<">>> Fourier transform 3D cube and compute 1D and 2D power spectra"<<endl;
   // do the FFT for spectrum analysis
   fluct3d->ReComputeFourier();
   PrtTim(">>>> End ReComputing spectrum");

   // Compute 1D spectrum
   cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
   double dkmin = fluct3d->GetKincMin(), 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");

   // Compute 2D spectrum
   cout<<"\n--- Computing final 2D spectrum"<<endl;
   double dktmin = fluct3d->GetKTincMin(), ktnyqmax = fluct3d->GetKTmax();
   double dkzmin = fluct3d->GetKinc()[2], kznyqmax = fluct3d->GetKnyq()[2];
   long nherrt = long(ktnyqmax/dktmin+0.5), 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");
 }

 // --- end of job, write objects in ppf
 pos.PutObject(dvlcor,"dvlcor");
 if(fluct3d) delete fluct3d;

 //----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

disp hmpc

# cross-correlation
disp nksi
set imil ${dvlcor.imil}
n/plot ksirv.val%n-${imil} ! ! "nsta cpts"
n/plot ksirvc.val%n-${imil} ! ! "nsta cpts same red"

n/plot ksirr.val%n-${imil} ! ! "nsta cpts"
n/plot ksirrc.val%n-${imil} ! ! "nsta cpts same red"

# cross-power spectrum
n/plot pkr.val%n ! ! "nsta cpts logx"
n/plot pkrc.val%n ! ! "nsta cpts logx same red"

n/plot pkvr.val%n ! ! "nsta cpts logx"

# reconstructed 1D power spectrum
n/plot hpkrec.val%x x>0 ! "nsta cpts logx"

# reconstructed 2D power spectrum
imag hpkrec2
addoval 0 0 0.05 0.05 "green" false
addoval 0 0 0.1 0.1 "green" false
addoval 0 0 0.25 0.25 "green" false
addoval 0 0 0.5 0.5 "green" false
x = ${hpkrec2.xmax} / 2.
addoval 0 0 $x $x "green" false
x = ${hpkrec2.ymax} / 2.
addoval 0 0 $x $x "green" false

# 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
 */