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