Index: trunk/Cosmo/SimLSS/genefluct3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3331)
+++ trunk/Cosmo/SimLSS/genefluct3d.cc	(revision 3349)
@@ -24,5 +24,5 @@
 #include "genefluct3d.h"
 
-#define FFTW_THREAD
+#define WITH_FFTW_THREAD
 
 #define MODULE2(_x_) ((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
@@ -31,75 +31,55 @@
 
 //-------------------------------------------------------
-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,18 +151,11 @@
 }
 
-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,7 +168,65 @@
  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,55 +1195,7 @@
 }
 
-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,24 +1215,27 @@
  }
 
- // 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,4 +1315,83 @@
 }
 
+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,77 +1488,3 @@
  
 }
-
-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
+*********************************************************************/
