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

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

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