#include "sopnamsp.h" #include "machdefs.h" #include #include #include #include #include #include #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 -o rvlos_${simul}.ppf -n 1,1 -z -1,0 -H 50 -2 -S 10000 -s 1. ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits void usage(void); void usage(void) { cout <<"cmvrvloscor [options] rho.fits vlos.fits"<=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: "< read vlos: "<Nz || nplanz<=0) nplanz = Nz; if(iplanz0<0)iplanz0 = 0; if(iplanz0>=Nz) {cout<<"Error: iplanz0="<Nz) nplanz = Nz - iplanz0; cout<<" process "< V4read(Nz); TVector R(nplanz), V(nplanz); TVector Rdis(nplanz); // --- Histo des deplacements Histo* hmpc = NULL; if(nhfilllos>0) { cout<<"...Fill Mpc displacement histo with "<0.) dmin *= fabs(scaldis); long nmax = max(Nx,max(Ny,Nz)); hmpc = new Histo(-dmin*nmax/10.,dmin*nmax/10.,nmax); cout<<" dmin="< fill one los every "< > FR, FV; TVector< complex > pkvr, FRdis; TVector pkr, pkrc; FFTWServer fftserv; fftserv.setNormalize(false); // pour accord avec GeneFluct3D if(npkmax>0) { cout<<"...compute R*conj(R) V*conj(V) V*conj(R) cross-power spectrum with "< fill one los every "<* rgen = NULL; if(docube) { cout<<"...Create and fill 3D cube"<Print(); rgen = &(fluct3d->GetRealArray()); *rgen = 0.; cout<<"rgen: size [1]="<SizeX() <<" [2]="<SizeY() <<" [3]="<SizeZ()<>> filling redshift distorted cube, scale displacement by "< M2d, M2dv, M2dc; if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) { M2d.ReSize(Ny,nplanz); M2d = 0.; M2dv.ReSize(Ny,nplanz); M2dv = 0.; M2dc.ReSize(Ny,nplanz); M2dc = 0.; } for(int j=0;jAdd(d); double lpd = (double)l + d/Dz; // valeur du deplacee if(dodistrib) { // 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=0 && l2=0 && l10) for(int l=0;l0 && nlosread%npkmax==0) { fftserv.FFTForward(V,FV); int nf = FV.Size(); if(pkvr.Size()<=0) { cout<<"...Create vector for cross-power spectrum computation sz="<(0.); } fftserv.FFTForward(R,FR); fftserv.FFTForward(Rdis,FRdis); for(int l=0;l0) { char str[64]; sprintf(str,"mx_%d",i); pos.PutObject(M2d,str); sprintf(str,"mxv_%d",i); pos.PutObject(M2dv,str); sprintf(str,"mxc_%d",i); pos.PutObject(M2dc,str); } } // fin boucle sur nx cout<<"Number of processed los: "<Show(); if(hmpc->NEntries()>0) pos.PutObject(*hmpc,"hmpc"); } if(npkav>0) { double nn = Dx*Dy*Dz/double(nplanz)/double(npkav); // moyenne + normalisation / hpkrec cout<<"Number of averaged spectra: npkav="<>>> 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"<ReComputeFourier(); PrtTim(">>>> End ReComputing spectrum"); // Compute 1D spectrum cout<GetKincMin(), knyqmax = fluct3d->GetKmax(); long nherr = long(knyqmax/dkmin+0.5); cout<<"\nFor HistoErr: d="<>>> End Computing final spectrum"); // Compute 2D spectrum cout<<"\n--- Computing final 2D spectrum"<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="<ComputeSpectrum2D(hpkrec2); pos.PutObject(hpkrec2,"hpkrec2"); PrtTim(">>>> End Computing final 2D spectrum"); } // --- end of job, write objects in ppf dvllos("nlosread") = (int_4)nlosread; dvllos("nlospred") = (int_4)nlospred; dvllos("Nx") = (int_4)Nx; dvllos("Ny") = (int_4)Ny; dvllos("Nz") = (int_4)Nz; dvllos("Dx") = (r_8)Dx; dvllos("Dy") = (r_8)Dy; dvllos("Dz") = (r_8)Dz; dvllos("nplanz") = (int_4)nplanz; dvllos("iplanz0") = (int_4)iplanz0; dvllos("nplany") = (int_4)nplany; dvllos("nplanx") = (int_4)nplanx; dvllos("npkav") = (int_4)npkav; dvllos("dkzpk") = (r_8)dkzpk; pos.PutObject(dvllos,"dvllos"); if(fluct3d != NULL) delete fluct3d; if(hmpc != NULL) delete hmpc; //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH } catch (PException& exc) { cerr<<"cmvrvloscor.cc catched PException"<0&&val>0 ! "nsta cpts logx logy" n/plot pkrc.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy same red" n/plot pkvr.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy" # reconstructed 1D power spectrum n/plot hpkrec.val%x x>0&&val>0 ! "nsta cpts logx logy" # 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 defscript tom del m2 c++exec TMatrix m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \ for(int i=0;i