Changeset 3771 in Sophya

Timestamp:

May 8, 2010, 12:59:59 AM (15 years ago)

Author:

cmv

Message:

add proportionnal filling in pixels, cmv 08/05/2010

File:

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

trunk/Cosmo/SimLSS/cmvrvloscor.cc

@@ -27 +27 @@
 {
  int nthread = 1;
- int_8 ntfill = 100000;
  SophyaInit();
  InitTim();
@@ -41 +40 @@
  cout<<"> read vlos: "<<arg[2]<<endl;
  FitsImg3DRead f3dv(arg[2],0,5);
- long naxis1 = f3dr.ReadKeyL("NAXIS1");
- long naxis2 = f3dr.ReadKeyL("NAXIS2");
- long naxis3 = f3dr.ReadKeyL("NAXIS3");
- cout<<"Naxis: 1="<<naxis1<<" 2="<<naxis2<<" 3="<<naxis3<<endl;
- long nx = f3dr.ReadKeyL("Nx");
- long ny = f3dr.ReadKeyL("Ny");;
- long nz = f3dr.ReadKeyL("Nz");;
- cout<<"N: x="<<nx<<" y="<<ny<<" z="<<nz<<endl;
- double dx = f3dr.ReadKey("Dx");
- double dy = f3dr.ReadKey("Dy");
- double dz = f3dr.ReadKey("Dz");
- cout<<"D: x="<<dx<<" y="<<dy<<" z="<<dz<<endl;
- double zref = f3dr.ReadKey("ZREF");
- double href = f3dr.ReadKey("HREF");
- cout<<"Zref="<<zref<<" Href="<<href<<endl;
+ long Nx = f3dr.ReadKeyL("Nx");
+ long Ny = f3dr.ReadKeyL("Ny");;
+ long Nz = f3dr.ReadKeyL("Nz");;
+ cout<<"N: x="<<Nx<<" y="<<Ny<<" z="<<Nz<<endl;
+ double Dx = f3dr.ReadKey("Dx");
+ double Dy = f3dr.ReadKey("Dy");
+ double Dz = f3dr.ReadKey("Dz");
+ cout<<"D: x="<<Dx<<" y="<<Dy<<" z="<<Dz<<endl;
+ double Zref = f3dr.ReadKey("ZREF");
+ double Href = f3dr.ReadKey("HREF");
+ cout<<"Zref="<<Zref<<" Href="<<Href<<endl;
 
  POutPersist pos("cmvrvloscor.ppf");
 
- GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
+ GeneFluct3D fluct3d(Nx,Ny,Nz,Dx,Dy,Dz,nthread,2);
  fluct3d.Print();
  TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
  rgen = 0.;
- TVector<r_8> R(nz), Rdis(nz);
-
- const int nvar = 3;
- const char *vname[nvar] = {"l","ll","d"};
- NTuple nt(nvar,vname);
- r_4 xnt[nvar];
- int_8 ntmod = fluct3d.NPix()/ntfill; if(ntmod==0) ntmod = 1;
+ TVector<r_8> R(Nz), Rdis(Nz);
 
  cout<<"> filling redshift distorted cube"<<endl;
- int_8 nf = 0;
- for(int i=0;i<nx;i++) {
-   if(i%(nx/10)==0) cout<<"i="<<i<<endl;
- for(int j=0;j<ny;j++) {
-   for(int l=0;l<nz;l++) R(l) = f3dr.Read(l,j,i);
+ for(int i=0;i<Nx;i++) {
+   if(i%(Nx/10)==0) cout<<"i="<<i<<endl;
+ for(int j=0;j<Ny;j++) {
+   for(int l=0;l<Nz;l++) R(l) = f3dr.Read(l,j,i);
    Rdis = 0.;
-   for(int l=0;l<nz;l++) {
-     double d = (1.+zref) / href * f3dv.Read(l,j,i);
-     int ll = int((double)l + d/dz + 0.5);
-     if(ll<0 || ll>=nz) continue;
-     Rdis(ll) += R(l);
-     if(ntfill>0 && nf%ntmod==0)
-       {xnt[0]=l; xnt[1]=ll; xnt[2]=d; nt.Fill(xnt);}
-     nf++;
+   for(int l=0;l<Nz;l++) {
+     double d = (1.+Zref) / Href * f3dv.Read(l,j,i);
+     double lpd = (double)l + d/Dz; // valeur du deplacee
+     // on repartit proportionnelement au recouvrement sur 2 pixels
+     int l1 = int(lpd); // pixel de droite
+     int l2 = l1 + 1;  // pixel de gauche
+     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;
    }
-   for(int l=0;l<nz;l++) rgen(l,j,i) += Rdis(l);
+   for(int l=0;l<Nz;l++) rgen(l,j,i) += Rdis(l);
  }
  }
  PrtTim(">>>> End filling redshift distorted cube");
- pos.PutObject(nt,"nt");
 
  fluct3d.ReComputeFourier();
@@ -146 +134 @@
 
 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 cpts"
-
  */