Changeset 3349 in Sophya

Timestamp:

Oct 11, 2007, 4:39:58 PM (18 years ago)

Author:

cmv

Message:

• gros changements dans la structure de la classe GeneFluct3D (constructeur et logique d'aloocation memoire, init_fftw etc...)
• suppression des valeurs de masse<0 mises a -999. directement mises a zero
• suppression de TurnMass2HIMass qui fait maintenant la meme chose que TurnMass2MeanNumber
• legere restructuration de cmvobserv3d.cc pour compat.

cmv 11/10/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (15 diffs)
• genefluct3d.cc (modified) (8 diffs)
• genefluct3d.h (modified) (5 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -45 +45 @@
      <<" -R schmassdist.ppf : read mass distribution for trials from file"<<endl
      <<"               instead of computing it (ONLY if \"-Q\" option is activated)"<<endl
-     <<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
-     <<"    AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
-     <<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
      <<" -W : write cube in FITS format (complex cube is coded as real cube)"<<endl
      <<" -P : write cube in PPF format"<<endl
@@ -54 +51 @@
      <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
      <<endl;
+     ////<<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
+     ////<<"    AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
+     ////<<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
 }
 
@@ -101 +101 @@
  int noise_evol = 0;
 
- // *** AGN
- bool do_agn = false;
- double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
- double powlaw_agn = 0.;
+ //// *** AGN
+ ////bool do_agn = false;
+ ////double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
+ ////double powlaw_agn = 0.;
 
  // *** type de generation
@@ -177 +177 @@
     schmassdistfile = optarg;
     break;
-  case 'A' :
-    do_agn = true;
-    sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
-    break;
+    ////  case 'A' :
+    ////do_agn = true;
+    ////sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
+    ////break;
   case 'V' :
     compvarreal = true;
@@ -224 +224 @@
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
  cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
- if(do_agn)
-   cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
-       <<" , powlaw="<<powlaw_agn<<endl;
+ ////if(do_agn)
+ ////  cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+ ////      <<" , powlaw="<<powlaw_agn<<endl;
 
  string tagobs = "cmvobserv3d.ppf";
@@ -264 +264 @@
  cout<<endl<<"\n--- Compute variance for top-hat R="<<R
      <<" at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_th(pkz,0);
- double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
+ VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
+ double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
  double pkmax_th = varpk_th(kfind_th);
  cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
  double k1=kmin, k2=kmax;
- int rc = varpk_th.FindLimits(R,pkmax_th/1.e4,k1,k2,eps);
+ int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
  cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
      <<k2<<" ("<<log10(k2)<<")"<<endl;
@@ -275 +275 @@
  double ldlk = (log10(k2)-log10(k1))/npt;
  varpk_th.SetInteg(0.01,ldlk,-1.,4);
- double sr2 = varpk_th.Variance(R,k1,k2);
+ double sr2 = varpk_th.Variance(k1,k2);
  cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
 
@@ -292 +292 @@
  //-----------------------------------------------------------------
  cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_int(pkz,2);
-
- double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
+ VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
+
+ double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
  double pkmax_int = varpk_int(kfind_int);
  cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
  double k1int=kmin, k2int=kmax;
- int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
+ int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
  cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
      <<k2int<<" ("<<log10(k2int)<<")"<<endl;
@@ -304 +304 @@
  double ldlkint = (log10(k2int)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int = varpk_int.Variance(R,k1int,k2int);
+ double sr2int = varpk_int.Variance(k1int,k2int);
  cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
  
@@ -355 +355 @@
  cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
 
- TArray< complex<r_8> > pkgen;
- GeneFluct3D fluct3d(pkgen);
- fluct3d.SetPrtLevel(2);
- fluct3d.SetNThread(nthread);
- fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
  fluct3d.SetObservator(zref,-nz/2.);
  fluct3d.SetCosmology(univ);
  fluct3d.SetGrowthFactor(growth);
  fluct3d.LosComRedshift(0.001,-1);
- //TArray<r_8>& rgen = fluct3d.GetRealArray();
+ TArray< complex<r_8> >& pkgen = fluct3d.GetComplexArray();
+ TArray<r_8>& rgen = fluct3d.GetRealArray();
  cout<<endl; fluct3d.Print();
  cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
@@ -391 +388 @@
  ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
+ double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
  cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
 
@@ -538 +535 @@
 
  if(no_poisson) {
-
    cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
        <<mass_by_mpc3*fluct3d.GetDVol()<<endl;
-   rm = fluct3d.TurnMass2HIMass(mass_by_mpc3);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
-   PrtTim(">>>> End Converting !!!DIRECTLY!!! mass into HI mass");
-
+   rm = fluct3d.TurnMass2MeanNumber(mass_by_mpc3);
  } else {
-
    cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
        <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
    rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
-   PrtTim(">>>> End Converting mass into galaxy number");
-
-   cout<<"\n--- Set negative and null pixels to BAD"<<endl;
-   nm = fluct3d.SetToVal(0.,1e+200,-999.);
-   PrtTim(">>>> End Set negative pixels to BAD etc...");
+ }
+ nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ PrtTim(">>>> End Converting mass into galaxy number or mass");
+
+ if( !no_poisson ) {
 
    cout<<"\n--- Apply poisson on galaxy number"<<endl;
    fluct3d.ApplyPoisson();
-     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2);
    double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
    PrtTim(">>>> End Apply poisson on galaxy number");
@@ -584 +574 @@
    }
    cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
-     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
      cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+     nm = fluct3d.MeanSigma2(rm,rs2);
    PrtTim(">>>> End Convert Galaxy number to HI mass");
-
-   cout<<"\n--- Set BAD pixels to Zero"<<endl;
-   nm = fluct3d.SetToVal(-998.,1e+200,0.);
-   PrtTim(">>>> End Set BAD pixels to Zero etc...");
 
  }
