Changeset 3781 in Sophya for trunk/Cosmo

Timestamp:

May 20, 2010, 11:19:52 AM (15 years ago)

Author:

cmv

Message:

ameliorations mineures, cmv 20/05/2010

Location:

trunk/Cosmo/SimLSS

Files:

• cmvginit3d.cc (modified) (1 diff)
• cmvrvloscor.cc (modified) (5 diffs)
• genefluct3d.cc (modified) (3 diffs)
• genefluct3d.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -507 +507 @@
    PrtTim(">>>> End Reading the Pk(vec(k)) cube");
 
+   if(compdspec&1) {
+     cout<<"\n--- Check realization spectra that has been re-read"<<endl;
+     HistoErr hpkgenf_rr(hpkgen); hpkgenf_rr.Zero();
+     fluct3d.ComputeSpectrum(hpkgenf_rr);
+     PrtTim(">>>> End Checking re-read realization 1D spectra");
+     posobs.PutObject(hpkgenf_rr,"hpkgenf_rr");
+   }
+
    //-----------------------------------------------------------------
    cout<<"\n--- Modifying cube for Transverse velocity"<<endl;

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -13 +13 @@
 #include "histos.h"
 #include "fabtcolread.h"
+#include "fftwserver.h"
 
 #include "constcosmo.h"
 #include "geneutils.h"
 #include "genefluct3d.h"
+// cmvrvloscorf -n 1,30 -K 75 -S ginit3d_6_0p0_100_r.fits ginit3d_6_0p0_100_rv.fits
 
 void usage(void);
 void usage(void)
 {
-  cout<<"cmvrvloscor rho.fits vlos.fits"<<endl;
+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 comsuming!!!)"<<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
+<<endl;
 }
 
 int main(int narg,char *arg[])
 {
- int nthread = 1;
- int_8 nhfill = 100000;
-
- if(narg<=2) {usage(); return -1;}
+ int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
+ bool docube=true, dopk = false;
+ 
+ // --- Decodage des arguments
+ char c;
+ while((c = getopt(narg,arg,"hn:K:SN")) != -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 'h' :
+  default :
+    usage(); return -1;
+  }
+ }
+ if(optind>=narg-1) {usage(); return -1;}
 
  //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
@@ -38 +71 @@
  InitTim();
 
- cout<<"> read rho: "<<arg[1]<<endl;
- FitsImg3DRead f3dr(arg[1],0,5);
- cout<<"> read vlos: "<<arg[2]<<endl;
- FitsImg3DRead f3dv(arg[2],0,5);
+ // --- 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");;
@@ -58 +92 @@
  cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
  Histo hmpc(-dmin*nmax,dmin*nmax,4.*nmax);
- nhfill = (int_8)Nx*(int_8)Ny*(int_8)Nz/nhfill; if(nhfill<=0) nhfill = 1;
 
  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<GEN3D_TYPE> R(Nz), V(Nz);
+ 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.;
+ }
+
+ // --- 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);
-
- cout<<"> filling redshift distorted cube"<<endl;
- int_8 nread = 0;
+ 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;
- for(int j=0;j<Ny;j++) {
+   if(i%(Nx/10)==0) cout<<"   i="<<i<<endl;
+ 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);
@@ -79 +153 @@
    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(nread%nhfill==0) hmpc.Add(d);
+     if(fhis) hmpc.Add(d);
      double lpd = (double)l + d/Dz; // valeur du deplacee
-     // on repartit proportionnelement au recouvrement sur 2 pixels
-     int l1 = int(lpd); // pixel de droite
-     int l2 = l1 + 1;  // pixel de gauche
+     // on repartit proportionellement au recouvrement sur 2 pixels
+     long l1 = long(lpd); // pixel de droite
+     long l2 = l1 + 1;  // pixel de gauche
      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;
    }
-   for(int l=0;l<Nz;l++) rgen(l,j,i) += 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++;
+ }
+ }
+
+ 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");
- pos.PutObject(hmpc,"hmpc");
-
- 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");
+
+ // --- 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
@@ -148 +291 @@
 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"
 
+# recosntructed 2D power spectrum
 imag hpkrec2
 addoval 0 0 0.05 0.05 "green" false

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -120 +120 @@
  NRtot_ = (int_8)Nx_*(int_8)Ny_*(int_8)Nz_; // nombre de pixels dans le survey
  NCz_ =  Nz_/2 +1;
+ NFcoef_ = (int_8)Nx_*(int_8)Ny_*(int_8)NCz_; // nombre de coeff de Fourier dans le survey
  NTz_ = 2*NCz_;
 
@@ -659 +660 @@
  cout<<"      Resol: dx="<<Dx_<<" dy="<<Dy_<<" dz="<<Dz_<<" Mpc"
      <<", dVol="<<dVol_<<", Vol="<<Vol_<<" Mpc^3"<<endl;
- cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<endl;
+ cout<<"Fourier Size : nx="<<Nx_<<" ny="<<Ny_<<" nz="<<NCz_<<" ncoeff="<<NFcoef_<<endl;
  cout<<"        Resol: dkx="<<Dkx_<<" dky="<<Dky_<<" dkz="<<Dkz_<<" Mpc^-1"
      <<", Dk3="<<Dk3_<<" Mpc^-3"<<endl;
@@ -840 +841 @@
  if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
 
- if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(NRtot_-nconj1-nconj2)-nreal<<endl;
+ if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(NFcoef_-nconj1-nconj2)-nreal<<endl;
 
  return nreal+nconj1+nconj2;

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -77 +77 @@
   vector<long> GetNpix(void) {return N_;}
   int_8 NPix(void) {return NRtot_;}
+  int_8 NCoeff(void) {return NFcoef_;}
   long Nx(void) {return Nx_;}
   long Ny(void) {return Ny_;}
@@ -176 +177 @@
   long Nx_,Ny_,Nz_;  vector<long> N_;
   long NCz_,NTz_;
-  int_8 NRtot_;
+  int_8 NRtot_,NFcoef_;
 
   double Dx_,Dy_,Dz_;  vector<double> D_;