Changeset 3790 in Sophya for trunk/Cosmo/SimLSS/cmvrvloscor.cc

Timestamp:

Jun 28, 2010, 2:12:24 PM (15 years ago)

Author:

cmv

Message:

etude des cubes generes par le GSM, cmv 28/06/2010

File:

• trunk/Cosmo/SimLSS/cmvrvloscor.cc (modified) (10 diffs)

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -20 +20 @@
 // set simul = 6_0p0_780
 // cmvrvloscorf -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
-// cmvrvloscorf -n 1,30 -K 75 -S ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
+// cmvrvloscorf -n 1,30 -K 75 -S -2 ginit3d_${simul}_r.fits ginit3d_${simul}_rv.fits
 
 void usage(void);
@@ -33 +33 @@
 <<"-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
+<<"-2: compute 2D projection fpr dRho and dRho(corrected) (def=no)"<<endl
 <<endl;
 }
@@ -39 +40 @@
 {
  int nthread = 1, nplany=1, nhfilllos = 25, npixcor = 0;
- bool docube=true, dopk = false;
+ bool docube=true, dopk = false, do2d = false;
  
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"hn:K:SN")) != -1) {
+ while((c = getopt(narg,arg,"hn:K:SN2")) != -1) {
   switch (c) {
   case 'n' :
@@ -58 +59 @@
   case 'N' :
     docube = false;
+    break;
+  case '2' :
+    do2d = true;
     break;
   case 'h' :
@@ -93 +97 @@
  double nmax = max(Nx,max(Ny,Nz));
  cout<<"dmin="<<dmin<<" nmax="<<nmax<<endl;
- Histo hmpc(-dmin*nmax,dmin*nmax,4.*nmax);
+ Histo hmpc(-dmin*nmax,dmin*nmax,4*nmax);
 
  POutPersist pos("cmvrvloscor.ppf");
@@ -107 +111 @@
    rgen = &(fluct3d->GetRealArray());
    *rgen = 0.;
+   cout<<"rgen:  size [1]="<<rgen->SizeX()
+                           <<" [2]="<<rgen->SizeY()
+                           <<" [3]="<<rgen->SizeZ()<<endl;
+   cout<<"pkgen: size [1]="<<fluct3d->GetComplexArray().SizeX()
+                           <<" [2]="<<fluct3d->GetComplexArray().SizeY()
+                           <<" [3]="<<fluct3d->GetComplexArray().SizeZ()<<endl;
  }
 
@@ -147 +157 @@
  for(int i=0;i<Nx;i++) {
    if(i%(Nx/10)==0) cout<<"   i="<<i<<endl;
+   TMatrix<r_4> M2d, M2dc;
+   if(do2d && (i==0 || i==Nx/2 || i==Nx-1)) {
+     M2d.ReSize(Ny,Nz); M2d = 0.;
+     M2dc.ReSize(Ny,Nz); M2dc = 0.;
+   }
  for(int j=0;j<Ny;j+=nplany) {
    bool fhis = false;
@@ -161 +176 @@
      double lpd = (double)l + d/Dz; // valeur du deplacee
      // on repartit proportionellement au recouvrement sur 2 pixels
-     long l1 = long(lpd); // pixel de droite
-     long l2 = l1 + 1;  // pixel de gauche
+     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<Nz) Rdis(l1) += R(l) * (1.-lpd);
      if(l2>=0 && l2<Nz) Rdis(l2) += R(l) * lpd;
    }
+   // On remplit eventuellement les matrices 2D
+   if(do2d && M2d.Size()>0)
+     for(int l=0;l<Nz;l++) {M2d(j,l) = R(l); M2dc(j,l) = 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);
@@ -203 +221 @@
    nlosread++;
  }
+   if(do2d && M2d.Size()>0) {
+     char str[64];
+     sprintf(str,"mx_%d",i);
+     pos.PutObject(M2d,str);
+     sprintf(str,"mxc_%d",i);
+     pos.PutObject(M2dc,str);
+   }
  }
 
@@ -324 +349 @@
 # 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
  */