@@ -609 +595 @@
 
  //-----------------------------------------------------------------
- if(do_agn) {
-   cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
-       <<" , powlaw="<<powlaw_agn<<endl;
-   fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
-     nm = fluct3d.MeanSigma2(rm,rs2);
-   PrtTim(">>>> End Add AGN");
- }
+ ////if(do_agn) {
+ ////  cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+ ////      <<" , powlaw="<<powlaw_agn<<endl;
+ ////  fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
+ ////    nm = fluct3d.MeanSigma2(rm,rs2);
+ ////   PrtTim(">>>> End Add AGN");
+ ////}
 
  //-----------------------------------------------------------------
@@ -627 +613 @@
  }
 
+
+ //-----------------------------------------------------------------
  if(scalecube!=0.) {  // Si scalecube==0 pas de normalisation
    cout<<"\n--- Scale cube rs2ref="<<rs2ref<<endl;

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -24 +24 @@
 #include "genefluct3d.h"
 
-#define FFTW_THREAD
+#define WITH_FFTW_THREAD
 
 #define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
@@ -31 +31 @@
 
 //-------------------------------------------------------
-GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : Nx_(0) , Ny_(0) , Nz_(0)
-  , lp_(0)
-  , array_allocated_(false) , T_(T)
-  , cosmo_(NULL) , growth_(NULL)
-  , redsh_ref_(-999.), kredsh_ref_(0.), dred_ref_(-999.)
-  , loscom_ref_(-999.), dtrc_ref_(-999.), dlum_ref_(-999.), dang_ref_(-999.)
-  , nu_ref_(-999.), dnu_ref_ (-999.)
-  , loscom_min_(-999.), loscom_max_(-999.)
-  , loscom2zred_min_(0.) , loscom2zred_max_(0.)
-{
+GeneFluct3D::GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz
+                        ,unsigned short nthread,int lp)
+{
+ init_default();
+
+ lp_ = lp;
+ nthread_ = nthread;
+
+ setsize(nx,ny,nz,dx,dy,dz);
+ setalloc();
+ setpointers(false);
+ init_fftw();
+}
+
+GeneFluct3D::GeneFluct3D(unsigned short nthread)
+{
+ init_default();
+ setsize(2,2,2,1.,1.,1.);
+ nthread_ = nthread;
+ setalloc();
+ setpointers(false);
+ init_fftw();
+}
+
+GeneFluct3D::~GeneFluct3D(void)
+{
+ delete_fftw();
+}
+
+//-------------------------------------------------------
+void GeneFluct3D::init_default(void)
+{
+ Nx_ = Ny_ = Nz_ = 0;
+ is_set_fftw_plan = false;
+ nthread_ = 0;
+ lp_ = 0;
+ array_allocated_ = false;
+ cosmo_ = NULL;
+ growth_ = NULL;
+ redsh_ref_ = -999.;
+ kredsh_ref_ = 0.;
+ dred_ref_ = -999.;
+ loscom_ref_ = -999.;
+ dtrc_ref_ = dlum_ref_ = dang_ref_ = -999.;
+ nu_ref_ = dnu_ref_ = -999.;
+ loscom_min_ = loscom_max_ = -999.;
+ loscom2zred_min_ = loscom2zred_max_ = 0.;
  xobs_[0] = xobs_[1] = xobs_[2] = 0.;
  zred_.resize(0);
  loscom_.resize(0);
  loscom2zred_.resize(0);
- SetNThread();
-}
-
-GeneFluct3D::~GeneFluct3D(void)
-{
- fftw_destroy_plan(pf_);
- fftw_destroy_plan(pb_);
-#ifdef FFTW_THREAD
- if(nthread_>0) fftw_cleanup_threads();
-#endif
-}
-
-//-------------------------------------------------------
-void GeneFluct3D::SetSize(long nx,long ny,long nz,double dx,double dy,double dz)
-{
-  setsize(nx,ny,nz,dx,dy,dz);
-  setalloc();
-  setpointers(false);
-  init_fftw();
-}
-
-void GeneFluct3D::SetObservator(double redshref,double kredshref)
-// L'observateur est au redshift z=0
-//               est situe sur la "perpendiculaire" a la face x,y
-//                         issue du centre de cette face
-// Il faut positionner le cube sur l'axe des z cad des redshifts:
-//     redshref  = redshift de reference
-//                 Si redshref<0 alors redshref=0
-//     kredshref = indice (en double) correspondant a ce redshift
-//                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
-// Exemple: redshref=1.5 kredshref=250.75
-//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
-{
- if(redshref<0.) redshref = 0.;
- if(kredshref<0.) {
-   if(Nz_<=0) {
-     char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 SetSize should be called first";
-     cout<<bla<<endl; throw ParmError(bla);
-   }
-   kredshref = Nz_/2.;
- }
- redsh_ref_  = redshref;
- kredsh_ref_ = kredshref;
- if(lp_>0)
-   cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
-}
-
-void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
-{
- cosmo_ = &cosmo;
- if(lp_>1) cosmo_->Print();
-}
-
-void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
-{
- growth_ = &growth;
 }
 
@@ -171 +151 @@
 }
 
