Changeset 3141 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.cc

Timestamp:

Jan 17, 2007, 6:44:04 PM (19 years ago)

Author:

cmv

Message:

chgt HProf->HistoErr + spectre 2D cmv 17/01/2007

File:

• trunk/Cosmo/SimLSS/cmvobserv3d.cc (modified) (8 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -10 +10 @@
 #include "ntuple.h"
 #include "matharr.h"
-#include "srandgen.h"
-#include "perandom.h"
-#include "fabtwriter.h"
-
-#include "arrctcast.h"
 
 #include "constcosmo.h"
@@ -208 +203 @@
  pkz.SetZ(zref);
  TArray< complex<r_8> > pkgen;
- GeneFluct3D fluct3d(pkgen,pkz);
+ GeneFluct3D fluct3d(pkgen);
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ TArray<r_8>& rgen = fluct3d.GetRealArray();
  fluct3d.Print();
+
+ double dkmin = fluct3d.GetKincMin();
  double knyqmax = fluct3d.GetKmax();
- double dkmin = fluct3d.GetKinc()[0];
+ long nherr = long(knyqmax/dkmin+0.5);
+ cout<<"For HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
+
+ double dktmin = fluct3d.GetKTincMin();
+ double ktnyqmax = fluct3d.GetKTmax();
+ long nherrt = long(ktnyqmax/dktmin+0.5);
+ double dkzmin = fluct3d.GetKinc()[2];
+ double kznyqmax = fluct3d.GetKnyq()[2];
+ long nherrz = long(kznyqmax/dkzmin+0.5);
+ cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
+     <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
 
  cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
@@ -219 +227 @@
  // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
  // Vc/Vs = 0.675   ->  keff = kmax * (0.675)^(1/3) = kmax * 0.877
- knyqmax *= 0.877;
- ldlkint = (log10(knyqmax)-log10(k1int))/npt;
+ double knyqmax_mod = 0.877*knyqmax;
+ ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax);
- cout<<"varpk_int(<"<<knyqmax<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
+ double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
+ cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
 
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
- if(computefourier0) fluct3d.ComputeFourier0();
-   else              fluct3d.ComputeFourier();
-
- cout<<"\n--- Checking realization spectra"<<endl;
- long nhprof = long(fluct3d.GetKmax()/dkmin+0.5);
- HProf hpkgen(0.,fluct3d.GetKmax(),nhprof);
- hpkgen.ReCenterBin();
- fluct3d.ComputeSpectrum(hpkgen);
- {
- tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
+ if(computefourier0) fluct3d.ComputeFourier0(pkz);
+   else              fluct3d.ComputeFourier(pkz);
+
+ if(1) {
+   cout<<"\n--- Checking realization spectra"<<endl;
+   HistoErr hpkgen(0.,knyqmax,nherr);
+   hpkgen.ReCenterBin(); hpkgen.Zero();
+   hpkgen.Show();
+   fluct3d.ComputeSpectrum(hpkgen);
+   {
+   tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
+   }
+ }
+
+ if(1) {
+   cout<<"\n--- Checking realization 2D spectra"<<endl;
+   Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkgen2.ReCenterBin(); hpkgen2.Zero();
+   hpkgen2.Show();
+   fluct3d.ComputeSpectrum2D(hpkgen2);
+   {
+   tagobs = "hpkgen2"; posobs.PutObject(hpkgen2,tagobs);
+   }
  }
 
@@ -241 +262 @@
    cout<<"\n--- Computing convolution by pixel shape"<<endl;
    fluct3d.FilterByPixel();
-
+ }
+
+ if(1) fluct3d.WriteFits("!cmvobserv3d_k.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_k.ppf",false);
+
+ if(1) {
    cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
-   HProf hpkgenf(hpkgen); hpkgenf.Zero();
+   HistoErr hpkgenf(0.,knyqmax,nherr);
+   hpkgenf.ReCenterBin(); hpkgenf.Zero();
+   hpkgenf.Show();
    fluct3d.ComputeSpectrum(hpkgenf);
    {
@@ -250 +278 @@
  }
 
- if(0) {
-   cout<<"\n--- Writing cmvobserv3dk.ppf"<<endl;
-   string tag = "cmvobserv3dk.ppf";
-   POutPersist pos(tag);
-   tag = "pkgen"; pos.PutObject(pkgen,tag);
+ if(1) {
+   cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
+   Histo2DErr hpkgenf2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkgenf2.ReCenterBin(); hpkgenf2.Zero();
+   hpkgenf2.Show();
+   fluct3d.ComputeSpectrum2D(hpkgenf2);
+   {
+   tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs);
+   }
  }
 
  cout<<"\n--- Computing a realization in real space"<<endl;
  fluct3d.ComputeReal();
- r_8 undouble=0.;
- TArray<r_8> rgen = ArrayCast(pkgen,undouble);
  double rmin,rmax; rgen.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
 
- cout<<"\n--- Check mean and variance in real space"<<endl;
- int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
- double rm,rs2; int_8 nm;
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
+ if(1) fluct3d.WriteFits("!cmvobserv3d_r0.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_r0.ppf",true);
+
+ int_8 nm;
+ double rm,rs2;
+ if(1) {
+   cout<<"\n--- Check mean and variance in real space"<<endl;
+   int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
+   nm = fluct3d.MeanSigma2(rm,rs2);
+   cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
+       <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
+ }
 
  if(1) {
@@ -331 +367 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- if(1) {
-   cout<<"\n---Writing Cube to Fits file"<<endl;
-   FitsImg3DWriter fwrt("!cmvobserv3dr.fits",FLOAT_IMG,5);
-   fwrt.WriteKey("ZREF",zref," redshift");
-   vector<long> n = fluct3d.GetNpix();
-   fwrt.WriteKey("NX",n[0]," axe transverse 1");
-   fwrt.WriteKey("NY",n[1]," axe transverse 2");
-   fwrt.WriteKey("NZ",n[2]," axe longitudinal (redshift)");
-   vector<r_8> d;
-   d = fluct3d.GetDinc();
-   fwrt.WriteKey("DX",d[0]," Mpc");
-   fwrt.WriteKey("DY",d[1]," Mpc");
-   fwrt.WriteKey("DZ",d[2]," Mpc");
-   d = fluct3d.GetKinc();
-   fwrt.WriteKey("DKX",d[0]," Mpc^-1");
-   fwrt.WriteKey("DKY",d[1]," Mpc^-1");
-   fwrt.WriteKey("DKZ",d[2]," Mpc^-1");
-   fwrt.WriteKey("SNOISE",snoise," Msol");
-   fwrt.Write(rgen);
- }
+ if(1) fluct3d.WriteFits("!cmvobserv3d_r.fits");
+ if(1) fluct3d.WritePPF("cmvobserv3d_r.ppf",true);
 
  cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
@@ -358 +376 @@
      <<rs2<<" -> "<<sqrt(rs2)<<endl;
 
- if(0) {
-   cout<<"\n--- Writing cmvobserv3dr.ppf"<<endl;
-   string tag = "cmvobserv3dr.ppf";
-   POutPersist pos(tag);
-   tag = "rgen"; pos.PutObject(rgen,tag);
- }
-
  //-----------------------------------------------------------------
  // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
+ 
  if(1) {
    cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
    fluct3d.ReComputeFourier();
-   HProf hpkrec(0.,fluct3d.GetKmax(),nhprof);
+ }
+
+ if(1) {
+   cout<<endl<<"\n--- Computing final spectrum"<<endl;
+   HistoErr hpkrec(0.,knyqmax,nherr);
    hpkrec.ReCenterBin();
+   hpkrec.Show();
    fluct3d.ComputeSpectrum(hpkrec);
    tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
  }
 
+ if(1) {
+   cout<<"\n--- Computing final 2D spectrum"<<endl;
+   Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkrec2.ReCenterBin(); hpkrec2.Zero();
+   hpkrec2.Show();
+   fluct3d.ComputeSpectrum2D(hpkrec2);
+   {
+   tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
+   }
+ }
+
  return 0;
 }
 
 /*
-openppf cmvobserv3dk.ppf
-openppf cmvobserv3dr.ppf
+######################################################
+readfits cmvobserv3d_k.fits
+readfits cmvobserv3d_r0.fits
+readfits cmvobserv3d_r.fits
+
+openppf cmvobserv3d_k.ppf
+openppf cmvobserv3d_r0.ppf
+openppf cmvobserv3d_r.ppf
+
+# pour le plot 2D d'une slice en Z du 3D: to2d nom_obj3D num_slice
+defscript to2d
+ objaoper $1 sliceyz $2
+ mv sliceyz_${2} ${1}_$2
+ disp ${1}_$2
+ echo display slice $2 of $1
+endscript 
+
+to2d $cobj 0
+
+######################################################
 openppf cmvobserv3d.ppf
 
-objaoper pkgen sliceyz 0
-
+zone
 n/plot hpkz.val%x ! ! "nsta connectpoints"
 n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
@@ -394 +439 @@
 n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
 
-defscript rgensl
- objaoper rgen sliceyz $1
- disp sliceyz_$1 
-endscript 
-
-rgensl 0
-
+zone 2 2
+imag hpkgen2
+imag hpkgenf2
+imag hpkrec2
+
+zone
 disp hmdndm
 disp tirhmdndm