Changeset 3120 in Sophya for trunk/Cosmo


Ignore:
Timestamp:
Dec 21, 2006, 4:45:09 PM (19 years ago)
Author:
cmv
Message:

Suite de la mise dans la base cvs cmv 21/12/2006

Location:
trunk/Cosmo/SimLSS
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/Makefile

    r3115 r3120  
    1919       $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
    2020       $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
    21        $(EXE)cmvtluc $(EXE)cmvtpoisson $(EXE)cmvconchprof
     21       $(EXE)cmvtpoisson $(EXE)cmvconchprof
     22#$(EXE)cmvtluc
    2223 
    2324PROGSOBJ = \
    … …  
    2526          $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
    2627          $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o \
    27           $(OBJ)cmvtluc.o $(OBJ)cmvtpoisson.o $(OBJ)cmvconchprof.o
     28          $(OBJ)cmvtpoisson.o $(OBJ)cmvconchprof.o $(OBJ)cmvtluc.o
    2829 
    2930LIBROBJ = \

  • trunk/Cosmo/SimLSS/cmvconchprof.cc

    r3115 r3120  
    8888
    8989disp hpconc
     90
    9091n/plot hpconc.val%log10(x) x>0 ! "connectpoints"
     92n/plot hpconc.err%log10(x) x>0 ! "connectpoints same red"
     93
     94n/plot hpconc.err/val%log10(x) x>0&&val>0 ! "connectpoints"
    9195*/

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3115 r3120  
    1212#include "srandgen.h"
    1313#include "perandom.h"
     14#include "fabtwriter.h"
    1415
    1516#include "arrctcast.h"
    … …  
    3132     <<" -z nz,dz : taille en z (npix,Mpc)"<<endl
    3233     <<" -Z zref : redshift median"<<endl
    33      <<" -s scale : sigma du bruit en Msol"<<endl
     34     <<" -s snoise : sigma du bruit en Msol"<<endl
    3435     <<endl;
    3536}
    … …  
    163164 pkz.SetZ(zref);
    164165
    165  Histo hpkz (lkmin,lkmax,npt); hpkz.ReCenterBin();
     166 Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
    166167 FuncToHisto(pkz,hpkz,true);
    167168 {
    … …  
    330331     <<rs2<<" -> "<<sqrt(rs2)<<endl;
    331332
     333 if(1) {
     334   cout<<"\n---Writing Cube to Fits file"<<endl;
     335   FitsImg3DWriter fwrt("!cmvobserv3dr.fits",FLOAT_IMG,5);
     336   fwrt.WriteKey("ZREF",zref," redshift");
     337   vector<long> n = fluct3d.GetNpix();
     338   fwrt.WriteKey("NX",n[0]," axe transverse 1");
     339   fwrt.WriteKey("NY",n[1]," axe transverse 2");
     340   fwrt.WriteKey("NZ",n[2]," axe longitudinal (redshift)");
     341   vector<r_8> d;
     342   d = fluct3d.GetDinc();
     343   fwrt.WriteKey("DX",d[0]," Mpc");
     344   fwrt.WriteKey("DY",d[1]," Mpc");
     345   fwrt.WriteKey("DZ",d[2]," Mpc");
     346   d = fluct3d.GetKinc();
     347   fwrt.WriteKey("DKX",d[0]," Mpc^-1");
     348   fwrt.WriteKey("DKY",d[1]," Mpc^-1");
     349   fwrt.WriteKey("DKZ",d[2]," Mpc^-1");
     350   fwrt.WriteKey("SNOISE",snoise," Msol");
     351   fwrt.Write(rgen);
     352 }
     353
    332354 cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
    333355 fluct3d.AddNoise2Real(snoise);
    … …  
    364386objaoper pkgen sliceyz 0
    365387
    366 n/plot hpkz.val%x x>0 ! "connectpoints"
    367 n/plot hpkgen.val%log10(x) x>0 ! "same red connectpoints"
    368 n/plot hpkgenf.val%log10(x) x>0 ! "same pink connectpoints"
    369 n/plot hpkrec.val%log10(x) x>0 ! "same blue connectpoints"
     388n/plot hpkz.val%x ! ! "nsta connectpoints"
     389n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
     390n/plot hpkgenf.val%log10(x) x>0 ! "nsta same orange connectpoints"
     391n/plot hpkrec.val%log10(x) x>0 ! "nsta same blue connectpoints"
    370392
    371393set k pow(10.,x)
    … …  
    379401rgensl 0
    380402
     403disp hmdndm
     404disp tirhmdndm
     405addline 0 1 20 1 "red"
     406
    381407 */

  • trunk/Cosmo/SimLSS/cmvtstpk.cc

    r3115 r3120  
    1919void usage(void);
    2020void usage(void) {
    21  cout<<"cmvtstpk -k npt,lkmin,lkmax -s scale zval"<<endl
     21 cout<<"cmvtstpk [options] z_redshift"<<endl
     22     <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
    2223     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
    2324     <<" -s scale : on multiplie pkz par scale"<<endl;
    … …  
    2627int main(int narg,char *arg[])
    2728{
    28  double ob0 = 0.0444356;
     29 double Ob0 = 0.0444356;
    2930 // -- WMAP
    30  double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
     31 double h100=0.71, Om0=0.267804, Ol0=0.73,w0=-1.;
    3132 // -- ouvert matter only
    32  //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
     33 //double h100=0.71, Om0=0.3, Ol0=0.,w0=-1.;
    3334 // -- plat matter only
    34  //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;
     35 //double h100=0.71, Om0=1., Ol0=0.,w0=-1.;
     36
    3537
    3638 double ns = 1., as = 1.;
    3739
    38  int npt = 1000;
    39  double lkmin = -4., lkmax=2.;
     40 int npt = 10000;
     41 double lkmin = -3., lkmax=2.;
    4042 double scale = 1.;
    4143
    4244 char c;
    43   while((c = getopt(narg,arg,"hk:s:")) != -1) {
     45  while((c = getopt(narg,arg,"hk:s:H:M:B:L:")) != -1) {
    4446  switch (c) {
     47  case 'H' :
     48    sscanf(optarg,"%lf",&h100);
     49    break;
     50  case 'M' :
     51    sscanf(optarg,"%lf",&Om0);
     52    break;
     53  case 'B' :
     54    sscanf(optarg,"%lf",&Ob0);
     55    break;
     56  case 'L' :
     57    sscanf(optarg,"%lf,%lf",&Ol0,&w0);
     58    break;
    4559  case 'k' :
    4660    sscanf(optarg,"%d,%lf,%lf",&npt,&lkmin,&lkmax);
    … …  
    6276 if(optind<narg) zval = atof(arg[optind]);
    6377
     78 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Ol0="<<Ol0<<" w0="<<w0<<endl;
    6479 cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
    6580 cout<<"scale="<<scale<<endl;
    … …  
    6984 InitialSpectrum pkini(ns,as);
    7085
    71  TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
     86 TransfertEisenstein tf(h100,Om0-Ob0,Ob0,T_CMB_Par,false);
    7287 cout<<"kpeak="<<tf.KPeak()<<endl;
    7388 TransfertEisenstein tfnosc2(tf); tfnosc2.SetNoOscEnv(2);
    7489 TransfertEisenstein tfnosc1(tf); tfnosc1.SetNoOscEnv(1);
    75  TransfertEisenstein tfnob(h100,om0,0.,T_CMB_Par,true);
    76 
    77  GrowthFactor d1(om0,ol0);
     90 TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
     91
     92 GrowthFactor d1(Om0,Ol0);
    7893 cout<<"GrowthFactor: "<<d1(zval)<<endl;
    7994
    … …  
    89104
    90105 //--------------------------
    91  Histo hd1(0.,1000.,10000); hd1.ReCenterBin();
     106 Histo hd1(0.,20.,10000); hd1.ReCenterBin();
    92107 FuncToHisto(d1,hd1,false);
    93108
    … …  
    99114 FunRan talea(hpkz,true);
    100115 Histo halea(hpkz); halea.Zero();
    101  int nalea = 1000000;
     116 int nalea = 100000;
    102117 for(int i=0;i<nalea;i++) halea.Add(talea.Random());
    103118 halea *= hpkz.Sum()/halea.Sum();
    … …  
    169184#### growth-factor
    170185zone
    171 n/plot hd1.log10(val)%log10(x) val>0.&&x>0. ! "nsta connectpoints"
     186n/plot hd1.val%x ! ! "nsta connectpoints"
    172187
    173188#### Spectre initial
    174189zone
    175 n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta crossmarker3"
     190n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta connectpoints"
     191
     192#### Gestion des abscisses
     193set xk log10(k)
     194set xk k
    176195
    177196#### fct de transfert
    178197zone
    179 n/plot nt.tf%log10(k) ! ! "nsta connectpoints"
    180 n/plot nt.tfnosc2%log10(k) ! ! "nsta connectpoints same red"
    181 n/plot nt.tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
    182 n/plot nt.tfnob%log10(k) ! ! "nsta connectpoints same green"
    183 
    184 n/plot nt.tf/tfnosc2%log10(k) ! ! "nsta connectpoints red"
    185 n/plot nt.tf/tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
    186 n/plot nt.tf/tfnob%log10(k) ! ! "nsta connectpoints same red"
     198n/plot nt.tf%$xk ! ! "nsta connectpoints"
     199n/plot nt.tfnosc2%$xk ! ! "nsta connectpoints same red"
     200n/plot nt.tfnosc1%$xk ! ! "nsta connectpoints same blue"
     201n/plot nt.tfnob%$xk ! ! "nsta connectpoints same green"
     202
     203n/plot nt.tf/tfnosc2%$xk ! ! "nsta connectpoints red"
     204n/plot nt.tf/tfnosc1%$xk ! ! "nsta connectpoints same blue"
     205n/plot nt.tf/tfnob%$xk ! ! "nsta connectpoints same green"
    187206addline -10 1 10 1
    188207
    189208#### Spectre a z=0
    190209zone
    191 n/plot nt.pk0%log10(k) ! ! "nsta connectpoints"
    192 n/plot nt.pk0nosc2%log10(k) ! ! "nsta connectpoints same red"
    193 n/plot nt.pk0nosc1%log10(k) ! ! "nsta connectpoints same blue"
    194 n/plot nt.pk0nob%log10(k) ! ! "nsta connectpoints same green"
     210n/plot nt.pk0%$xk ! ! "nsta connectpoints"
     211n/plot nt.pk0nosc2%$xk ! ! "nsta connectpoints same red"
     212n/plot nt.pk0nosc1%$xk ! ! "nsta connectpoints same blue"
     213n/plot nt.pk0nob%$xk ! ! "nsta connectpoints same green"
    195214
    196215# Check
    197216zone 2 2
    198 n/plot nt.pk0/pkini-tf*tf%log10(k) pkini>0 ! "nsta crossmarker3"
    199 n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%log10(k) pkini>0 ! "nsta crossmarker3"
    200 n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%log10(k) pkini>0 ! "nsta crossmarker3"
    201 n/plot nt.pk0nob/pkini-tfnob*tfnob%log10(k) pkini>0 ! "nsta crossmarker3"
     217n/plot nt.pk0/pkini-tf*tf%$xk pkini>0 ! "nsta crossmarker3"
     218n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%$xk pkini>0 ! "nsta crossmarker3"
     219n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%$xk pkini>0 ! "nsta crossmarker3"
     220n/plot nt.pk0nob/pkini-tfnob*tfnob%$xk pkini>0 ! "nsta crossmarker3"
    202221
    203222#### Spectre a z
    204223zone
    205 n/plot nt.pk%log10(k) ! ! "nsta connectpoints"
    206 n/plot nt.pknosc2%log10(k) ! ! "nsta connectpoints same red"
    207 n/plot nt.pknosc1%log10(k) ! ! "nsta connectpoints same blue"
    208 n/plot nt.pknob%log10(k) ! ! "nsta connectpoints same green"
    209 
    210 n/plot nt.pk/pknosc2%log10(k) ! ! "nsta connectpoints red"
    211 n/plot nt.pk/pknosc1%log10(k) ! ! "nsta connectpoints same blue"
    212 n/plot nt.pk/pknob%log10(k) ! ! "nsta connectpoints same green"
     224n/plot nt.pk%$xk ! ! "nsta connectpoints"
     225n/plot nt.pknosc2%$xk ! ! "nsta connectpoints same red"
     226n/plot nt.pknosc1%$xk ! ! "nsta connectpoints same blue"
     227n/plot nt.pknob%$xk ! ! "nsta connectpoints same green"
     228
     229n/plot nt.pk/pknosc2%$xk ! ! "nsta connectpoints red"
     230n/plot nt.pk/pknosc1%$xk ! ! "nsta connectpoints same blue"
     231n/plot nt.pk/pknob%$xk ! ! "nsta connectpoints same green"
    213232addline -10 1 10 1
    214233
    215234#### Le spectre version Delta^2
    216235set D2 k*k*k*pk/(2*M_PI*M_PI)
    217 n/plot nt.$D2%log10(k)  !  ! "nsta crossmarker3 connectpoints"
     236n/plot nt.$D2%$xk  !  ! "nsta crossmarker3 connectpoints"
    218237
    219238#### Test des transferts dans Histo et TVector

  • trunk/Cosmo/SimLSS/cmvtstsch.cc

    r3115 r3120  
    9494 for(int i=0;i<=5;i++) {
    9595   double num = IntegrateFuncLog(sch,lnx1-i,lnx2+i,perc,dlxinc,dlxmax,glorder);
    96    cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated number = "<<num<<endl;
     96   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
    9797 }
    9898 Histo hdndm(lnx1,lnx2,npt); hdndm.ReCenterBin();
    … …  
    105105 for(int i=0;i<=5;i++) {
    106106   double sum = IntegrateFuncLog(sch,lnx1-i,lnx2+i,perc,dlxinc,dlxmax,glorder);
    107    cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated mass = "<<sum<<" Msol"<<endl;
     107   cout<<"["<<lnx1-i<<","<<lnx2+i<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
    108108 }
    109109 Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();

  • trunk/Cosmo/SimLSS/constcosmo.h

    r3115 r3120  
    1010static const double k_Boltzman_Cst = 1.3806503e-23;  // Boltzman constant SI "J / K"
    1111static const double ProtonMass_Cst = 1.67262171e-24; // Proton mass in "g"
    12 static const double SolarMass_Cst = 1.98844e+33; // Proton mass in "g"
     12static const double SolarMass_Cst = 1.98844e+33; // Solar mass in "g"
    1313//            Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
    1414static const double StefBoltz_Cst = 5.670400e-8; // constante de Stefan-Boltzman in "W/m^2/K^4"

  • trunk/Cosmo/SimLSS/genefluct3d.cc

    r3115 r3120  
    344344void GeneFluct3D::FilterByPixel(void)
    345345// Filtrage par la fonction fenetre du pixel (parallelepipede)
    346 // TF = 1/(dx*dy*dz)*Int[{0,dx},{0,dy},{0,dz}]
     346// TF = 1/(dx*dy*dz)*Int[{-dx/2,dx/2},{-dy/2,dy/2},{-dz/2,dz/2}]
    347347//                   e^(ik_x*x) e^(ik_y*y) e^(ik_z*z) dxdydz
    348 // |TF|^2 = 2*(1-cos(k_x*dx)/dx^2/k_x^2 * (idem dy) * (idem dz)
    349 // et on traite la divergence en K_x = 0:
    350 // 2*(1-cos(k*d)/d^2/k^2 ~= 1 - (k*d)^2/12 *[1 - (k*d)^2/30]
     348//    = 2/(k_x*dx) * sin(k_x*dx/2)  * (idem y) * (idem z)
     349// Gestion divergence en 0: sin(y)/y = 1 - y^2/6*(1-y^2/20)
     350//                          avec y = k_x*dx/2
    351351{
    352352 int lp=2;
    … …  
    355355 for(int i=0;i<Nx_;i++) {
    356356   int ii = (i>Nx_/2) ? Nx_-i : i;
    357    double kx = ii*Dkx_ *Dx_;
    358    double pkx2 = pixelfilter(kx);
     357   double kx = ii*Dkx_ *Dx_/2;
     358   double pkx = pixelfilter(kx);
    359359   for(int j=0;j<Ny_;j++) {
    360360     int jj = (j>Ny_/2) ? Ny_-j : j;
    361      double ky = jj*Dky_ *Dy_;
    362      double pky2 =  pixelfilter(ky);
     361     double ky = jj*Dky_ *Dy_/2;
     362     double pky =  pixelfilter(ky);
    363363     for(int l=0;l<NCz_;l++) {
    364        double kz = l*Dkz_ *Dz_;
    365        double pkz2 =  pixelfilter(kz);
    366        T_(l,j,i) *= pkx2*pky2*pkz2;
     364       double kz = l*Dkz_ *Dz_/2;
     365       double pkz =  pixelfilter(kz);
     366       T_(l,j,i) *= pkx*pky*pkz;
    367367     }
    368368   }

  • trunk/Cosmo/SimLSS/genefluct3d.h

    r3115 r3120  
    3535  vector<r_8> GetKnyq(void)
    3636    {vector<r_8> kn; kn.push_back(Knyqx_); kn.push_back(Knyqy_); kn.push_back(Knyqz_); return kn;}
     37  vector<r_8> GetDinc(void)
     38    {vector<r_8> d; d.push_back(Dx_); d.push_back(Dy_); d.push_back(Dz_); return d;}
     39  vector<long> GetNpix(void)
     40    {vector<long> d; d.push_back(Nx_); d.push_back(Ny_); d.push_back(Nz_); return d;}
    3741
    3842  void ComputeFourier0(void);
    … …  
    6064  int manage_coefficients(void);
    6165  double compute_power_carte(void);
    62   inline double pixelfilter(double kx)
    63     {return (kx<0.05) ? 1.-kx*kx/12.*(1.-kx*kx/30.): 2.*(1.-cos(kx))/(kx*kx);}
     66  inline double pixelfilter(double x)
     67    {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
    6468
    6569

  • trunk/Cosmo/SimLSS/geneutils.h

    r3115 r3120  
    1616  InterpFunc(double xmin,double xmax,vector<double>& y);
    1717  virtual ~InterpFunc(void) { }
     18
     19  double XMin(void) {return _xmin;}
     20  double XMax(void) {return _xmax;}
    1821
    1922  //! Retourne l'element le plus proche de f donnant y=f(x)
Note: See TracChangeset for help on using the changeset viewer.