-void GeneFluct3D::check_array_alloc(void)
-// Pour tester si le tableau T_ est alloue
-{
- if(array_allocated_) return;
- char bla[90];
- sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
- cout<<bla<<endl; throw ParmError(bla);
-}
-
 void GeneFluct3D::init_fftw(void)
 {
+ if( is_set_fftw_plan ) delete_fftw();
+
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
  if(lp_>1) cout<<"--- GeneFluct3D::init_fftw ---"<<endl;
-#ifdef FFTW_THREAD
+#ifdef WITH_FFTW_THREAD
  if(nthread_>0) {
    cout<<"...Computing with "<<nthread_<<" threads"<<endl;
@@ -195 +168 @@
  if(lp_>1) cout<<"...backward plan"<<endl;
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
+ is_set_fftw_plan = true;
+}
+
+void GeneFluct3D::delete_fftw(void)
+{
+ if( !is_set_fftw_plan ) return;
+ fftw_destroy_plan(pf_);
+ fftw_destroy_plan(pb_);
+#ifdef WITH_FFTW_THREAD
+ if(nthread_>0) fftw_cleanup_threads();
+#endif
+ is_set_fftw_plan = false;
+}
+
+void GeneFluct3D::check_array_alloc(void)
+// Pour tester si le tableau T_ est alloue
+{
+ if(array_allocated_) return;
+ char bla[90];
+ sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
+ cout<<bla<<endl; throw ParmError(bla);
 }
 
 //-------------------------------------------------------
+void GeneFluct3D::SetObservator(double redshref,double kredshref)
+// L'observateur est au redshift z=0
+//               est situe sur la "perpendiculaire" a la face x,y
+//                         issue du centre de cette face
+// Il faut positionner le cube sur l'axe des z cad des redshifts:
+//     redshref  = redshift de reference
+//                 Si redshref<0 alors redshref=0
+//     kredshref = indice (en double) correspondant a ce redshift
+//                 Si kredshref<0 alors kredshref=nz/2 (milieu du cube)
+// Exemple: redshref=1.5 kredshref=250.75
+//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
+{
+ if(redshref<0.) redshref = 0.;
+ if(kredshref<0.) {
+   if(Nz_<=0) {
+     char *bla = "GeneFluct3D::SetObservator_Error: for kredsh_ref<0 define cube geometry first";
+     cout<<bla<<endl; throw ParmError(bla);
+   }
+   kredshref = Nz_/2.;
+ }
+ redsh_ref_  = redshref;
+ kredsh_ref_ = kredshref;
+ if(lp_>0)
+   cout<<"--- GeneFluct3D::SetObservator zref="<<redsh_ref_<<" kref="<<kredsh_ref_<<endl;
+}
+
+void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
+{
+ cosmo_ = &cosmo;
+ if(lp_>1) cosmo_->Print();
+}
+
+void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
+{
+ growth_ = &growth;
+}
+
 long GeneFluct3D::LosComRedshift(double zinc,long npoints)
 // Given a position of the cube relative to the observer
@@ -1164 +1195 @@
 }
 
