source: Sophya/trunk/Cosmo/SimLSS/cmvtvarspec.cc@ 3989

Last change on this file since 3989 was 3805, checked in by cmv, 15 years ago

Refonte totale de la structure des classes InitialSpectrum, TransfertFunction, GrowthFactor, PkSpectrum et classes tabulate pour lecture des fichiers CAMB, cmv 24/07/2010

File size: 5.4 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <iostream>
4#include <stdlib.h>
5#include <stdio.h>
6#include <string.h>
7#include <math.h>
8#include <unistd.h>
9#include "timing.h"
10#include "ntuple.h"
11
12#include "constcosmo.h"
13#include "pkspectrum.h"
14
15void usage(void);
16void usage(void)
17{
18 cout<<"cmvtvarspec -z zref -R R -n as,ns -k kmin,kmax -i eps,npt,fracmax"<<endl;
19}
20
21
22int main(int narg,char *arg[])
23{
24 double ob0 = 0.0444356;
25 // -- WMAP
26 double h100=0.71, om0=0.267804, ol0=0.73;
27 // -- ouvert matter only
28 //double h100=0.71, om0=0.3, ol0=0.;
29 // -- plat matter only
30 //double h100=0.71, om0=1., ol0=0.;
31
32 double ns = 1., as = 1.;
33
34 double zref = 0.5;
35 double R=8./h100;
36
37 double kmin=1e-5,kmax=1000.;
38 int npt = 10000;
39 double eps=1.e-3, fracmax=1e-4;
40
41 char c;
42 while((c = getopt(narg,arg,"hz:R:n:k:i:")) != -1) {
43 switch (c) {
44 case 'z' :
45 sscanf(optarg,"%lf",&zref);
46 break;
47 case 'R' :
48 sscanf(optarg,"%lf",&R);
49 if(R<0.) R = -R/h100;
50 break;
51 case 'n' :
52 sscanf(optarg,"%lf,%lf",&as,&ns);
53 break;
54 case 'k' :
55 sscanf(optarg,"%lf,%lf",&kmin,&kmax);
56 break;
57 case 'i' :
58 sscanf(optarg,"%lf,%d,%lf",&eps,&npt,&fracmax);
59 break;
60 case 'h' :
61 default :
62 usage(); return -1;
63 }
64 }
65
66 double Rg=R/sqrt(5.);
67 double lkmin=log10(kmin), lkmax=log10(kmax), dlk=(lkmax-lkmin)/npt;
68
69 cout<<"zref="<<zref<<endl;
70 cout<<"as="<<as<<" ns="<<ns<<endl;
71 cout<<"R="<<R<<", Rg="<<Rg<<endl;
72 cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<")"
73 <<", dlk="<<dlk<<" npt="<<npt<<endl;
74 cout<<"eps="<<eps<<" fracmax="<<fracmax<<endl;
75
76 //-----------------------------------------------------------------
77 InitialPowerLaw pkini(ns,as);
78
79 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
80 //tf.SetNoOscEnv(2);
81
82 GrowthEisenstein d1(om0,ol0);
83
84 PkSpecCalc pkz(pkini,tf,d1,zref);
85
86 //-----------------------------------------------------------------
87 double k1,k2;
88 int rc = 0;
89 //----
90 cout<<endl<<"Filtrage top_hat"<<endl;
91 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
92 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
93 double pkmax_th = varpk_th(kfind_th);
94 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
95 k1=kmin, k2=kmax;
96 rc = varpk_th.FindLimits(pkmax_th*fracmax,k1,k2,eps);
97 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
98
99 varpk_th.SetInteg(0.01,dlk,-1.,4);
100 cout<<"varpk_th="<<varpk_th.Variance(k1,k2)<<endl;
101
102 //----
103 cout<<endl<<"Filtrage gaussien"<<endl;
104 VarianceSpectrum varpk_ga(pkz,Rg,VarianceSpectrum::GAUSSIAN);
105 double kfind_ga = varpk_ga.FindMaximum(kmin,kmax,eps);
106 double pkmax_ga = varpk_ga(kfind_ga);
107 cout<<"kfind_ga = "<<kfind_ga<<" -> "<<log10(kfind_ga)<<", integrand="<<pkmax_ga<<endl;
108 k1=kmin, k2=kmax;
109 rc = varpk_ga.FindLimits(pkmax_ga*fracmax,k1,k2,eps);
110 cout<<"limit_ga: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
111
112 varpk_ga.SetInteg(0.01,dlk,-1.,4);
113 cout<<"varpk_ga="<<varpk_ga.Variance(k1,k2)<<endl;
114
115 //----
116 cout<<endl<<"Filtrage 1 (integrale du spectre)"<<endl;
117 VarianceSpectrum varpk_int(pkz,Rg,VarianceSpectrum::NOFILTER);
118 double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
119 double pkmax_int = varpk_int(kfind_int);
120 cout<<"kfind_int = "<<kfind_int<<" -> "<<log10(kfind_int)<<", integrand="<<pkmax_int<<endl;
121 k1=kmin, k2=kmax;
122 rc = varpk_int.FindLimits(pkmax_int*fracmax,k1,k2,eps);
123 cout<<"limit_int: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;
124
125 varpk_int.SetInteg(0.01,dlk,-1.,4);
126 cout<<"varpk_int="<<varpk_int.Variance(k1,k2)<<endl;
127
128//-----------------------------------------------------------------
129 const int n = 2;
130 const char *vname[n] = {"k","pk"};
131 NTuple nt(n,vname);
132 double xnt[n];
133
134 for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
135 double k = pow(10.,lk);
136 xnt[0] = k;
137 xnt[1] = pkz(k);
138 nt.Fill(xnt);
139 }
140
141 cout<<">>>> Ecriture"<<endl;
142 string tag = "cmvtvarspec.ppf"; POutPersist pos(tag); tag = "nt"; pos.PutObject(nt,tag);
143
144 return 0;
145}
146
147/*
148openppf cmvtvarspec.ppf
149
150#### Spectre a z
151n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"
152
153# L'approximation pour j1(x) (top-hat) et la gaussienne equivalente
154set l1 1e-6
155set l2 10
156set j1sx pow(3*(sin(x)-x*cos(x))/(x*x*x),2.)
157set j1sxa3 (1.-x*x/5.)
158set j1sxa5 (1.-x*x/5.*(1.-3.*x*x/35.))
159set j1sxa7 (1.-x*x/5.*(1.-3.*x*x/35.*(1.-4.*x*x/81.)))
160func $j1sx $l1 $l2 1000
161func $j1sxa3 $l1 $l2 1000 "red same"
162func $j1sxa5 $l1 $l2 1000 "blue same"
163func $j1sxa7 $l1 $l2 1000 "orange same"
164set gauss pow(exp(-(x/sqrt(5.))*(x/sqrt(5.))/2.),2.)
165func $gauss $l1 $l2 1000 "green same"
166
167func $j1sx-$j1sxa3 $l1 $l2 1000 "red"
168func $j1sx-$j1sxa5 $l1 $l2 1000 "blue same"
169func $j1sx-$j1sxa7 $l1 $l2 1000 "orange same"
170func $j1sx-$gauss $l1 $l2 1000 "green same"
171addline $l1 0 $l2 0
172
173# la fonction "Pk" a integrer
174zone
175set D (k*k*pk/(2.*M_PI*M_PI))
176n/plot nt.$D%log10(k) ! ! "nsta crossmarker3"
177
178# La fonction a integrer en dk lineaire (top-hat et gaussien)
179set R 10
180set RG ($R/sqrt(5.))
181set x (k*$R)
182set j1sx pow(3*(sin($x)-$x*cos($x))/($x*$x*$x),2.)
183set y (k*$RG)
184set gauss pow(exp(-$y*$y/2.),2.)
185
186n/plot nt.$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
187n/plot nt.$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"
188
189# La fonction a integrer en d(log(x)) (top-hat et gaussien)
190# c'est ce que fait Variance()
191n/plot nt.k*$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
192n/plot nt.k*$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"
193n/plot nt.k*$D%log10(k) ! ! "nsta crossmarker3 same blue"
194
195*/
Note: See TracBrowser for help on using the repository browser.