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

Last change on this file since 3989 was 3821, checked in by cmv, 15 years ago

type d'interpolation pour les spectres tabules, amelioration de cmvrvloscor.cc, cmv 30/07/2010
File size: 12.6 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <iostream>
4#include <stdlib.h>
5#include <stdio.h>
6#include <string.h>
7#include <math.h>
8#include <unistd.h>
9
10#include "sophyainit.h"
11#include "timing.h"
12#include "dvlist.h"
13#include "histos.h"
14#include "fabtcolread.h"
15#include "fftwserver.h"
16
17#include "constcosmo.h"
18#include "geneutils.h"
19#include "genefluct3d.h"
20// set simul = 6_0p0_780
21// 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
22
23void usage(void);
24void usage(void)
25{
26cout
27<<"cmvrvloscor [options] rho.fits vlos.fits"<<endl
28<<"-n nplany,nplanx: process one Y (X) plane every \"nplany\" (\"nplanx\") (def:1,1(all))"<<endl
29<<"-z nplanz,iplanz0: process only \"nplanz\" Z plane beginnig at plane \"iplanz0\" (def:all Z planes)"<<endl
30<<"-H nhfill: fill displacement histo with \"nhfill\" los (def:25)"<<endl
31<<"-s scaldis: scale the displacement by factor \"scaldis\" (def: 1.)"<<endl
32<<"-p: do not distribute dRho/Rho in pixel, just take the nearest (def: no, distribute rho)"<<endl
33<<"-S npkmax: compute cross-power spectrum of V*conj(R) with \"npkmax\" los (def: no)"<<endl
34<<"-2: compute 2D projection for dRho and dRho(corrected) for 3 Y-Z plans (def=no)"<<endl
35<<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
36<<"-o ppfout.ppf: output ppf file name (def=\"cmvrvloscor.ppf\")"<<endl
37<<endl;
38}
39
40int main(int narg,char *arg[])
41{
42 int nthread = 1, nplany=1, nplanx=1, nplanz=-1, iplanz0=0, nhfilllos = 25, npkmax = 0;
43 double scaldis = 1.;
44 bool docube=true, do2d = false, dodistrib=true;
45 string fnppf = "cmvrvloscor.ppf";
46
47 // --- Decodage des arguments
48 char c;
49 while((c = getopt(narg,arg,"hn:H:S:Np2o:s:z:")) != -1) {
50 switch (c) {
51 case 'n' :
52 sscanf(optarg,"%d,%d",&nplany,&nplanx);
53 if(nplany<=0) nplany = 1;
54 if(nplanx<=0) nplanx = 1;
55 break;
56 case 'z' :
57 sscanf(optarg,"%d,%d",&nplanz,&iplanz0);
58 break;
59 case 'H' :
60 nhfilllos = atoi(optarg);
61 if(nhfilllos<=0) nhfilllos = 0;
62 break;
63 case 'S' :
64 npkmax = atoi(optarg);
65 break;
66 case 'N' :
67 docube = false;
68 break;
69 case 'p' :
70 dodistrib = false;
71 break;
72 case '2' :
73 do2d = true;
74 break;
75 case 'o' :
76 fnppf = optarg;
77 break;
78 case 's' :
79 scaldis = atof(optarg);
80 break;
81 case 'h' :
82 default :
83 usage(); return -1;
84 }
85 }
86 if(optind>=narg-1) {usage(); return -1;}
87
88 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
89 try {
90 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
91
92 SophyaInit();
93 InitTim();
94
95 // --- open FITS files (dRho/Rho and Vlos)
96 cout<<"> read rho: "<<arg[optind]<<endl;
97 FitsImg3DRead f3dr(arg[optind],0,5);
98 cout<<"> read vlos: "<<arg[optind+1]<<endl;
99 FitsImg3DRead f3dv(arg[optind+1],0,5);
100 long Nx = f3dr.ReadKeyL("Nx");
101 long Ny = f3dr.ReadKeyL("Ny");;
102 long Nz = f3dr.ReadKeyL("Nz");;
103 cout<<"N: x="<<Nx<<" y="<<Ny<<" z="<<Nz<<endl;
104 double Dx = f3dr.ReadKey("Dx");
105 double Dy = f3dr.ReadKey("Dy");
106 double Dz = f3dr.ReadKey("Dz");
107 cout<<"D: x="<<Dx<<" y="<<Dy<<" z="<<Dz<<endl;
108 double Zref = f3dr.ReadKey("ZREF");
109 double Href = f3dr.ReadKey("HREF");
110 cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
111
112 POutPersist pos(fnppf.c_str());
113 DVList dvllos;
114
115 // --- Read and process data
116 if(nplany>Ny || nplany<=0) nplany = Ny;
117 if(nplanx>Nx || nplanx<=0) nplanx = Nx;
118 cout<<"...Will read one Y plane every "<<nplany<<endl;
119 cout<<" one X plane every "<<nplanx<<endl;
120 int nlospred = Nx*Ny/nplany/nplanx;
121 cout<<" predicted number of los to be read "<<nlospred<<endl;
122 if(nplanz>Nz || nplanz<=0) nplanz = Nz;
123 if(iplanz0<0)iplanz0 = 0;
124 if(iplanz0>=Nz) {cout<<"Error: iplanz0="<<iplanz0<<" out of range"<<endl; return -2;}
125 if(iplanz0+nplanz>Nz) nplanz = Nz - iplanz0;
126 cout<<" process "<<nplanz<<" Z planes in ["<<iplanz0<<","<<iplanz0+nplanz<<"["<<endl;
127 TVector<r_4> V4read(Nz);
128 TVector<r_4> R(nplanz), V(nplanz);
129 TVector<r_8> Rdis(nplanz);
130
131 // --- Histo des deplacements
132 Histo* hmpc = NULL;
133 if(nhfilllos>0) {
134 cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
135 double dmin = min(Dx,min(Dy,Dz)); if(fabs(scaldis)>0.) dmin *= fabs(scaldis);
136 long nmax = max(Nx,max(Ny,Nz));
137 hmpc = new Histo(-dmin*nmax/10.,dmin*nmax/10.,nmax);
138 cout<<" dmin="<<dmin<<" nmax="<<nmax<<endl;
139 nhfilllos = int((double)nlospred/nhfilllos + 0.5);
140 if(nhfilllos<=0) nhfilllos = 1;
141 cout<<" -> fill one los every "<<nhfilllos<<endl;
142 }
143
144 // --- Vector for PK cross-correlation computation
145 if(npkmax<0) npkmax = Nx*Ny;
146 TVector< complex<r_4> > FR, FV;
147 TVector< complex<r_8> > pkvr, FRdis;
148 TVector<r_8> pkr, pkrc;
149 FFTWServer fftserv;
150 fftserv.setNormalize(false); // pour accord avec GeneFluct3D
151 if(npkmax>0) {
152 cout<<"...compute R*conj(R) V*conj(V) V*conj(R) cross-power spectrum with "<<npkmax<<" los"<<endl;
153 npkmax = int((double)nlospred/npkmax + 0.5);
154 if(npkmax<=0) npkmax = 1;
155 cout<<" -> fill one los every "<<npkmax<<endl;
156 }
157 int npkav = 0;
158 double dkzpk = 2.*M_PI/(nplanz*Dz); // Mpc^-1
159
160 // --- Create a Cube for analysis
161 GeneFluct3D *fluct3d = NULL;
162 TArray<GEN3D_TYPE>* rgen = NULL;
163 if(docube) {
164 cout<<"...Create and fill 3D cube"<<endl;
165 fluct3d = new GeneFluct3D(Nx,Ny,nplanz,Dx,Dy,Dz,nthread,2);
166 fluct3d->Print();
167 rgen = &(fluct3d->GetRealArray());
168 *rgen = 0.;
169 cout<<"rgen: size [1]="<<rgen->SizeX()
170 <<" [2]="<<rgen->SizeY()
171 <<" [3]="<<rgen->SizeZ()<<endl;
172 cout<<"pkgen: size [1]="<<fluct3d->GetComplexArray().SizeX()
173 <<" [2]="<<fluct3d->GetComplexArray().SizeY()
174 <<" [3]="<<fluct3d->GetComplexArray().SizeZ()<<endl;
175 }
176
177 // --- Do the jobs !!!
178 cout<<">>> filling redshift distorted cube, scale displacement by "<<scaldis<<endl;
179 int nlosread = 0;
180 int lpri = nlospred/25; if(lpri==0) lpri = 1;
181 for(int i=0;i<Nx;i+=nplanx) {
182 TMatrix<r_4> M2d, M2dv, M2dc;
183 if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
184 M2d.ReSize(Ny,nplanz); M2d = 0.;
185 M2dv.ReSize(Ny,nplanz); M2dv = 0.;
186 M2dc.ReSize(Ny,nplanz); M2dc = 0.;
187 }
188 for(int j=0;j<Ny;j+=nplany) {
189 if(nlosread%lpri==0) cout<<" nlosread="<<nlosread<<" (i="<<i<<",j="<<j<<")"<<endl;
190 bool fhis = (hmpc != NULL && nlosread%nhfilllos==0) ? true: false;
191 //for(int l=0;l<nplanz;l++) R(l) = f3dr.Read(iplanz0+l,j,i);
192 //for(int l=0;l<nplanz;l++) V(l) = f3dv.Read(iplanz0+l,j,i);
193 f3dr.Read(j,i,V4read);
194 for(int l=0;l<nplanz;l++) R(l) = V4read(iplanz0+l);
195 f3dv.Read(j,i,V4read);
196 for(int l=0;l<nplanz;l++) V(l) = V4read(iplanz0+l);
197 Rdis = 0.;
198 // Calcul du champ R redshift distordu
199 for(int l=0;l<nplanz;l++) {
200 double d = scaldis * V(l) / Href;
201 if(fhis) hmpc->Add(d);
202 double lpd = (double)l + d/Dz; // valeur du deplacee
203 if(dodistrib) {
204 // on repartit proportionellement au recouvrement sur 2 pixels
205 long l1 = long(lpd); // pixel de gauche
206 long l2 = l1 + 1; // pixel de droite
207 lpd -= (double)l1; // recouvrement du pixel du dessus
208 if(l1>=0 && l1<nplanz) Rdis(l1) += R(l) * (1.-lpd);
209 if(l2>=0 && l2<nplanz) Rdis(l2) += R(l) * lpd;
210 } else { // take nearest pixel
211 long l1 = long(lpd + 0.5);
212 if(l1>=0 && l1<nplanz) Rdis(l1) += R(l);
213 }
214 }
215 // On remplit eventuellement les matrices 2D
216 if(do2d && M2d.Size()>0)
217 for(int l=0;l<nplanz;l++) {M2d(j,l) = R(l); M2dv(j,l) = V(l); M2dc(j,l) = Rdis(l);}
218 // On remplit le cube avec le champ R redshift distordu
219 if(fluct3d) for(int l=0;l<nplanz;l++) (*rgen)(l,j,i) += Rdis(l);
220 // Cross-power spectrum computation
221 if(npkmax>0 && nlosread%npkmax==0) {
222 fftserv.FFTForward(V,FV);
223 int nf = FV.Size();
224 if(pkvr.Size()<=0) {
225 cout<<"...Create vector for cross-power spectrum computation sz="<<nf<<endl;
226 pkr.ReSize(nf); pkr = 0.;
227 pkrc.ReSize(nf); pkrc = 0.;
228 pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
229 }
230 fftserv.FFTForward(R,FR);
231 fftserv.FFTForward(Rdis,FRdis);
232 for(int l=0;l<nf;l++) {
233 pkr(l) += norm(FR(l));
234 pkrc(l) += norm(FRdis(l));
235 pkvr(l) += FV(l)*conj(FR(l));
236 }
237 npkav++;
238 }
239 nlosread++;
240 } // fin boucle sur ny
241 if(do2d && M2d.Size()>0) {
242 char str[64];
243 sprintf(str,"mx_%d",i);
244 pos.PutObject(M2d,str);
245 sprintf(str,"mxv_%d",i);
246 pos.PutObject(M2dv,str);
247 sprintf(str,"mxc_%d",i);
248 pos.PutObject(M2dc,str);
249 }
250 } // fin boucle sur nx
251
252 cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
253 if(hmpc != NULL) {
254 hmpc->Show();
255 if(hmpc->NEntries()>0) pos.PutObject(*hmpc,"hmpc");
256 }
257 if(npkav>0) {
258 double nn = Dx*Dy*Dz/double(nplanz)/double(npkav); // moyenne + normalisation / hpkrec
259 cout<<"Number of averaged spectra: npkav="<<npkav<<" norm="<<nn<<endl;
260 pkr *= nn;
261 pkrc *= nn;
262 pkvr *= nn;
263 pos.PutObject(pkr,"pkr");
264 pos.PutObject(pkrc,"pkrc");
265 pos.PutObject(pkvr,"pkvr");
266 }
267 PrtTim(">>>> End filling redshift distorted cube");
268
269 // --- Fourier transform 3D cube and compute 1D and 2D power spectra
270 if(fluct3d) {
271 cout<<">>> Fourier transform 3D cube and compute 1D and 2D power spectra"<<endl;
272 // do the FFT for spectrum analysis
273 fluct3d->ReComputeFourier();
274 PrtTim(">>>> End ReComputing spectrum");
275
276 // Compute 1D spectrum
277 cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
278 double dkmin = fluct3d->GetKincMin(), knyqmax = fluct3d->GetKmax();
279 long nherr = long(knyqmax/dkmin+0.5);
280 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
281 HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
282 hpkrec.ReCenterBin(); hpkrec.Show();
283 fluct3d->ComputeSpectrum(hpkrec);
284 pos.PutObject(hpkrec,"hpkrec");
285 PrtTim(">>>> End Computing final spectrum");
286
287 // Compute 2D spectrum
288 cout<<"\n--- Computing final 2D spectrum"<<endl;
289 double dktmin = fluct3d->GetKTincMin(), ktnyqmax = fluct3d->GetKTmax();
290 double dkzmin = fluct3d->GetKinc()[2], kznyqmax = fluct3d->GetKnyq()[2];
291 long nherrt = long(ktnyqmax/dktmin+0.5), nherrz = long(kznyqmax/dkzmin+0.5);
292 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
293 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
294 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
295 hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
296 fluct3d->ComputeSpectrum2D(hpkrec2);
297 pos.PutObject(hpkrec2,"hpkrec2");
298 PrtTim(">>>> End Computing final 2D spectrum");
299 }
300
301 // --- end of job, write objects in ppf
302 dvllos("nlosread") = (int_4)nlosread;
303 dvllos("nlospred") = (int_4)nlospred;
304 dvllos("Nx") = (int_4)Nx;
305 dvllos("Ny") = (int_4)Ny;
306 dvllos("Nz") = (int_4)Nz;
307 dvllos("Dx") = (r_8)Dx;
308 dvllos("Dy") = (r_8)Dy;
309 dvllos("Dz") = (r_8)Dz;
310 dvllos("nplanz") = (int_4)nplanz;
311 dvllos("iplanz0") = (int_4)iplanz0;
312 dvllos("nplany") = (int_4)nplany;
313 dvllos("nplanx") = (int_4)nplanx;
314 dvllos("npkav") = (int_4)npkav;
315 dvllos("dkzpk") = (r_8)dkzpk;
316 pos.PutObject(dvllos,"dvllos");
317 if(fluct3d != NULL) delete fluct3d;
318 if(hmpc != NULL) delete hmpc;
319
320 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
321 } catch (PException& exc) {
322 cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
323 return 77;
324 } catch (std::exception& sex) {
325 cerr << "cmvrvloscor.cc std::exception :"
326 << (string)typeid(sex).name() << "\n msg= "
327 << sex.what() << endl;
328 return 78;
329 } catch (...) {
330 cerr << "cmvrvloscor.cc catched unknown (...) exception " << endl;
331 return 79;
332 }
333 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
334
335 return 0;
336}
337
338/*
339openppf cmvrvloscor.ppf
340
341c++exec dvllos.Print();
342set dkz ${dvllos.dkzpk}
343
344disp hmpc
345
346# cross-power spectrum
347n/plot pkr.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy"
348n/plot pkrc.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy same red"
349
350n/plot pkvr.val%n*${dkz} n>0&&val>0 ! "nsta cpts logx logy"
351
352# reconstructed 1D power spectrum
353n/plot hpkrec.val%x x>0&&val>0 ! "nsta cpts logx logy"
354
355# reconstructed 2D power spectrum
356imag hpkrec2
357addoval 0 0 0.05 0.05 "green" false
358addoval 0 0 0.1 0.1 "green" false
359addoval 0 0 0.25 0.25 "green" false
360addoval 0 0 0.5 0.5 "green" false
361x = ${hpkrec2.xmax} / 2.
362addoval 0 0 $x $x "green" false
363x = ${hpkrec2.ymax} / 2.
364addoval 0 0 $x $x "green" false
365
366# proj selon kT (black), selon kZ (red)
367n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
368
369# les matrices 2D
370set n 0
371disp mx_$n
372newwin
373disp mxc_$n
374
375defscript tom
376del m2
377c++exec TMatrix<r_8> m2(hpkrec2.NBinX(),hpkrec2.NBinY()); \
378for(int i=0;i<hpkrec2.NBinX();i++) \
379for(int j=0;j<hpkrec2.NBinY();j++) m2(i,j) = hpkrec2(i,j); \
380KeepObj(m2); cout<<"OK"<<endl;
381endscript
382
383###### Comparaison avec le spectre initial
384cd /home
385mkdir ok
386mkdir gene
387cd /home
388
389cd /ok
390openppf cmvrvloscor.ppf
391tom
392cd /home
393
394cd /gene
395openppf ginit3d_?_?p?_???.ppf
396tom
397cd /home
398
399h/oper / /ok/m2 /gene/m2 mdiff
400disp mdiff
401 */
Note: See TracBrowser for help on using the repository browser.