-double  GeneFluct3D::TurnMass2HIMass(double m_by_mpc3)
-{
- if(lp_>0) cout<<"--- TurnMass2HIMass : "<<m_by_mpc3<<" Msol/Mpc^3"<<endl;
-
- double mall=0., mgood=0.;
- int_8  nall=0,  ngood=0;
- for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = IndexR(i,j,l);
-   mall += data_[ip]; nall++;
-   if(data_[ip]>0.) {mgood += data_[ip]; ngood++;}
- }
- if(ngood>0) mgood /= (double)ngood;
- if(nall>0)  mall  /= (double)nall;
- if(lp_>0) cout<<"...ngood="<<ngood<<" mgood="<<mgood
-               <<",  nall="<<nall<<" mall="<<mall<<endl;
- if(ngood<=0 || mall<=0.) {
-   cout<<"TurnMass2HIMass_Error: ngood="<<ngood<<" <=0 || mall="<<mall<<" <=0"<<endl;
-   throw RangeCheckError("TurnMass2HIMass_Error: ngood<=0 || mall<=0");
- }
-
- // On doit mettre m*Vol masse de HI dans notre survey
- // On en met uniquement dans les pixels de masse >0.
- // On MET a zero les pixels <0
- // On renormalise sur les pixels>0 pour qu'on ait m*Vol masse de HI
- //   comme on ne prend que les pixels >0, on doit normaliser
- //   a la moyenne de <1+d_rho/rho> sur ces pixels
- //   (rappel sur tout les pixels <1+d_rho/rho>=1)
- // masse de HI a mettre ds 1 px:
- double dm = m_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
- if(lp_>0) cout<<"...HI mass density move from "
-               <<m_by_mpc3*Vol_/double(NRtot_)<<" to "<<dm<<" / pixel"<<endl;
-
- double sum = 0.;
- for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
-   int_8 ip = IndexR(i,j,l);
-   if(data_[ip]<=0.) data_[ip] = 0.;
-   else {
-     data_[ip] *= dm; // par coherence on multiplie aussi les <=0
-     sum += data_[ip];  // mais on ne les compte pas
-   }
- }
-
- if(lp_>0) cout<<sum<<"...mass HI put into survey / "<<m_by_mpc3*Vol_<<endl;
-
- return sum;
-}
-
-double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
-// do NOT treate negative or nul values
-{
- if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<n_by_mpc3<<" gal/Mpc^3"<<endl;
+double GeneFluct3D::TurnMass2MeanNumber(double val_by_mpc3)
+{
+ if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;
 
  double mall=0., mgood=0.;
@@ -1232 +1215 @@
  }
 
- // On doit mettre n*Vol galaxies dans notre survey
+ // On doit mettre n*Vol galaxies (ou Msol) dans notre survey
  // On en met uniquement dans les pixels de masse >0.
- // On NE met PAS a zero les pixels <0
- // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies
+ // On MET a zero les pixels <0
+ // On renormalise sur les pixels>0 pour qu'on ait n*Vol galaxies (ou Msol)
  //   comme on ne prend que les pixels >0, on doit normaliser
  //   a la moyenne de <1+d_rho/rho> sur ces pixels
  //   (rappel sur tout les pixels <1+d_rho/rho>=1)
- // nb de gal a mettre ds 1 px:
- double dn = n_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
- if(lp_>0) cout<<"...galaxy density move from "
-               <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
+ // nb de gal (ou Msol) a mettre ds 1 px:
+ double dn = val_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
+ if(lp_>0) cout<<"...quantity density move from "
+               <<val_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
 
  double sum = 0.;
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
    int_8 ip = IndexR(i,j,l);
-   data_[ip] *= dn; // par coherence on multiplie aussi les <=0
-   if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
- }
-
- if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
+   if(data_[ip]<=0.) data_[ip] = 0.;
+   else {
+     data_[ip] *= dn;
+     sum += data_[ip];
+   }
+ }
+
+ if(lp_>0) cout<<sum<<"...quantity put into survey / "<<val_by_mpc3*Vol_<<endl;
 
  return sum;
