Changeset 3768 in Sophya for trunk/Cosmo/SimLSS


Ignore:
Timestamp:
May 3, 2010, 4:08:34 PM (15 years ago)
Author:
cmv
Message:
  • refonte du code pour creer uniquement des conditions initiales
  • introduction du tirage des vitesse LOS pour les redshift-distortion

cmv 03/05/2010

Location:
trunk/Cosmo/SimLSS
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/Makefile

    r3746 r3768  
    1111
    1212MYEXTINC = ${EXTLIBDIR}/${MACH}/include
    13 MYLIB  = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lcmvgenfluc -lfftw3_threads -lfftw3 -lm
    14 MYLIB4 = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lcmvgenfluc4 -lfftw3f_threads -lfftw3f -lm
     13MYLIB  = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvgenfluc -lcmvsimbao -lfftw3_threads -lfftw3 -lm
     14MYLIB4 = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvgenfluc4 -lcmvsimbao -lfftw3f_threads -lfftw3f -lm
    1515
    1616#--------------------------------------------------------------------------
    … …  
    1818PROGS = \
    1919       $(EXE)cmvobserv3d $(EXE)cmvobserv3df \
     20       $(EXE)cmvginit3d $(EXE)cmvginit3df \
    2021       $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
    2122       $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
    … …  
    2526PROGSOBJ = \
    2627          $(OBJ)cmvobserv3d.o $(OBJ)cmvobserv3df.o \
     28          $(OBJ)cmvginit3d.o $(OBJ)cmvginit3df.o \
    2729          $(OBJ)cmvtuniv.o $(OBJ)cmvtransf.o $(OBJ)cmvtgrowth.o $(OBJ)cmvtstpk.o \
    2830          $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
    … …  
    189191
    190192##############################################################################
     193cmvginit3d: $(EXE)cmvginit3d
     194        echo $@ " done"
     195$(EXE)cmvginit3d: $(OBJ)cmvginit3d.o $(LIBR) $(LIBG)
     196        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvginit3d.o $(MYLIB)
     197$(OBJ)cmvginit3d.o: cmvginit3d.cc
     198        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvginit3d.cc
     199
     200cmvginit3df: $(EXE)cmvginit3df
     201        echo $@ " done"
     202$(EXE)cmvginit3df: $(OBJ)cmvginit3df.o $(LIBR) $(LIBG4)
     203        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvginit3df.o $(MYLIB4)
     204$(OBJ)cmvginit3df.o: cmvginit3d.cc
     205        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -DGEN3D_FLOAT -o $@ cmvginit3d.cc
     206
     207##############################################################################
    191208cmvtintfun: $(EXE)cmvtintfun
    192209        echo $@ " done"

  • trunk/Cosmo/SimLSS/cmvhshsns.cc

    r3653 r3768  
    33// > cmvhshsns -n -d 0.105 -g 2,0.105,0. -i 2,0.41,0. -t 50,180,0 -f 1420 1410 1430
    44// > cmvhshsns -n -D -1 -d 10 -g 1 -i 10,10,0. -t 50,180,0 -f 1420 1410 1430
     5// - Pittsburgh Sun 24 Nov
     6//   > cmvhshsns -n -d 0.147 -D 0 -g 1 -i 8,0.147,61.06 -t 360,90,0 -f 1420
    57#include "sopnamsp.h"
    68#include "machdefs.h"
    … …  
    2830  <<"  -D : +1 lobe du dipole approximation de Hertz"<<endl
    2931  <<"       -1 le lobe du dipole est remplace par le sinc() du cylindre (etude E-W)"<<endl
     32  <<"        0 lobe du dipole a 2 brins de longueur \"L\""<<endl
    3033  <<"  -g N_g,D_g,Theta_g : regroupement de dipoles"<<endl
    3134  <<"                       N_g<=1 pas de regroupement"<<endl

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3746 r3768  
    1717#include "constcosmo.h"
    1818#include "cosmocalc.h"
    19 #include "schechter.h"
    2019#include "geneutils.h"
    2120#include "genefluct3d.h"
    … …  
    141140 string rootnameout = "cmvobserv3d";
    142141
     142 // *** Number of entries in control NTuple
    143143 unsigned long ntnent = 10000;  // 0 = do not fill NTuple
    144144
    … …  
    415415
    416416 GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
    417  fluct3d.SetObservator(zref,-nz/2.);
    418417 fluct3d.SetCosmology(univ);
    419418 fluct3d.SetGrowthFactor(growth);
     419 fluct3d.SetObservator(zref,-nz/2.);
    420420 fluct3d.LosComRedshift(0.001,-1);
    421421 //TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();

  • trunk/Cosmo/SimLSS/cmvtgrowth.cc

    r3572 r3768  
    1414
    1515void usage(void);
    16 void usage(void) {cout<<"cmvtgrowth z1,z2,dz [Omatter0,Lambda0]"<<endl;}
    17 
    18 void tstprint(CosmoCalc& univ,double z1,double z2,double dz);
    19 void tstspeed(CosmoCalc& univ,double z1,double z2,double dz);
    20 void tstntuple(CosmoCalc& univ,double z1,double z2,double dz);
     16void usage(void) {cout<<"cmvtgrowth z1,z2,nz [Omatter0,Lambda0]"<<endl;}
    2117
    2218int main(int narg,char *arg[])
    2319{
    2420
    25  double z1=0., z2=10., dz=1.;
    26  if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);
    27  cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
     21 double z1=0., z2=10.;
     22 int nz = 100;;
     23 if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&z1,&z2,&nz);
     24 if(nz<=0) nz = 100;
     25 double dz = (z2-z1)/nz;
     26 cout<<"z1="<<z1<<"  z2="<<z2<<"  nz="<<nz<<"  dz="<<dz<<endl;
    2827
    2928 double om0=0.267804, ol0=0.73;
    … …  
    3433 cout<<"D1(z=0) = "<<growth(0.)<<endl;
    3534
    36  const int n = 2;
    37  const char *vname[n] = {"z","d1"};
     35 const int n = 3;
     36 const char *vname[n] = {"z","d1","d1dz"};
    3837 NTuple nt(n,vname);
    3938 double xnt[n];
    … …  
    4140   xnt[0]  = z;
    4241   xnt[1]  = growth(z);
     42   xnt[2]  = growth.DsDz(z,dz/2.);
    4343   nt.Fill(xnt);
    4444 }
    … …  
    5959set cut z<5
    6060
     61# --- growth
    6162n/plot nt.d1%z $cut ! "nsta connectpoints"
    6263n/plot nt.1./(1.+z)%z $cut ! "nsta connectpoints red same"
     64
     65# --- growth'/growth
     66n/plot nt.d1dz/d1%z $cut ! "nsta connectpoints"
     67n/plot nt.-1./(1.+z)%z $cut ! "nsta connectpoints red same"
     68
     69# --- d(growth)/dz
     70zmin = 0.
     71zmax = 10.
     72set npt ${dv.size}
     73dd = ($zmax-$zmin)/(${npt}-1.)
     74exptovec dv nt d1
     75c++exec for(int i=1;i<dv.Size(); i++) dv(i-1) = dv(i)-dv(i-1);
     76
     77n/plot nt.-1./pow(1.+z,2.)%z $cut ! "nsta connectpoints red"
     78n/plot nt.d1dz%z $cut ! "nsta connectpoints same"
     79n/plot dv.val/${dd}%${zmin}+(n+1)*$dd ! ! "nsta plusmarker5 green same"
     80
    6381*/

  • trunk/Cosmo/SimLSS/cosmocalc.cc

    r3749 r3768  
    2121
    2222//----------------------------------------------------------
    23 // flat = 0 : no compute Otot from other
     23// flat = 0 : do not set Otot=1, compute Otot from others
    2424//      = 1 : change Omatter to get Otot=1
    2525//      = 2 : change Olambda to get Otot=1
    … …  
    9696 if(_xspl) delete [] _xspl; _xspl = NULL;
    9797 if(_yspl) delete [] _yspl; _yspl = NULL;
     98}
     99
     100//----------------------------------------------------------
     101// On change le zmax de la cosmologie
     102void CosmoCalc::ChangeZmax(double zmax,double dzinc,double dzmax)
     103{
     104  _zmax = zmax;
     105  if(dzinc<=0. ) dzinc = _dzinc;
     106  if(dzmax<=0.)  dzmax = _dzmax;
     107  cout<<"CosmoCalc::ChangeZmax: zmax="<<_zmax<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<endl;
     108  SetInteg(_dperc,dzinc,_dzmax,_glorder);
    98109}
    99110

  • trunk/Cosmo/SimLSS/cosmocalc.h

    r3749 r3768  
    2121
    2222  double ZMax(void) { return _zmax;}
     23  void ChangeZmax(double zmax,double dzinc=-1.,double dzmax=-1.);
    2324  void SetInteg(double dperc=-1.,double dzinc=-1.,double dzmax=-1.,unsigned short order=4);
    2425  void PrtInteg(void);

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3746 r3768  
    1 #include "sopnamsp.h"
    21#include "machdefs.h"
    32#include <iostream>
    … …  
    8887 cosmo_ = NULL;
    8988 growth_ = NULL;
     89 good_dzinc_ = 0.01;
     90 compute_pk_redsh_ref_ = -999.;
    9091 redsh_ref_ = -999.;
    9192 kredsh_ref_ = 0.;
    … …  
    262263//
    263264{
     265 if(zinc<=0.) zinc = good_dzinc_;
    264266 if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<" , npoints="<<npoints<<endl;
    265267
    266  if(cosmo_ == NULL || redsh_ref_<0.) {
    267    const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first";
     268 if(cosmo_ == NULL || growth_==NULL || redsh_ref_<0.) {
     269   const char *bla = "GeneFluct3D::LosComRedshift_Error: set Observator, Cosmology and Growth first";
    268270   cout<<bla<<endl; throw ParmError(bla);
    269271 }
    … …  
    275277 dlum_ref_ = cosmo_->Dlum(redsh_ref_);
    276278 dang_ref_ = cosmo_->Dang(redsh_ref_);
     279 h_ref = cosmo_->H(redsh_ref_);
     280 growth_ref_ = (*growth_)(redsh_ref_);
     281 dsdz_growth_ref_ = growth_->DsDz(redsh_ref_,zinc);
    277282 nu_ref_   = Fr_HyperFin_Par/(1.+redsh_ref_); // GHz
    278283 dnu_ref_  = Fr_HyperFin_Par *dred_ref_/pow(1.+redsh_ref_,2.); // GHz
    … …  
    282287       <<", nuref="<<nu_ref_ <<" GHz"
    283288       <<", dnuref="<<dnu_ref_ <<" GHz"<<endl
     289       <<"   H="<<h_ref<<" km/s/Mpc"
     290       <<", D="<<growth_ref_
     291       <<", dD/dz="<<dsdz_growth_ref_<<endl
    284292       <<"   dlosc="<<loscom_ref_<<" Mpc com"
    285293       <<", dtrc="<<dtrc_ref_<<" Mpc com"
    … …  
    335343
    336344 // Fill the corresponding vectors for loscom and zred
    337  // Be shure to have one dlc <loscom_min and one >loscom_max
    338  if(zinc<=0.) zinc = 0.01;
     345 // Be shure to have one dlc<loscom_min and one dlc>loscom_max
    339346 double zmin = 0., dlcmin=0.;
    340347 while(1) {
    … …  
    397404 }
    398405
     406
     407 // Compute the table for D'(z)/D(z) where D'(z)=dD/dEta (Eta is conformal time)
     408 zredmin_dpsd_ = zred_[0];
     409 zredmax_dpsd_ = zred_[zred_.size()-1];
     410 long nptd = long(sqrt(Nx_*Nx_ + Ny_*Ny_ +Nz_*Nz_));
     411 if(nptd<10) nptd = 10;
     412 good_dzinc_ = (zredmax_dpsd_ - zredmin_dpsd_)/nptd;
     413 if(lp_>0) cout<<"...good_dzinc changed to "<<good_dzinc_<<endl;
     414 dpsdfrzred_.resize(0);
     415 double dz = (zredmax_dpsd_ - zredmin_dpsd_)/(nptd-1);
     416 int nmod = nptd/5; if(nmod==0) nmod = 1;
     417 if(lp_>0) cout<<"...Compute table D'/D on "<<nptd<<" pts for z["
     418               <<zredmin_dpsd_<<","<<zredmax_dpsd_<<"]"<<endl;
     419 for(int i=0;i<nptd;i++) {
     420   double z = zredmin_dpsd_ + i*dz;
     421   double v = -cosmo_->H(z) * growth_->DsDz(z,good_dzinc_) / (*growth_)(z);
     422   dpsdfrzred_.push_back(v);
     423   if(lp_ && (i%nmod==0 || i==nptd-1)) cout<<"    z="<<z<<"  D'/D="<<v<<" km/s/Mpc"<<endl;
     424 }
     425
    399426 return zred_.size();
    400427}
    … …  
    415442   fwrt.WriteKey("DKY",Dky_," Mpc^-1");
    416443   fwrt.WriteKey("DKZ",Dkz_," Mpc^-1");
     444   fwrt.WriteKey("ZREFPK",compute_pk_redsh_ref_," Pk computed redshift");
    417445   fwrt.WriteKey("ZREF",redsh_ref_," reference redshift");
    418446   fwrt.WriteKey("KZREF",kredsh_ref_," reference redshift on z axe");
     447   fwrt.WriteKey("Growth",Dref()," growth at reference redshift");
     448   fwrt.WriteKey("GrowthPk",DrefPk()," growth at Pk computed redshift");
    419449   fwrt.Write(R_);
    420450 } catch (PThrowable & exc) {
    … …  
    440470   double dy = fimg.ReadKey("DY");
    441471   double dz = fimg.ReadKey("DZ");
     472   double pkzref = fimg.ReadKey("ZREFPK");
    442473   double zref = fimg.ReadKey("ZREF");
    443474   double kzref = fimg.ReadKey("KZREF");
    … …  
    447478   SetObservator(zref,kzref);
    448479   array_allocated_ = true;
     480   compute_pk_redsh_ref_ = pkzref;
    449481 } catch (PThrowable & exc) {
    450482   cout<<"Exception : "<<(string)typeid(exc).name()
    … …  
    468500   R_.Info()["DY"] = (r_8)Dy_;
    469501   R_.Info()["DZ"] = (r_8)Dz_;
     502   R_.Info()["ZREFPK"] = (r_8)compute_pk_redsh_ref_;
    470503   R_.Info()["ZREF"] = (r_8)redsh_ref_;
    471504   R_.Info()["KZREF"] = (r_8)kredsh_ref_;
     505   R_.Info()["Growth"] = (r_8)Dref();
     506   R_.Info()["GrowthPk"] = (r_8)DrefPk();
    472507   POutPersist pos(cfname.c_str());
    473508   if(write_real) pos << PPFNameTag("rgen")  << R_;
    … …  
    506541   r_8 dy = R_.Info()["DY"];
    507542   r_8 dz = R_.Info()["DZ"];
     543   r_8 pkzref = R_.Info()["ZREFPK"];
    508544   r_8 zref = R_.Info()["ZREF"];
    509545   r_8 kzref = R_.Info()["KZREF"];
    … …  
    512548   SetObservator(zref,kzref);
    513549   array_allocated_ = true;
     550   compute_pk_redsh_ref_ = pkzref;
    514551 } catch (PThrowable & exc) {
    515552   cout<<"Exception : "<<(string)typeid(exc).name()
    … …  
    630667
    631668//-------------------------------------------------------
    632 void GeneFluct3D::ComputeFourier0(GenericFunc& pk_at_z)
     669void GeneFluct3D::ComputeFourier0(PkSpectrumZ& pk_at_z)
    633670// cf ComputeFourier() mais avec autre methode de realisation du spectre
    634671//    (attention on fait une fft pour realiser le spectre)
    635672{
    636 
     673  compute_pk_redsh_ref_ = pk_at_z.GetZ();
    637674 // --- realisation d'un tableau de tirage gaussiens
    638  if(lp_>0) cout<<"--- ComputeFourier0: before gaussian filling ---"<<endl;
     675  if(lp_>0) cout<<"--- ComputeFourier0 at z="<<compute_pk_redsh_ref_
     676                <<": before gaussian filling ---"<<endl;
    639677 // On tient compte du pb de normalisation de FFTW3
    640678 double sntot = sqrt((double)NRtot_);
    … …  
    651689 if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
    652690 T_(0,0,0) = complex<GEN3D_TYPE>(0.);  // on coupe le continue et on l'initialise
    653  long lmod = Nx_/10; if(lmod<1) lmod=1;
     691 long lmod = Nx_/20; if(lmod<1) lmod=1;
    654692 for(long i=0;i<Nx_;i++) {
    655    long ii = (i>Nx_/2) ? Nx_-i : i;
    656    double kx = ii*Dkx_;  kx *= kx;
    657    if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
     693   double kx = Kx(i);  kx *= kx;
     694   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
    658695   for(long j=0;j<Ny_;j++) {
    659      long jj = (j>Ny_/2) ? Ny_-j : j;
    660      double ky = jj*Dky_;  ky *= ky;
     696     double ky = Ky(j);  ky *= ky;
    661697     for(long l=0;l<NCz_;l++) {
    662        double kz = l*Dkz_;  kz *= kz;
     698       double kz = Kz(l);  kz *= kz;
    663699       if(i==0 && j==0 && l==0) continue; // Suppression du continu
    664700       double k = sqrt(kx+ky+kz);
    … …  
    680716
    681717//-------------------------------------------------------
    682 void GeneFluct3D::ComputeFourier(GenericFunc& pk_at_z)
     718void GeneFluct3D::ComputeFourier(PkSpectrumZ& pk_at_z)
    683719// Calcule une realisation du spectre "pk_at_z"
    684720// Attention: dans TArray le premier indice varie le + vite
    … …  
    692728 // --- RaZ du tableau
    693729 T_ = complex<GEN3D_TYPE>(0.);
     730 compute_pk_redsh_ref_ = pk_at_z.GetZ();
    694731
    695732 // --- On remplit avec une realisation
    696  if(lp_>0) cout<<"--- ComputeFourier ---"<<endl;
    697  long lmod = Nx_/10; if(lmod<1) lmod=1;
     733 if(lp_>0) cout<<"--- ComputeFourier at z="<<compute_pk_redsh_ref_<<" ---"<<endl;
     734 long lmod = Nx_/20; if(lmod<1) lmod=1;
    698735 for(long i=0;i<Nx_;i++) {
    699    long ii = (i>Nx_/2) ? Nx_-i : i;
    700    double kx = ii*Dkx_;  kx *= kx;
    701    if(lp_>0 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
     736   double kx = Kx(i);  kx *= kx;
     737   if(lp_>0 && i%lmod==0) cout<<"i="<<i<<"\t nx-i="<<Nx_-i<<"\t kx="<<kx<<"\t pk="<<pk_at_z(kx)<<endl;
    702738   for(long j=0;j<Ny_;j++) {
    703      long jj = (j>Ny_/2) ? Ny_-j : j;
    704      double ky = jj*Dky_;  ky *= ky;
     739     double ky = Ky(j);  ky *= ky;
    705740     for(long l=0;l<NCz_;l++) {
    706        double kz = l*Dkz_;  kz *= kz;
     741       double kz = Kz(l);  kz *= kz;
    707742       if(i==0 && j==0 && l==0) continue; // Suppression du continu
    708743       double k = sqrt(kx+ky+kz);
    … …  
    728763// Take into account the real and complexe conjugate coefficients
    729764// because we want a realization of a real data in real space
     765// On ecrit que: P(k_x,k_y,k_z) = P(-k_x,-k_y,-k_z)
     766//               avec k_x = i, -k_x = N_x - i  etc...
    730767{
    731768 if(lp_>1) cout<<"...managing coefficients"<<endl;
    … …  
    800837 if(lp_>1) cout<<"Number of forced conjugate hors cont+nyq ="<<nconj2<<endl;
    801838
    802  if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-8<<endl;
     839 if(lp_>1) cout<<"Check: ddl= "<<NRtot_<<" =?= "<<2*(Nx_*Ny_*NCz_-nconj1-nconj2)-nreal<<endl;
    803840
    804841 return nreal+nconj1+nconj2;
    … …  
    840877
    841878 for(long i=0;i<Nx_;i++) {
    842    long ii = (i>Nx_/2) ? Nx_-i : i;
    843    double kx = ii*Dkx_ *Dx_/2;
     879   double kx = Kx(i) *Dx_/2;
    844880   double pk_x = pixelfilter(kx);
    845881   for(long j=0;j<Ny_;j++) {
    846      long jj = (j>Ny_/2) ? Ny_-j : j;
    847      double ky = jj*Dky_ *Dy_/2;
     882     double ky = Ky(j) *Dy_/2;
    848883     double pk_y = pixelfilter(ky);
    849884     for(long l=0;l<NCz_;l++) {
    850        double kz = l*Dkz_ *Dz_/2;
     885       double kz = Kz(l) *Dz_/2;
    851886       double pk_z =  pixelfilter(kz);
    852887       T_(l,j,i) *= pk_x*pk_y*pk_z;
    … …  
    855890 }
    856891
     892}
     893
     894//-------------------------------------------------------------------
     895void GeneFluct3D::ToVelTrans(void)
     896{
     897  cout<<"-------------------- TO BE VERIFIED ToVelTrans"<<endl;
     898
     899 double zpk = compute_pk_redsh_ref_;
     900 double dpsd = -cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
     901 if(lp_>0) cout<<"--- ToVelTrans --- at z="<<zpk<<", D'/D="<<dpsd
     902               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
     903 check_array_alloc();
     904
     905 for(long i=0;i<Nx_;i++) {
     906   double kx = Kx(i);
     907   for(long j=0;j<Ny_;j++) {
     908     double ky = Ky(j);
     909     double kt2 = kx*kx + ky*ky;
     910     for(long l=0;l<NCz_;l++) {
     911       double kz = Kz(l);
     912       double k2 = kt2 + kz*kz;
     913       if(k2<=0.) continue;
     914       T_(l,j,i) *= complex<double>(0.,dpsd*kz/k2);
     915     }
     916   }
     917 }
    857918}
    858919
    … …  
    879940
    880941 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
    881  //unsigned short ok;
    882 
    883  //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
     942
    884943 for(long i=0;i<Nx_;i++) {
    885944   double dx2 = DXcom(i); dx2 *= dx2;
    … …  
    890949       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
    891950         else dz = fabs(dz); // tous les plans Z au meme redshift
    892        double z = interpinv(dz);
    893        //CHECK: hgr.Add(z);
    894        double dzgr = (*growth_)(z);   // interpolation par morceau
    895        //double dzgr = growth_->Linear(z,ok);  // interpolation lineaire
    896        //double dzgr = growth_->Parab(z,ok);  // interpolation parabolique
     951       double z = interpinv(dz);   // interpolation par morceau
     952       double dzgr = (*growth_)(z);
    897953       int_8 ip = IndexR(i,j,l);
    898954       data_[ip] *= dzgr;
    … …  
    901957 }
    902958
    903  //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
     959}
     960
     961//-------------------------------------------------------------------
     962void GeneFluct3D::ApplyDerGrowthFactor(int type_evol)
     963// Apply Conformal derivative of Growth to real space for transverse velocity cube
     964{
     965  cout<<"-------------------- TO BE VERIFIED ApplyDerGrowthFactor"<<endl;
     966
     967 if(lp_>0) cout<<"--- ApplyDerGrowthFactor: evol="<<type_evol<<endl;
     968 check_array_alloc();
     969
     970 if(growth_ == NULL) {
     971   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first";
     972   cout<<bla<<endl; throw ParmError(bla);
     973 }
     974 if(type_evol<1 || type_evol>2) {
     975   const char *bla = "GeneFluct3D::ApplyGrowthFactor_Error: bad type_evol value";
     976   cout<<bla<<endl; throw ParmError(bla);
     977 }
     978
     979 double zpk = compute_pk_redsh_ref_;
     980 double dpsd_orig = - cosmo_->H(zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
     981
     982 InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
     983 InterpFunc interdpd(zredmin_dpsd_,zredmax_dpsd_,dpsdfrzred_);
     984
     985 for(long i=0;i<Nx_;i++) {
     986   double dx2 = DXcom(i); dx2 *= dx2;
     987   for(long j=0;j<Ny_;j++) {
     988     double dy2 = DYcom(j); dy2 *= dy2;
     989     for(long l=0;l<Nz_;l++) {
     990       double dz = DZcom(l);
     991       if(type_evol==1) dz = sqrt(dx2+dy2+dz*dz);
     992         else dz = fabs(dz); // tous les plans Z au meme redshift
     993       double z = interpinv(dz);   // interpolation par morceau
     994       double dpsd = interdpd(z);
     995       int_8 ip = IndexR(i,j,l);
     996       data_[ip] *= dpsd / dpsd_orig;
     997     }
     998   }
     999 }
    9041000
    9051001}
    … …  
    9481044 // Attention a l'ordre
    9491045 for(long i=0;i<Nx_;i++) {
    950    long ii = (i>Nx_/2) ? Nx_-i : i;
    951    double kx = ii*Dkx_;  kx *= kx;
     1046   double kx = Kx(i);  kx *= kx;
    9521047   for(long j=0;j<Ny_;j++) {
    953      long jj = (j>Ny_/2) ? Ny_-j : j;
    954      double ky = jj*Dky_;  ky *= ky;
     1048     double ky = Ky(j);  ky *= ky;
    9551049     for(long l=0;l<NCz_;l++) {
    956        double kz = l*Dkz_;
     1050       double kz = Kz(l);
    9571051       double k = sqrt(kx+ky+kz*kz);
    9581052       double pk = MODULE2(T_(l,j,i));
    … …  
    9801074 // Attention a l'ordre
    9811075 for(long i=0;i<Nx_;i++) {
    982    long ii = (i>Nx_/2) ? Nx_-i : i;
    983    double kx = ii*Dkx_;  kx *= kx;
     1076   double kx = Kx(i);  kx *= kx;
    9841077   for(long j=0;j<Ny_;j++) {
    985      long jj = (j>Ny_/2) ? Ny_-j : j;
    986      double ky = jj*Dky_;  ky *= ky;
     1078     double ky = Ky(j);  ky *= ky;
    9871079     double kt = sqrt(kx+ky);
    9881080     for(long l=0;l<NCz_;l++) {
    989        double kz = l*Dkz_;
     1081       double kz = Kz(l);
    9901082       double pk = MODULE2(T_(l,j,i));
    9911083       herr.Add(kt,kz,pk);
    … …  
    10261118 // Attention a l'ordre
    10271119 for(long i=0;i<Nx_;i++) {
    1028    long ii = (i>Nx_/2) ? Nx_-i : i;
    1029    double kx = ii*Dkx_;
     1120   double kx = Kx(i);
    10301121   double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
    10311122   kx *= kx; fx *= fx;
    10321123   for(long j=0;j<Ny_;j++) {
    1033      long jj = (j>Ny_/2) ? Ny_-j : j;
    1034      double ky = jj*Dky_;
     1124     double ky = Ky(j);
    10351125     double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
    10361126     ky *= ky; fy *= fy;
    10371127     for(long l=0;l<NCz_;l++) {
    1038        double kz = l*Dkz_;
     1128       double kz = Kz(l);
    10391129       double k = sqrt(kx+ky+kz*kz);
    10401130       double pk = MODULE2(T_(l,j,i)) - sigma2;
    … …  
    10731163 // Attention a l'ordre
    10741164 for(long i=0;i<Nx_;i++) {
    1075    long ii = (i>Nx_/2) ? Nx_-i : i;
    1076    double kx = ii*Dkx_;
     1165   double kx = Kx(i);
    10771166   double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
    10781167   kx *= kx; fx *= fx;
    10791168   for(long j=0;j<Ny_;j++) {
    1080      long jj = (j>Ny_/2) ? Ny_-j : j;
    1081      double ky = jj*Dky_;
     1169     double ky = Ky(j);
    10821170     double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
    10831171     ky *= ky; fy *= fy;
    10841172     double kt = sqrt(kx+ky);
    10851173     for(long l=0;l<NCz_;l++) {
    1086        double kz = l*Dkz_;
     1174       double kz = Kz(l);
    10871175       double pk = MODULE2(T_(l,j,i)) - sigma2;
    10881176       double fz = (pixcor) ? vfz(l): 1.;

  • trunk/Cosmo/SimLSS/genefluct3d.h

    r3662 r3768  
    1818#include "cosmocalc.h"
    1919#include "pkspectrum.h"
     20#include "schechter.h"
    2021
    2122#define WITH_FFTW_THREAD
    2223
    23 //NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE #define GEN3D_FLOAT
    24 
     24//NE PAS DECOMMENTER, UTILISEZ LA MAKEFILE ("g++ -DGEN3D_FLOAT ...")
    2525#if defined(GEN3D_FLOAT)
    2626#define GEN3D_TYPE r_4
    … …  
    7777  vector<long> GetNpix(void) {return N_;}
    7878  int_8 NPix(void) {return NRtot_;}
    79   long GetNx(void) {return Nx_;}
    80   long GetNy(void) {return Ny_;}
    81   long GetNz(void) {return Nz_;}
     79  long Nx(void) {return Nx_;}
     80  long Ny(void) {return Ny_;}
     81  long Nz(void) {return Nz_;}
     82  double Dx(void) {return Dx_;}
     83  double Dy(void) {return Dy_;}
     84  double Dz(void) {return Dz_;}
    8285
    8386  // Return |K_i| module relative to pixel indices
    … …  
    101104  double GetKTincMax(void) {return max(Dk_[0],Dk_[1]);}
    102105
    103   void ComputeFourier0(GenericFunc& pk_at_z);
    104   void ComputeFourier(GenericFunc& pk_at_z);
     106  double Zref(void) {return redsh_ref_;}
     107  double ZrefPk(void) {return compute_pk_redsh_ref_;}
     108  double Dref(void) {return growth_ref_;}
     109  double DrefPk(void)
     110    {return (growth_!=NULL && compute_pk_redsh_ref_>=0.) ? (*growth_)(compute_pk_redsh_ref_): -999.;}
     111
     112  void ComputeFourier0(PkSpectrumZ& pk_at_z);
     113  void ComputeFourier(PkSpectrumZ& pk_at_z);
    105114  void FilterByPixel(void);
     115  void ToVelTrans(void);
    106116
    107117  void ComputeReal(void);
    108118  void ApplyGrowthFactor(int type_evol=1);
     119  void ApplyDerGrowthFactor(int type_evol=1);
    109120
    110121  void ReComputeFourier(void);
    … …  
    188199  CosmoCalc *cosmo_;
    189200  GrowthFactor *growth_;
     201  double good_dzinc_;
    190202  double redsh_ref_,kredsh_ref_,dred_ref_;
    191   double loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
     203  double compute_pk_redsh_ref_;
     204  double h_ref, loscom_ref_,dtrc_ref_, dlum_ref_, dang_ref_;
     205  double growth_ref_, dsdz_growth_ref_;
    192206  double nu_ref_, dnu_ref_ ;
    193207  double xobs_[3];
     208
    194209  double loscom_min_, loscom_max_;
    195210  vector<double> zred_, loscom_;
    … …  
    197212  vector<double> loscom2zred_;
    198213
     214  double zredmin_dpsd_, zredmax_dpsd_;
     215  vector<double> dpsdfrzred_;
     216
    199217};
    200218

  • trunk/Cosmo/SimLSS/pkspectrum.cc

    r3662 r3768  
    399399}
    400400
     401double GrowthFactor::DsDz(double z,double dzinc)
     402// y-y0 = a*(x-x0)^2 + b*(x-x0)
     403// dy = a*dx^2 + b*dx   avec  dx = x-x0
     404// dy'(dx) = 2*a*dx + b -> pour x=x0 on a dy'(0) = b
     405{
     406 if(z<0. || dzinc<=0.) {
     407   cout<<"GrowthFactor::DsDz: z<0 or dzinc<=0. !"<<endl;
     408   throw ParmError("GrowthFactor::DsDz: z<0 or dzinc<=0. !");
     409 }
     410
     411 double z1, z2;
     412 if(z>dzinc/2.) {
     413   // cas ou z est suffisamment loin de zero
     414   // resolution avec 2 points encadrant x0=z
     415   z1 = z - dzinc; if(z1<0.) z1 = 0.;
     416   z2 = z + dzinc;
     417 } else {
     418   // cas ou z est proche de zero
     419   // resolution avec 2 points au dessus, x0=z1
     420   z1 = z + dzinc;
     421   z2 = z + 2.*dzinc;
     422 }
     423
     424 double gz  = (*this)(z);
     425 double dgz1 = (*this)(z1) - gz;
     426 double dgz2 = (*this)(z2) - gz;
     427
     428 z1 -= z;
     429 z2 -= z;
     430 return (dgz2*z1*z1 - dgz1*z2*z2)/(z1*z2*(z1-z2));
     431}
     432
    401433void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
    402434{

  • trunk/Cosmo/SimLSS/pkspectrum.h

    r3381 r3768  
    7777  virtual ~GrowthFactor(void);
    7878  virtual double operator() (double z);
     79  virtual double DsDz(double z,double dzinc=0.01);
    7980  void SetParTo(double OmegaMatter0,double OmegaLambda0);
    8081  bool SetParTo(double OmegaMatter0);
Note: See TracChangeset for help on using the changeset viewer.