Changeset 3199 in Sophya for trunk/Cosmo/SimLSS


Ignore:
Timestamp:
Apr 5, 2007, 7:19:06 PM (18 years ago)
Author:
cmv
Message:

Intro raffinement dans simul AGN, variation dNu Da Dlum selon Oz cmv 05/04/2007

Location:
trunk/Cosmo/SimLSS
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvdefsurv.cc

    r3196 r3199  
    348348 // AGN
    349349 double flux_agn = pow(10.,lflux_agn);
    350  cout<<"AGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="<<flux_agn<<" Jy"<<endl;
     350 double mass_agn = FluxHI2Msol(flux_agn*Jansky2Watt_cst,dlum);
     351 cout<<"AGN: log10(S_agn)="<<lflux_agn<<" -> S_agn="
     352     <<flux_agn<<" Jy -> "<<mass_agn<<" equiv. Msol/Hz"<<endl;
    351353 double flux_agn_pix = flux_agn*(dnuhiz*1e9);
    352354 double mass_agn_pix = FluxHI2Msol(flux_agn_pix*Jansky2Watt_cst,dlum);

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3196 r3199  
    3131     <<" -2 : compute 2D spectrum"<<endl
    3232     <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
    33      <<" -A <log10(S_agn en Msol)>,sigma : distribution des AGN par pixel"<<endl
     33     <<" -A <log10(S_agn en Jy)>,sigma,powlaw : distribution des AGN par pixel"<<endl
     34     <<" -K : on utilise le calcul de spectre bricole pour les AGN (bidon!)"<<endl
    3435     <<" -W : write cube in FITS format"<<endl
    3536     <<" -P : write cube in PPF format"<<endl
    … …  
    7980 // *** AGN
    8081 bool do_agn = false;
    81  double lmsol_agn=-99., lsigma_agn=0.;   // en equivalent MSol
     82 double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
     83 double powlaw_agn = 0.;
     84 bool killkz = false;
    8285
    8386 // *** type de generation
    … …  
    9598
    9699 char c;
    97  while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:M:A:")) != -1) {
     100 while((c = getopt(narg,arg,"ha0PWV2GKx:y:z:s:Z:M:A:")) != -1) {
    98101  switch (c) {
    99102  case 'a' :
    … …  
    129132  case 'A' :
    130133    do_agn = true;
    131     sscanf(optarg,"%lf,%lf",&lmsol_agn,&lsigma_agn);
     134    sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
     135    break;
     136  case 'K' :
     137    killkz = true;
    132138    break;
    133139  case 'V' :
    … …  
    161167     <<", schnpt="<<schnpt<<endl;
    162168 cout<<"snoise="<<snoise<<" equivalent Msol"<<endl;
    163  if(do_agn) cout<<"AGN: <log10(Msol)>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
     169 if(do_agn)
     170   cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
     171       <<" , powlaw="<<powlaw_agn<<endl;
    164172
    165173 //-----------------------------------------------------------------
    … …  
    271279 fluct3d.SetCosmology(univ);
    272280 fluct3d.SetGrowthFactor(growth);
    273  fluct3d.LosComRedshift(0.001);
     281 fluct3d.LosComRedshift(0.001,-1);
    274282 TArray<r_8>& rgen = fluct3d.GetRealArray();
    275283 cout<<endl; fluct3d.Print();
    … …  
    383391   cout<<"\n--- Apply Growth factor"<<endl;
    384392   cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
    385    fluct3d.ApplyGrowthFactor(-1);
     393   fluct3d.ApplyGrowthFactor();
    386394   rmin,rmax; rgen.MinMax(rmin,rmax);
    387395   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
    … …  
    486494
    487495 if(do_agn) {
    488    cout<<"\n--- Add AGN: <mass>="<<lmsol_agn<<" , sigma="<<lsigma_agn<<endl;
    489    fluct3d.AddAGN(lmsol_agn,lsigma_agn);
     496   cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
     497       <<" , powlaw="<<powlaw_agn<<endl;
     498   fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
    490499   nm = fluct3d.MeanSigma2(rm,rs2);
    491500   cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
    … …  
    503512 }
    504513
     514 if(wfits) {
     515   fluct3d.WriteFits("!cmvobserv3d_rf.fits");
     516   PrtTim(">>>> End WriteFits");
     517 }
     518 if(wppf) {
     519   fluct3d.WritePPF("cmvobserv3d_rf.ppf",true);
     520   PrtTim(">>>> End WritePPF");
     521 }
     522
    505523 //-----------------------------------------------------------------
    506524 // -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
    … …  
    526544   hpkrec.ReCenterBin();
    527545   hpkrec.Show();
    528    fluct3d.ComputeSpectrum(hpkrec);
     546   if(killkz) fluct3d.ComputeSpectrum_bricolo(hpkrec);
     547     else     fluct3d.ComputeSpectrum(hpkrec);
    529548   tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
    530549   PrtTim(">>>> End Computing final spectrum");
    … …  
    536555   hpkrec2.ReCenterBin(); hpkrec2.Zero();
    537556   hpkrec2.Show();
    538    fluct3d.ComputeSpectrum2D(hpkrec2);
     557   if(killkz) fluct3d.ComputeSpectrum2D_bricolo(hpkrec2);
     558     else     fluct3d.ComputeSpectrum2D(hpkrec2);
    539559   {
    540560   tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
    … …  
    553573readfits cmvobserv3d_r0.fits
    554574readfits cmvobserv3d_r.fits
     575readfits cmvobserv3d_rf.fits
    555576
    556577openppf cmvobserv3d_k0.ppf
    … …  
    558579openppf cmvobserv3d_r0.ppf
    559580openppf cmvobserv3d_r.ppf
    560 
    561 # pour le plot 2D d'une slice en Z du 3D: to2d nom_obj3D num_slice
    562 defscript to2d
     581openppf cmvobserv3d_rf.ppf
     582
     583# pour le plot 2D d'une slice en Z du 3D: xy2d nom_obj3D num_slice
     584defscript xy2d
    563585 objaoper $1 sliceyz $2
    564  mv sliceyz_${2} ${1}_$2
    565  disp ${1}_$2
    566  echo display slice $2 of $1
     586 mv sliceyz_${2} ${1}_Z_$2
     587 disp ${1}_Z_$2
     588 echo display slice $2 of $1 name is ${1}_Z_$2
     589endscript
     590 
     591# pour le plot 2D d'une slice en Y du 3D: xz2d nom_obj3D num_slice
     592defscript xz2d
     593 objaoper $1 slicexy $2
     594 mv slicexy_${2} ${1}_Y_$2
     595 disp ${1}_Y_$2
     596 echo display slice $2 of $1 name is ${1}_Y_$2
    567597endscript
    568 
    569 to2d $cobj 0
     598 
     599# pour le plot 2D d'une slice en X du 3D: yz2d nom_obj3D num_slice
     600defscript yz2d
     601 objaoper $1 slicexz $2
     602 mv slicexz_${2} ${1}_X_$2
     603 disp ${1}_X_$2
     604 echo display slice $2 of $1 name is ${1}_X_$2
     605endscript
     606
     607xy2d $cobj 0
     608xz2d $cobj 0
     609yz2d $cobj 0
    570610
    571611######################################################
    … …  
    575615set k pow(10.,x)
    576616n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
     617
     618echo ${hpkgen.sum}
     619echo ${hpkgenf.sum}
     620echo ${hpkrec.sum}
    577621
    578622zone

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3196 r3199  
    2121#include "constcosmo.h"
    2222#include "geneutils.h"
     23#include "schechter.h"
    2324
    2425#include "genefluct3d.h"
    … …  
    3334  , redshref_(-999.) , kredshref_(0.) , cosmo_(NULL) , growth_(NULL)
    3435  , loscom_ref_(-999.), loscom_min_(-999.), loscom_max_(-999.)
    35 
    36 
     36  , loscom2zred_min_(0.) , loscom2zred_max_(0.)
    3737{
    3838 xobs_[0] = xobs_[1] = xobs_[2] = 0.;
    3939 zred_.resize(0);
    4040 loscom_.resize(0);
     41 loscom2zred_.resize(0);
    4142 SetNThread();
    4243}
    … …  
    7677 redshref_  = redshref;
    7778 kredshref_ = kredshref;
     79 if(lp_>0)
     80   cout<<"--- GeneFluct3D::SetObservator zref="<<redshref_<<" kref="<<kredshref_<<endl;
    7881}
    7982
    … …  
    9497                <<" D="<<dx<<","<<dy<<","<<dz<<endl;
    9598 if(nx<=0 || dx<=0.) {
    96    cout<<"GeneFluct3D::setsize_Error: bad value(s)"<<endl;
    97    throw ParmError("GeneFluct3D::setsize_Error: bad value(s)");
     99   char *bla = "GeneFluct3D::setsize_Error: bad value(s";
     100   cout<<bla<<endl; throw ParmError(bla);
    98101 }
    99102
    … …  
    161164 char bla[90];
    162165 sprintf(bla,"GeneFluct3D::check_array_alloc_Error: array is not allocated");
    163  cout<<bla<<endl;
    164  throw ParmError(bla);
     166 cout<<bla<<endl; throw ParmError(bla);
    165167}
    166168
    … …  
    183185
    184186//-------------------------------------------------------
    185 long GeneFluct3D::LosComRedshift(double zinc)
     187long GeneFluct3D::LosComRedshift(double zinc,long npoints)
    186188// Given a position of the cube relative to the observer
    187189// and a cosmology
    … …  
    190192//   the vector "zred_" of scanned redshift (by zinc increments)
    191193//   the vector "loscom_" of corresponding los comoving distance
     194// -- Input:
     195// zinc : redshift increment for computation
     196// npoints : number of points required for inverting loscom -> zred
    192197//
    193198{
    194  if(zinc<=0.) zinc = 0.01;
    195  if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<endl;
     199 if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
    196200
    197201 if(cosmo_ == NULL || redshref_<0.) {
    198    cout<<"GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first"<<endl;
    199    throw ParmError("");
    200  }
    201 
    202  // On calcule les coordonnees de l'observateur
    203  // Il est sur un axe centre sur le milieu de la face Oxy
     202   char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
     203   cout<<bla<<endl; throw ParmError(bla);
     204 }
     205
     206 // On calcule les coordonnees de l'observateur dans le repere du cube
     207 // cad dans le repere ou l'origine est au centre du pixel i=j=l=0.
     208 // L'observateur est sur un axe centre sur le milieu de la face Oxy
    204209 double loscom_ref_ = cosmo_->Dloscom(redshref_);
    205210 xobs_[0] = Nx_/2.*Dx_;
    … …  
    239244 }
    240245
    241  // Fill the corresponding vectors
     246 // Fill the corresponding vectors for loscom and zred
     247 if(zinc<=0.) zinc = 0.01;
    242248 for(double z=0.; ; z+=zinc) {
    243249   double dlc = cosmo_->Dloscom(z);
    … …  
    246252   loscom_.push_back(dlc);
    247253   z += zinc;
    248    if(dlc>loscom_max_) break; // on break apres avoir stoque un dlc>dlcmax
    249  }
    250 
    251  long n = zred_.size();
     254   if(dlc>loscom_max_) break; // on sort apres avoir stoque un dlc>dlcmax
     255 }
     256
    252257 if(lp_>0) {
    253    cout<<"...n="<<n;
     258   long n = zred_.size();
     259   cout<<"...zred/loscom tables[zinc="<<zinc<<"]: n="<<n;
    254260   if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
    255261   if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
    … …  
    257263 }
    258264
    259  return n;
     265 // Compute the parameters and tables needed for inversion loscom->zred
     266 if(npoints<3) npoints = zred_.size();
     267 InverseFunc invfun(zred_,loscom_);
     268 invfun.ComputeParab(npoints,loscom2zred_);
     269 loscom2zred_min_ = invfun.YMin();
     270 loscom2zred_max_ = invfun.YMax();
     271
     272 if(lp_>0) {
     273   long n = loscom2zred_.size();
     274   cout<<"...loscom -> zred[npoints="<<npoints<<"]: n="<<n
     275       <<" los_min="<<loscom2zred_min_
     276       <<" los_max="<<loscom2zred_max_
     277       <<" -> zred=[";
     278   if(n>0) cout<<loscom2zred_[0];
     279   cout<<",";
     280   if(n>1) cout<<loscom2zred_[n-1];
     281   cout<<"]"<<endl;
     282   if(lp_>1 && n>0)
     283     for(int i=0;i<n;i++)
     284       if(i==0 || abs(i-n/2)<2 || i==n-1)
     285         cout<<"    "<<i<<"  "<<loscom2zred_[i]<<endl;
     286 }
     287
     288 return zred_.size();
    260289}
    261290
    … …  
    624653
    625654//-------------------------------------------------------------------
    626 void GeneFluct3D::ApplyGrowthFactor(long npoints)
     655void GeneFluct3D::ApplyGrowthFactor(void)
    627656// Apply Growth to real space
    628657// Using the correspondance between redshift and los comoving distance
    629658// describe in vector "zred_" "loscom_"
    630659{
    631  if(npoints<3) npoints = zred_.size();
    632  if(lp_>0) cout<<"--- ApplyGrowthFactor ---  npoints="<<npoints<<endl;
     660 if(lp_>0) cout<<"--- ApplyGrowthFactor ---"<<endl;
    633661 check_array_alloc();
    634662
    635663 if(growth_ == NULL) {
    636    cout<<"GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first"<<endl;
    637    throw ParmError("GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first");
    638  }
    639 
    640  long n = zred_.size();
    641  InverseFunc invfun(zred_,loscom_);
    642  vector<double> invlos;
    643  invfun.ComputeParab(npoints,invlos);
    644 
    645  InterpFunc interpinv(invfun.YMin(),invfun.YMax(),invlos);
     664   char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
     665   cout<<bla<<endl; throw ParmError(bla);
     666 }
     667
     668 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
    646669 unsigned short ok;
    647670
    … …  
    969992}
    970993
    971 void GeneFluct3D::AddAGN(double lmsol,double lsigma)
     994void GeneFluct3D::AddAGN(double lfjy,double lsigma,double powlaw)
    972995// Add AGN flux into simulation:
    973996// --- Procedure:
    974997// 1. lancer "cmvdefsurv" avec les parametres du survey
     998//        (au redshift de reference du survey)
    975999//    et recuperer l'angle solide "angsol sr" du pixel elementaire
    9761000//    au centre du cube.
    … …  
    9781002// 3. regarder l'histo "hlfang" et en deduire un equivalent gaussienne
    9791003//    cad une moyenne <log10(S)> et un sigma "sig"
    980 //    Attention: la distribution n'etant pas gaussienne
    981 // 4. re-lancer "cmvdefsurv" en ajoutant l'info sur les AGN
    982 //    "cmvdefsurv ... -G <log10(S)> ..."
    983 //    et recuperer le log10(masse solaire equivalente)
     1004//    Attention: la distribution n'est pas gaussienne les "mean,sigma"
     1005//               de l'histo ne sont pas vraiment ce que l'on veut
    9841006// --- Limitations actuelle du code:
    985 // a. l'angle solide du pixel est pris au centre du cube
    986 //    et ne varie pas avec la distance a l'interieur du cube
    987 // b. la taille en dNu des pixels (selon z) est supposee constante
    988 //    et egale a celle calculee pour le centre du cube
    989 // c. les AGN sont supposes plats en flux
    990 // d. le flux des AGN est mis dans une colonne Oz (indice k)
    991 //    et pas sur la ligne de visee
    992 // e. la distribution est approximee a une gaussienne
    993 // .. C'est une approximation pour un observateur loin du centre du cube
    994 //    et pour un cube peu epais / distance observateur
     1007// . les AGN sont supposes en loi de puissance IDENTIQUE pour tout theta,phi
     1008// . le flux des AGN est mis dans une colonne Oz (indice k) et pas sur la ligne de visee
     1009// . la distribution est approximee a une gaussienne
     1010// ... C'est une approximation pour un observateur loin du centre du cube
     1011//     et pour un cube peu epais / distance observateur
    9951012// --- Parametres de la routine:
    996 // lmsol : c'est le <log10(Msol)> qui est la conversion en masse solaire
    997 //         du flux depose dans un pixel par les AGN
     1013// llfy : c'est le <log10(S)> du flux depose dans un pixel par les AGN
    9981014// lsigma : c'est le sigma de la distribution
     1015// powlaw : c'est la pente de ls distribution cad que le flux "lmsol"
     1016//          et considere comme le flux a 1.4GHz et qu'on suppose une loi
     1017//             F(nu) = (1.4GHz/nu)^powlaw * F(1.4GHz)
    9991018// - Comme on est en echelle log10():
    10001019//   on tire log10(Msol) + X
    … …  
    10031022// - Pas de probleme de pixel negatif car on a une multiplication!
    10041023{
    1005   if(lp_>0) cout<<"--- AddAGN: lmsol = "<<lmsol<<" , sigma = "<<lsigma<<endl;
     1024  if(lp_>0) cout<<"--- AddAGN: <log10(S Jy)> = "<<lfjy<<" , sigma = "<<lsigma<<endl;
    10061025  check_array_alloc();
    10071026
    1008   double m = pow(10.,lmsol);
    1009 
    1010   double sum=0., sum2=0., nsum=0.;
    1011   for(long l=0;l<Nz_;l++) {
    1012     double a = lsigma*NorRand();
    1013     a = m*pow(10.,a);
    1014     // On met le meme tirage le long de Oz (indice k)
    1015     for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
    1016       int_8 ip = IndexR(i,j,l);
    1017       data_[ip] += a;
    1018     }
    1019     sum += a; sum2 += a*a; nsum += 1.;
    1020   }
    1021 
    1022  if(nsum>1.) {
     1027 if(cosmo_ == NULL || redshref_<0.| loscom2zred_.size()<1) {
     1028   char *bla = "GeneFluct3D::AddAGN_Error: set Observator and Cosmology first";
     1029   cout<<bla<<endl; throw ParmError(bla);
     1030 }
     1031
     1032 // La distance angulaire/luminosite/Dnu au centre
     1033 double dangref = cosmo_->Dang(redshref_);
     1034 double dlumref = cosmo_->Dlum(redshref_);
     1035 double dredref  = Dz_/(cosmo_->Dhubble()/cosmo_->E(redshref_));
     1036 double dnuref = Fr_HyperFin_Par *dredref/pow(1.+redshref_,2.); // GHz
     1037 double fagnref = pow(10.,lfjy)*(dnuref*1.e9); // Jy.Hz
     1038 double magnref = FluxHI2Msol(fagnref*Jansky2Watt_cst,dlumref); // Msol
     1039 if(lp_>0) {
     1040   cout<<"Au centre: z="<<redshref_<<", dredref="<<dredref<<", dnuref="<<dnuref<<" GHz"<<endl
     1041       <<" dang="<<dangref<<" Mpc, dlum="<<dlumref<<" Mpc,"
     1042       <<" fagnref="<<fagnref<<" Jy.Hz (a 1.4GHz), magnref="<<magnref<<" Msol"<<endl;
     1043 }
     1044
     1045 if(powlaw!=0.) {
     1046   // F(nu) = (nu GHz/1.4 Ghz)^p * F(1.4GHz) et nu = 1.4 GHz / (1+z)
     1047   magnref *= pow(1/(1.+redshref_),powlaw);
     1048   if(lp_>0) cout<<" powlaw="<<powlaw<<"  -> change magnref to "<<magnref<<" Msol"<<endl;
     1049 }
     1050
     1051 // Les infos en fonction de l'indice "l" selon Oz
     1052 vector<double> correction;
     1053 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
     1054 for(long l=0;l<Nz_;l++) {
     1055   double z = fabs(xobs_[2] - l*Dz_);
     1056   double zred = interpinv(z);
     1057   double dang = cosmo_->Dang(zred);  // pour variation angle solide
     1058   double dlum = cosmo_->Dlum(zred);  // pour variation conversion mass HI
     1059   double dred = Dz_/(cosmo_->Dhubble()/cosmo_->E(zred));
     1060   double dnu  = Fr_HyperFin_Par *dred/pow(1.+zred,2.); // pour variation dNu
     1061   double corr = dnu/dnuref*pow(dangref/dang*dlum/dlumref,2.);
     1062   if(powlaw!=0.) corr *= pow((1.+redshref_)/(1.+zred),powlaw);
     1063   correction.push_back(corr);
     1064   if(lp_>0 && (l==0 || abs(l-Nz_/2)<2 || l==Nz_-1)) {
     1065     cout<<"l="<<l<<" z="<<z<<" red="<<zred
     1066         <<" da="<<dang<<" dlu="<<dlum<<" dred="<<dred
     1067         <<" dnu="<<dnu<<" -> corr="<<corr<<endl;
     1068   }
     1069 }
     1070
     1071 double sum=0., sum2=0., nsum=0.;
     1072 for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) {
     1073   double a = lsigma*NorRand();
     1074   a = magnref*pow(10.,a);
     1075   // On met le meme tirage le long de Oz (indice k)
     1076   for(long l=0;l<Nz_;l++) {
     1077     int_8 ip = IndexR(i,j,l);
     1078     data_[ip] += a*correction[l];
     1079   }
     1080   sum += a; sum2 += a*a; nsum += 1.;
     1081 }
     1082
     1083 if(lp_>0 && nsum>1.) {
    10231084   sum /= nsum;
    10241085   sum2 = sum2/nsum - sum*sum;
    … …  
    10391100 }
    10401101}
     1102
     1103
     1104
     1105//-------------------------------------------------------------------
     1106//-------------------------------------------------------------------
     1107//--------------------- BRICOLO A NE PAS GARDER ---------------------
     1108//-------------------------------------------------------------------
     1109//-------------------------------------------------------------------
     1110int GeneFluct3D::ComputeSpectrum_bricolo(HistoErr& herr)
     1111// Compute spectrum from "T" and fill HistoErr "herr"
     1112// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
     1113// cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
     1114{
     1115 if(lp_>0) cout<<"--- ComputeSpectrum_bricolo ---"<<endl;
     1116 check_array_alloc();
     1117
     1118 if(herr.NBins()<0) return -1;
     1119 herr.Zero();
     1120
     1121 // Attention a l'ordre
     1122 for(long i=0;i<Nx_;i++) {
     1123   long ii = (i>Nx_/2) ? Nx_-i : i;
     1124   double kx = ii*Dkx_;  kx *= kx;
     1125   for(long j=0;j<Ny_;j++) {
     1126     if(i==0 && j==0) continue;
     1127     long jj = (j>Ny_/2) ? Ny_-j : j;
     1128     double ky = jj*Dky_;  ky *= ky;
     1129     for(long l=1;l<NCz_;l++) {
     1130       double kz = l*Dkz_;  kz *= kz;
     1131       double k = sqrt(kx+ky+kz);
     1132       double pk = MODULE2(T_(l,j,i));
     1133       herr.Add(k,pk);
     1134     }
     1135   }
     1136 }
     1137 herr.ToVariance();
     1138
     1139 // renormalize to directly compare to original spectrum
     1140 double norm = Vol_;
     1141 herr *= norm;
     1142
     1143 return 0;
     1144}
     1145
     1146int GeneFluct3D::ComputeSpectrum2D_bricolo(Histo2DErr& herr)
     1147{
     1148 if(lp_>0) cout<<"--- ComputeSpectrum2D_bricolo ---"<<endl;
     1149 check_array_alloc();
     1150
     1151 if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
     1152 herr.Zero();
     1153
     1154 // Attention a l'ordre
     1155 for(long i=0;i<Nx_;i++) {
     1156   long ii = (i>Nx_/2) ? Nx_-i : i;
     1157   double kx = ii*Dkx_;  kx *= kx;
     1158   for(long j=0;j<Ny_;j++) {
     1159     if(i==0 && j==0) continue;
     1160     long jj = (j>Ny_/2) ? Ny_-j : j;
     1161     double ky = jj*Dky_;  ky *= ky;
     1162     double kt = sqrt(kx+ky);
     1163     for(long l=1;l<NCz_;l++) {
     1164       double kz = l*Dkz_;
     1165       double pk = MODULE2(T_(l,j,i));
     1166       herr.Add(kt,kz,pk);
     1167     }
     1168   }
     1169 }
     1170 herr.ToVariance();
     1171
     1172 // renormalize to directly compare to original spectrum
     1173 double norm = Vol_;
     1174 herr *= norm;
     1175
     1176 return 0;
     1177}

  • trunk/Cosmo/SimLSS/genefluct3d.h

    r3196 r3199  
    3131  void SetCosmology(CosmoCalc& cosmo);
    3232  void SetGrowthFactor(GrowthFactor& growth);
    33   long LosComRedshift(double zinc=0.001);
     33  long LosComRedshift(double zinc=0.001,long npoints=-1);
    3434
    3535  TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
    … …  
    7171  int  ComputeSpectrum(HistoErr& herr);
    7272  int  ComputeSpectrum2D(Histo2DErr& herr);
    73   void ApplyGrowthFactor(long npoints=-1);
     73  void ApplyGrowthFactor(void);
    7474
    7575  int_8 VarianceFrReal(double R,double& var);
    … …  
    8383  double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
    8484  double TurnMass2Flux(void);
    85   void AddAGN(double lmsol,double lsigma);
     85  void AddAGN(double lfjy,double lsigma,double powlaw=0.);
    8686  void AddNoise2Real(double snoise);
    8787
    … …  
    9494  void SetPrtLevel(int lp=0) {lp_ = lp;}
    9595  void Print(void);
     96
     97//-------------------------------------------------------------------
     98//--------------------- BRICOLO A NE PAS GARDER ---------------------
     99//-------------------------------------------------------------------
     100  int  ComputeSpectrum_bricolo(HistoErr& herr);
     101  int  ComputeSpectrum2D_bricolo(Histo2DErr& herr);
     102//-------------------------------------------------------------------
    96103
    97104protected:
    … …  
    138145  double loscom_ref_, loscom_min_, loscom_max_;
    139146  vector<double> zred_, loscom_;
     147  double loscom2zred_min_, loscom2zred_max_;
     148  vector<double> loscom2zred_;
     149
    140150};
    141151

  • trunk/Cosmo/SimLSS/geneutils.cc

    r3196 r3199  
    9292
    9393int InverseFunc::ComputeLinear(long n,vector<double>& xfcty)
     94// Compute table "xfcty" by linear interpolation of "x" versus "y"
     95//   on "n" points from "ymin" to "ymax":
     96// xfcty[i] = interpolation of function "x" for "ymin+i*(ymax-ymin)/(n-1.)"
    9497{
    9598  if(n<3) return -1;
Note: See TracChangeset for help on using the changeset viewer.