@@ -1329 +1315 @@
 }
 
+void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
+// add noise to every pixels (meme les <=0 !)
+// type_evol = 0 : pas d'evolution de la puissance du bruit
+//             1 : evolution de la puissance du bruit avec la distance a l'observateur
+//             2 : evolution de la puissance du bruit avec la distance du plan Z
+//                 (tous les plans Z sont mis au meme redshift z de leur milieu)
+{
+ if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
+ check_array_alloc();
+
+ if(type_evol<0) type_evol = 0;
+ if(type_evol>2) {
+   char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
+   cout<<bla<<endl; throw ParmError(bla);
+ }
+
+ vector<double> correction;
+ InterpFunc *intercor = NULL;
+
+ if(type_evol>0) {
+   // Sigma_Noise(en mass) :
+   //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
+   //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
+   //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
+   //                                    nlobe ~ 1/Dtrcom^2
+   //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
+   if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
+     char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
+     cout<<bla<<endl; throw ParmError(bla);
+   }
+   InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
+   long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
+   for(long i=0;i<nsz;i++) {
+     double d = interpinv.X(i);
+     double zred = interpinv(d);
+     double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
+     double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
+     double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
+     double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
+     double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
+     if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
+       cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
+           <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
+     correction.push_back(corr);
+   }
+   intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
+ } 
+
+ double corrlim[2] = {1.,1.};
+ for(long i=0;i<Nx_;i++) {
+   double dx2 = DXcom(i); dx2 *= dx2;
+   for(long j=0;j<Ny_;j++) {
+     double dy2 = DYcom(j); dy2 *= dy2;
+     for(long l=0;l<Nz_;l++) {
+       double corr = 1.;
+       if(type_evol>0) {
+         double dz = DZcom(l);
+         if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
+           else dz = fabs(dz); // tous les plans Z au meme redshift
+         corr = (*intercor)(dz);
+         if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
+       }
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] += snoise*corr*NorRand();
+     }
+   }
+ }
+ if(type_evol>0)
+   cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
+
+ if(intercor!=NULL) delete intercor;
+}
+
+}  // Fin namespace SOPHYA
+
+
+
+
+/*********************************************************************
 void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
 // Add AGN flux into simulation:
@@ -1423 +1488 @@
  
 }
-
-void GeneFluct3D::AddNoise2Real(double snoise,int type_evol)
-// add noise to every pixels (meme les <=0 !)
-// type_evol = 0 : pas d'evolution de la puissance du bruit
-//             1 : evolution de la puissance du bruit avec la distance a l'observateur
-//             2 : evolution de la puissance du bruit avec la distance du plan Z
-//                 (tous les plans Z sont mis au meme redshift z de leur milieu)
-{
- if(lp_>0) cout<<"--- AddNoise2Real: snoise = "<<snoise<<" evol="<<type_evol<<endl;
- check_array_alloc();
-
- if(type_evol<0) type_evol = 0;
- if(type_evol>2) {
-   char *bla = "GeneFluct3D::AddNoise2Real_Error: bad type_evol value";
-   cout<<bla<<endl; throw ParmError(bla);
- }
-
- vector<double> correction;
- InterpFunc *intercor = NULL;
-
- if(type_evol>0) {
-   // Sigma_Noise(en mass) :
-   //      Slim ~ 1/sqrt(DNu) * sqrt(nlobe)   en W/m^2Hz
-   //      Flim ~ sqrt(DNu) * sqrt(nlobe)   en W/m^2
-   //      Mlim ~ sqrt(DNu) * (Dlum)^2 * sqrt(nlobe)  en Msol
-   //                                    nlobe ~ 1/Dtrcom^2
-   //      Mlim ~ sqrt(DNu) * (Dlum)^2 / Dtrcom
-   if(cosmo_ == NULL || redsh_ref_<0.| loscom2zred_.size()<1) {
-     char *bla = "GeneFluct3D::AddNoise2Real_Error: set Observator and Cosmology first";
-     cout<<bla<<endl; throw ParmError(bla);
-   }
-   InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
-   long nsz = loscom2zred_.size(), nszmod=((nsz>10)? nsz/10: 1);
-   for(long i=0;i<nsz;i++) {
-     double d = interpinv.X(i);
-     double zred = interpinv(d);
-     double dtrc = cosmo_->Dtrcom(zred);  // pour variation angle solide
-     double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
-     double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
-     double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
-     double corr = sqrt(dnu/dnu_ref_) * pow(dlum/dlum_ref_,2.) * dtrc_ref_/dtrc;
-     if(lp_>0 && (i==0 || i==nsz-1 || i%nszmod==0))
-       cout<<"i="<<i<<" d="<<d<<" red="<<zred<<" dred="<<dred<<" dnu="<<dnu
-           <<" dtrc="<<dtrc<<" dlum="<<dlum<<" -> cor="<<corr<<endl;
-     correction.push_back(corr);
-   }
-   intercor = new InterpFunc(loscom2zred_min_,loscom2zred_max_,correction);
- } 
-
- double corrlim[2] = {1.,1.};
- for(long i=0;i<Nx_;i++) {
-   double dx2 = DXcom(i); dx2 *= dx2;
-   for(long j=0;j<Ny_;j++) {
-     double dy2 = DYcom(j); dy2 *= dy2;
-     for(long l=0;l<Nz_;l++) {
-       double corr = 1.;
-       if(type_evol>0) {
-         double dz = DZcom(l);
-         if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
-           else dz = fabs(dz); // tous les plans Z au meme redshift
-         corr = (*intercor)(dz);
-         if(corr<corrlim[0]) corrlim[0]=corr; else if(corr>corrlim[1]) corrlim[1]=corr;
-       }
-       int_8 ip = IndexR(i,j,l);
-       data_[ip] += snoise*corr*NorRand();
-     }
-   }
- }
- if(type_evol>0)
-   cout<<"correction factor range: ["<<corrlim[0]<<","<<corrlim[1]<<"]"<<endl;
-
- if(intercor!=NULL) delete intercor;
-}
-
-}  // Fin namespace SOPHYA
+*********************************************************************/

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -24 +24 @@
 class GeneFluct3D {
 public:
-  GeneFluct3D(TArray< complex<r_8 > >& T);
+  GeneFluct3D(long nx,long ny,long nz,double dx,double dy,double dz,unsigned short nthread=0,int lp=0);  // Mpc
+  GeneFluct3D(unsigned short nthread=0);
   virtual ~GeneFluct3D(void);
 
-  void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
-  void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc
   // Distance los comobile a l'observateur
   void SetObservator(double redshref=0.,double kredshref=0.);
@@ -111 +110 @@
 
   void TurnFluct2Mass(void);
-  double TurnMass2HIMass(double m_by_mpc3);
-  double TurnMass2MeanNumber(double n_by_mpc3);
+  double TurnMass2MeanNumber(double val_by_mpc3);
   double ApplyPoisson(void);
   double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
   double TurnNGal2MassQuick(SchechterMassDist& schmdist);
   double TurnMass2Flux(void);
-  void AddAGN(double lfjy,double lsigma,double powlaw=0.);
+  //void AddAGN(double lfjy,double lsigma,double powlaw=0.);
   void AddNoise2Real(double snoise,int type_evol=0);
 
@@ -133 +131 @@
 
 protected:
+  void init_default(void);
   void setsize(long nx,long ny,long nz,double dx,double dy,double dz);
   void setalloc(void);
   void setpointers(bool from_real);
   void init_fftw(void);
+  void delete_fftw(void);
   long manage_coefficients(void);
   double compute_power_carte(void);
@@ -157 +157 @@
 
   // la gestion de la FFT
+  bool is_set_fftw_plan;
   fftw_plan pf_,pb_;
   unsigned short nthread_;
@@ -163 +164 @@
   // le stockage du Cube de donnees et les pointeurs
   bool array_allocated_;  // true if array has been allocated
-  TArray< complex<r_8> >& T_;
+  TArray< complex<r_8> > T_;
   fftw_complex *fdata_;
   TArray<r_8> R_;