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

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

cmv, 29/06/2010

File size: 10.8 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 -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
22// cmvrvloscorf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
23// cmvrvloscorf -n 1,30 -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
24
25void usage(void);
26void usage(void)
27{
28cout
29<<"cmvrvloscor [options] rho.fits vlos.fits"<<endl
30<<"-n nplany,nhfill: process one Y plane every \"nplany\" (def:1(all))"<<endl
31<<" fill histos with \"nhfill\" los (def:25)"<<endl
32<<"-K npix: compute correlation R*V at +/- npix pixels (def: no)"<<endl
33<<" (very time consuming!!!)"<<endl
34<<"-S: compute cross-power spectrum of V*conj(R) (def: no)"<<endl
35<<"-N: do not create 3D cube and recompute 1D and 2D spectra (def: no do-it !)"<<endl
36<<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
37<<endl;
38}
39
40int main(int narg,char *arg[])
41{
42 int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
43 bool docube=true, dopk = false, do2d = false;
44
45 // --- Decodage des arguments
46 char c;
47 while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
48 switch (c) {
49 case 'n' :
50 sscanf(optarg,"%d,%d",&nplany,&nhfilllos);
51 if(nplany<=0) nplany = 1;
52 if(nhfilllos<=0) nhfilllos = 0;
53 break;
54 case 'K' :
55 npixcor = atoi(optarg);
56 break;
57 case 'S' :
58 dopk = true;
59 break;
60 case 'N' :
61 docube = false;
62 break;
63 case '2' :
64 do2d = true;
65 break;
66 case 'h' :
67 default :
68 usage(); return -1;
69 }
70 }
71 if(optind>=narg-1) {usage(); return -1;}
72
73 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
74 try {
75 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
76
77 SophyaInit();
78 InitTim();
79
80 // --- open FITS files (dRho/Rho and Vlos)
81 cout<<"> read rho: "<<arg[optind]<<endl;
82 FitsImg3DRead f3dr(arg[optind],0,5);
83 cout<<"> read vlos: "<<arg[optind+1]<<endl;
84 FitsImg3DRead f3dv(arg[optind+1],0,5);
85 long Nx = f3dr.ReadKeyL("Nx");
86 long Ny = f3dr.ReadKeyL("Ny");;
87 long Nz = f3dr.ReadKeyL("Nz");;
88 cout<<"N: x="<<Nx<<" y="<<Ny<<" z="<<Nz<<endl;
89 double Dx = f3dr.ReadKey("Dx");
90 double Dy = f3dr.ReadKey("Dy");
91 double Dz = f3dr.ReadKey("Dz");
92 cout<<"D: x="<<Dx<<" y="<<Dy<<" z="<<Dz<<endl;
93 double Zref = f3dr.ReadKey("ZREF");
94 double Href = f3dr.ReadKey("HREF");
95 cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
96
97 double dmin = min(Dx,min(Dy,Dz));
98 long nmax = max(Nx,max(Ny,Nz));
99 cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
100 Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);
101
102 POutPersist pos("cmvrvloscor.ppf");
103 DVList dvlcor;
104
105 // --- Create a Cube for analysis
106 GeneFluct3D *fluct3d = NULL;
107 TArray<GEN3D_TYPE>* rgen = NULL;
108 if(docube) {
109 cout<<"...Create and fill 3D cube"<<endl;
110 fluct3d = new GeneFluct3D(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
111 fluct3d->Print();
112 rgen = &(fluct3d->GetRealArray());
113 *rgen = 0.;
114 cout<<"rgen: size [1]="<<rgen->SizeX()
115 <<" [2]="<<rgen->SizeY()
116 <<" [3]="<<rgen->SizeZ()<<endl;
117 cout<<"pkgen: size [1]="<<fluct3d->GetComplexArray().SizeX()
118 <<" [2]="<<fluct3d->GetComplexArray().SizeY()
119 <<" [3]="<<fluct3d->GetComplexArray().SizeZ()<<endl;
120 }
121
122 // --- Vector for real-space correlation computation
123 int imil = Nz-1;
124 dvlcor("imil") = (int_4)imil;
125 TVector<int_4> nKsi;
126 TVector<r_8> Ksirv, Ksirvc, Ksirr, Ksirrc;
127 if(npixcor>0) {
128 Ksirv.ReSize(2*Nz-1); Ksirv = 0.;
129 Ksirvc.ReSize(2*Nz-1); Ksirvc = 0.;
130 Ksirr.ReSize(2*Nz-1); Ksirr = 0.;
131 Ksirrc.ReSize(2*Nz-1); Ksirrc = 0.;
132 nKsi.ReSize(2*Nz-1); nKsi = 0;
133 cout<<"...Compute R*V correlation on +/-"<<npixcor<<" px"<<endl;
134 }
135
136 // --- Vector for PK cross-correlation computation
137 int npk = 0;
138 TVector< complex<r_4> > FR, FV;
139 TVector< complex<r_8> > pkvr, FRdis;
140 TVector<r_8> pkr, pkrc;
141 FFTWServer fftserv;
142 if(dopk) cout<<"...compute V*conj(R) cross-power spectrum"<<endl;
143
144 // --- Read and process data
145 TVector<r_4> R(Nz), V(Nz);
146 TVector<r_8> Rdis(Nz);
147 if(nplany>Ny) nplany = Ny;
148 cout<<"...Will read one Y plane every "<<nplany<<endl;
149 if(nhfilllos) {
150 cout<<"...Fill Mpc displacement histo with "<<nhfilllos<<" los"<<endl;
151 nhfilllos = int((double)Nx*Ny/nplany/nhfilllos + 0.5);
152 if(nhfilllos<=0) nhfilllos = 1;
153 cout<<" -> fill one los every "<<nhfilllos<<endl;
154 }
155
156 cout<<">>> filling redshift distorted cube"<<endl;
157 int nlosread = 0;
158 for(int i=0;i<Nx;i++) {
159 if(i%(Nx/10)==0) cout<<" i="<<i<<endl;
160 TMatrix<r_4> M2d, M2dc;
161 if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
162 M2d.ReSize(Ny,Nz); M2d = 0.;
163 M2dc.ReSize(Ny,Nz); M2dc = 0.;
164 }
165 for(int j=0;j<Ny;j+=nplany) {
166 bool fhis = false;
167 if(nhfilllos) if(nlosread%nhfilllos==0) fhis = true;
168 //for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
169 //for(int l=0;l<Nz;l++) V(l) = f3dv.Read(l,j,i);
170 f3dr.Read(j,i,R);
171 f3dv.Read(j,i,V);
172 Rdis = 0.;
173 // Calcul du champ R redshift distordu
174 for(int l=0;l<Nz;l++) {
175 double d = (1.+Zref) / Href * V(l);
176 if(fhis) hmpc.Add(d);
177 double lpd = (double)l + d/Dz; // valeur du deplacee
178 // on repartit proportionellement au recouvrement sur 2 pixels
179 long l1 = long(lpd); // pixel de gauche
180 long l2 = l1 + 1; // pixel de droite
181 lpd -= (double)l1; // recouvrement du pixel du dessus
182 if(l1>=0 && l1<Nz) Rdis(l1) += R(l) * (1.-lpd);
183 if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
184 }
185 // On remplit eventuellement les matrices 2D
186 if(do2d && M2d.Size()>0)
187 for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dc(j,l) = Rdis(l);}
188 // On remplit le cube avec le champ R redshift distordu
189 if(fluct3d) for(int l=0;l<Nz;l++) (*rgen)(l,j,i) += Rdis(l);
190 // Calcul eventuel de la fonction de correlation R*V
191 if(npixcor>0) {
192 for(long l1=0;l1<Nz;l1++) {
193 for(long l2=max(0L,l1-npixcor);l2<min(Nz,l1+npixcor);l2++) {
194 int lc = imil+(l2-l1);
195 Ksirr(lc) += R(l1)*R(l2);
196 Ksirrc(lc) += Rdis(l1)*R(l2);
197 Ksirv(lc) += R(l1)*V(l2);
198 Ksirvc(lc) += Rdis(l1)*V(l2);
199 nKsi(lc)++;
200 }
201 }
202 }
203 // Cross-power spectrum computation
204 if(dopk) {
205 fftserv.FFTForward(V,FV);
206 int nf = FV.Size();
207 if(pkvr.Size()<=0) {
208 cout<<"...Create vector for cross-power spectrum computation"<<endl;
209 pkvr.ReSize(nf); pkvr = complex<r_8>(0.);
210 pkr.ReSize(nf); pkr = 0.;
211 pkrc.ReSize(nf); pkrc = 0.;
212 }
213 fftserv.FFTForward(R,FR);
214 for(int l=0;l<nf;l++) {
215 pkvr(l) += FV(l)*conj(FR(l));
216 pkr(l) += norm(FR(l));
217 }
218 fftserv.FFTForward(Rdis,FRdis);
219 for(int l=0;l<nf;l++) pkrc(l) += norm(FRdis(l));
220 npk++;
221 }
222 nlosread++;
223 }
224 if(do2d && M2d.Size()>0) {
225 char str[64];
226 sprintf(str,"mx_%d",i);
227 pos.PutObject(M2d,str);
228 sprintf(str,"mxc_%d",i);
229 pos.PutObject(M2dc,str);
230 }
231 }
232
233 cout<<"Number of processed los: "<<nlosread<<" / "<<Nx*Ny<<endl;
234 dvlcor("nlosread") = (int_4)nlosread;
235 if(hmpc.NEntries()>0) {
236 hmpc.Show();
237 pos.PutObject(hmpc,"hmpc");
238 }
239 if(Ksirr.Size()>0) {
240 for(int l=0;l<Ksirr.Size();l++) if(nKsi(l)>0) {
241 Ksirr(l) /= nKsi(l);
242 Ksirrc(l) /= nKsi(l);
243 Ksirv(l) /= nKsi(l);
244 Ksirvc(l) /= nKsi(l);
245 }
246 pos.PutObject(Ksirr,"ksirr");
247 pos.PutObject(Ksirrc,"ksirrc");
248 pos.PutObject(Ksirv,"ksirv");
249 pos.PutObject(Ksirvc,"ksirvc");
250 pos.PutObject(nKsi,"nksi");
251 }
252 if(npk>0) {
253 pkvr /= (double)npk;
254 pkr /= (double)npk;
255 pkrc /= (double)npk;
256 pos.PutObject(pkvr,"pkvr");
257 pos.PutObject(pkr,"pkr");
258 pos.PutObject(pkrc,"pkrc");
259 }
260 PrtTim(">>>> End filling redshift distorted cube");
261
262 // --- Fourier transform 3D cube and compute 1D and 2D power spectra
263 if(fluct3d) {
264 cout<<">>> Fourier transform 3D cube and compute 1D and 2D power spectra"<<endl;
265 // do the FFT for spectrum analysis
266 fluct3d->ReComputeFourier();
267 PrtTim(">>>> End ReComputing spectrum");
268
269 // Compute 1D spectrum
270 cout<<endl<<"\n--- Computing final 1D spectrum"<<endl;
271 double dkmin = fluct3d->GetKincMin(), knyqmax = fluct3d->GetKmax();
272 long nherr = long(knyqmax/dkmin+0.5);
273 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
274 HistoErr hpkrec(0.,knyqmax,nherr); hpkrec.Zero();
275 hpkrec.ReCenterBin(); hpkrec.Show();
276 fluct3d->ComputeSpectrum(hpkrec);
277 pos.PutObject(hpkrec,"hpkrec");
278 PrtTim(">>>> End Computing final spectrum");
279
280 // Compute 2D spectrum
281 cout<<"\n--- Computing final 2D spectrum"<<endl;
282 double dktmin = fluct3d->GetKTincMin(), ktnyqmax = fluct3d->GetKTmax();
283 double dkzmin = fluct3d->GetKinc()[2], kznyqmax = fluct3d->GetKnyq()[2];
284 long nherrt = long(ktnyqmax/dktmin+0.5), nherrz = long(kznyqmax/dkzmin+0.5);
285 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
286 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
287 Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
288 hpkrec2.ReCenterBin(); hpkrec2.Zero(); hpkrec2.Show();
289 fluct3d->ComputeSpectrum2D(hpkrec2);
290 pos.PutObject(hpkrec2,"hpkrec2");
291 PrtTim(">>>> End Computing final 2D spectrum");
292 }
293
294 // --- end of job, write objects in ppf
295 pos.PutObject(dvlcor,"dvlcor");
296 if(fluct3d) delete fluct3d;
297
298 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
299 } catch (PException& exc) {
300 cerr<<"cmvrvloscor.cc catched PException"<<exc.Msg()<<endl;
301 return 77;
302 } catch (std::exception& sex) {
303 cerr << "cmvrvloscor.cc std::exception :"
304 << (string)typeid(sex).name() << "\n msg= "
305 << sex.what() << endl;
306 return 78;
307 } catch (...) {
308 cerr << "cmvrvloscor.cc catched unknown (...) exception " << endl;
309 return 79;
310 }
311 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
312
313 return 0;
314}
315
316/*
317openppf cmvrvloscor.ppf
318
319disp hmpc
320
321# cross-correlation
322disp nksi
323set imil ${dvlcor.imil}
324n/plot ksirv.val%n-${imil} ! ! "nsta cpts"
325n/plot ksirvc.val%n-${imil} ! ! "nsta cpts same red"
326
327n/plot ksirr.val%n-${imil} ! ! "nsta cpts"
328n/plot ksirrc.val%n-${imil} ! ! "nsta cpts same red"
329
330# cross-power spectrum
331n/plot pkr.val%n ! ! "nsta cpts logx"
332n/plot pkrc.val%n ! ! "nsta cpts logx same red"
333
334n/plot pkvr.val%n ! ! "nsta cpts logx"
335
336# reconstructed 1D power spectrum
337n/plot hpkrec.val%x x>0 ! "nsta cpts logx"
338
339# reconstructed 2D power spectrum
340imag hpkrec2
341addoval 0 0 0.05 0.05 "green" false
342addoval 0 0 0.1 0.1 "green" false
343addoval 0 0 0.25 0.25 "green" false
344addoval 0 0 0.5 0.5 "green" false
345x = ${hpkrec2.xmax} / 2.
346addoval 0 0 $x $x "green" false
347x = ${hpkrec2.ymax} / 2.
348addoval 0 0 $x $x "green" false
349
350# proj selon kT (black), selon kZ (red)
351n/plot hpkrec2.val%sqrt(x*x+y*y) ! ! "nsta crossmarker3 logx"
352
353# les matrices 2D
354set n 0
355disp mx_$n
356newwin
357disp mxc_$n
358 */
Note: See TracBrowser for help on using the repository browser.