source: Sophya/trunk/Cosmo/SimLSS/cmvtstpk.cc@ 4061

Last change on this file since 4061 was 4061, checked in by cmv, 13 years ago

cmv, 17/04/2012
File size: 8.1 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 "histos.h"
11#include "tvector.h"
12#include "ntuple.h"
13#include "perandom.h"
14
15#include "constcosmo.h"
16#include "pkspectrum.h"
17#include "geneutils.h"
18
19void usage(void);
20void usage(void) {
21cout
22<<"cmvtstpk [options] z_redshift"<<endl
23<<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
24<<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
25<<" -s scale : on normalise le spectre avec scale"<<endl
26<<" -w : write spectra on ASCII file"<<endl
27<<" -F filename : read also CAMB file"<<endl;
28}
29
30int main(int narg,char *arg[])
31{
32 // -- WMAP
33 double h100=0.71, Om0=0.267804, Ob0 = 0.0444356, Ol0=0.73,w0=-1.;
34
35 double ns = 1., as = 1.;
36 int npt = 10000;
37 double lkmin = -3., lkmax=2.;
38 double scale = 2.54499e+07; // normalisation du spectre a z=0 selon SDSS
39 bool wrascii = false;
40 string fcmbfile = "";
41
42 char c;
43 while((c = getopt(narg,arg,"hwk:s:H:M:B:L:F:")) != -1) {
44 switch (c) {
45 case 'H' :
46 sscanf(optarg,"%lf",&h100);
47 break;
48 case 'M' :
49 sscanf(optarg,"%lf",&Om0);
50 break;
51 case 'B' :
52 sscanf(optarg,"%lf",&Ob0);
53 break;
54 case 'L' :
55 sscanf(optarg,"%lf,%lf",&Ol0,&w0);
56 break;
57 case 'k' :
58 sscanf(optarg,"%d,%lf,%lf",&npt,&lkmin,&lkmax);
59 if(npt<=0) npt=1000;
60 if(lkmax<lkmin) {lkmin=-4.; lkmax=2.;}
61 break;
62 case 's' :
63 sscanf(optarg,"%lf",&scale);
64 break;
65 case 'w' :
66 wrascii = true;
67 break;
68 case 'F' :
69 fcmbfile = optarg;
70 break;
71 case 'h' :
72 default :
73 usage(); return -1;
74 }
75 }
76 double dlk=(lkmax-lkmin)/npt;
77
78 // Fill z value into list
79 double zval = 0.;
80 if(optind<narg) zval = atof(arg[optind]);
81
82 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Ocdm="<<Om0-Ob0<<" Ol0="<<Ol0<<" w0="<<w0<<endl;
83 cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
84 cout<<"scale="<<scale<<endl;
85 cout<<"zval="<<zval<<endl;
86 cout<<"fcmbfile="<<fcmbfile<<endl;
87
88 //--------------------------
89 InitialPowerLaw pkini(ns,as);
90 pkini.SetNorm(scale);
91
92 TransfertEisenstein tf(h100,Om0-Ob0,Ob0,T_CMB_Par,false);
93 double tf0 = tf(0);
94 cout<<"kpeak="<<tf.KPeak()<<endl;
95 TransfertEisenstein tfnosc2(tf); tfnosc2.SetNoOscEnv(2);
96 double tfnosc20 = tfnosc2(0);
97 TransfertEisenstein tfnosc1(tf); tfnosc1.SetNoOscEnv(1);
98 double tfnosc10 = tfnosc1(0);
99 TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
100 double tfnob0 = tfnob(0);
101 cout<<"tf(0)="<<tf0<<" tfnosc2(0)="<<tfnosc20<<" tfnosc1(0)="<<tfnosc10<<" tfnob(0)="<<tfnob0<<endl;
102
103 GrowthEisenstein d1(Om0,Ol0);
104 cout<<"GrowthFactor: D1(0)="<<d1(0.)<<" D1("<<zval<<") = "<<d1(zval)<<endl;
105
106 PkSpecCalc pk0(pkini,tf,d1,0.);
107 PkSpecCalc pk0nosc2(pkini,tfnosc2,d1,0.);
108 PkSpecCalc pk0nosc1(pkini,tfnosc1,d1,0.);
109 PkSpecCalc pk0nob(pkini,tfnob,d1,0.);
110
111 PkSpecCalc pkz(pkini,tf,d1,zval);
112 PkSpecCalc pkznosc2(pkini,tfnosc2,d1,zval);
113 PkSpecCalc pkznosc1(pkini,tfnosc1,d1,zval);
114 PkSpecCalc pkznob(pkini,tfnob,d1,zval);
115
116 PkTabulate* pkcamb = NULL;
117 if(fcmbfile.size()>0) {
118 cout<<endl<<"PkTabulate for CAMB"<<endl;
119 pkcamb = new PkTabulate;
120 pkcamb->ReadCAMB(fcmbfile,h100,0.);
121 pkcamb->SetInterpTyp(2);
122 if(pkcamb->NPoints()==0) {delete pkcamb; pkcamb = NULL;}
123 }
124
125 //--------------------------
126 Histo hd1(0.,20.,10000); hd1.ReCenterBin();
127 FuncToHisto(d1,hd1,false);
128
129 Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
130 FuncToHisto(pkz,hpkz,true);
131 TVector<r_8> vpkz(npt);
132 FuncToVec(pkz,vpkz,lkmin,lkmax,true);
133
134 FunRan talea(hpkz,true);
135 Histo halea(hpkz); halea.Zero();
136 int nalea = 100000;
137 for(int i=0;i<nalea;i++) halea.Add(talea.Random());
138 halea *= hpkz.Sum()/halea.Sum();
139
140 //--------------------------
141 const int n = 15;
142 const char *vname[n] = {
143 "k","pkini",
144 "tf","pk0","pk",
145 "tfnosc2","pk0nosc2","pknosc2",
146 "tfnosc1","pk0nosc1","pknosc1",
147 "tfnob","pk0nob","pknob","pktab"
148 };
149 NTuple nt(n,vname);
150 double xnt[n];
151
152 for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
153 double k = pow(10.,lk);
154 xnt[0] = k;
155 xnt[1] = pkini(k);
156 xnt[2] = tf(k);
157 xnt[3] = pk0(k);
158 xnt[4] = pkz(k,zval);
159 xnt[5] = tfnosc2(k);
160 xnt[6] = pk0nosc2(k);
161 xnt[7] = pkznosc2(k,zval);
162 xnt[8] = tfnosc1(k);
163 xnt[9] = pk0nosc1(k);
164 xnt[10] = pkznosc1(k,zval);
165 xnt[11] = tfnob(k);
166 xnt[12] = pk0nob(k);
167 xnt[13] = pkznob(k,zval);
168 if(pkcamb != NULL) xnt[14] = (*pkcamb)(k); else xnt[14] = 0;
169 nt.Fill(xnt);
170 }
171
172 if(pkcamb != NULL) {delete pkcamb; pkcamb = NULL;}
173
174 //--------------------------
175 if(wrascii) {
176 cout<<">>>> ASCII"<<endl;
177 FILE * fdata = fopen("cmvtstpk.data","w");
178 fprintf(fdata,"# z= %g : k(Mpc^-1) pkini, tf pk0 pkz, tfnosc2 pk0nosc2 pkznosc2, ",zval);
179 fprintf(fdata,"tfnosc1 pk0nosc1 pkznosc1, tfnob pk0nob pkznob\n");
180 int n = 0;
181 for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
182 double k = pow(10.,lk);
183 fprintf(fdata,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
184 n,k,pkini(k),
185 tf(k),pk0(k),pkz(k,zval),
186 tfnosc2(k),pk0nosc2(k),pkznosc2(k,zval),
187 tfnosc1(k),pk0nosc1(k),pkznosc1(k,zval),
188 tfnob(k),pk0nob(k),pkznob(k,zval));
189 n++;
190 }
191 fclose(fdata);
192 }
193
194 //--------------------------
195 cout<<">>>> Ecriture"<<endl;
196 string tag = "cmvtstpk.ppf";
197 POutPersist pos(tag);
198 tag = "nt"; pos.PutObject(nt,tag);
199 tag = "hd1"; pos.PutObject(hd1,tag);
200 tag = "hpkz"; pos.PutObject(hpkz,tag);
201 tag = "vpkz"; pos.PutObject(vpkz,tag);
202 tag = "halea"; pos.PutObject(halea,tag);
203 return 0;
204}
205
206/*
207openppf cmvtstpk.ppf
208
209#### growth-factor
210zone
211n/plot hd1.val%x ! ! "nsta connectpoints"
212n/plot hd1.1./(1.+x)%x ! ! "nsta connectpoints same red"
213
214#### Spectre initial
215zone
216n/plot nt.log10(pkini)%log10(k) pkini>0.&&k>0. ! "nsta connectpoints"
217
218#### fct de transfert
219zone
220n/plot nt.tf%k k>0&&tf>0 ! "nsta connectpoints logx logy"
221n/plot nt.tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints same red logx logy"
222n/plot nt.tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx logy"
223n/plot nt.tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx logy"
224
225n/plot nt.tf/tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints red logx"
226n/plot nt.tf/tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx"
227n/plot nt.tf/tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx"
228addline -10 1 10 1
229
230#### Spectre a z=0
231zone
232n/plot nt.pk0%k k>0&&pk0>0 ! "nsta connectpoints logx logy"
233n/plot nt.pk0nosc2%k k>0&&pk0nosc2>0 ! "nsta connectpoints same red logx logy"
234n/plot nt.pk0nosc1%k k>0&&pk0nosc1>0 ! "nsta connectpoints same blue logx logy"
235n/plot nt.pk0nob%k k>0&&pk0nob>0 ! "nsta connectpoints same green logx logy"
236
237# Check
238zone 2 2
239n/plot nt.pk0/pkini-tf*tf%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
240n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
241n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
242n/plot nt.pk0nob/pkini-tfnob*tfnob%k k>0&&pkini>0 ! "nsta crossmarker3 logx"
243
244#### Spectre a z
245zone
246n/plot nt.pk%k k>0&&pk>0 ! "nsta connectpoints logx logy"
247n/plot nt.pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx logy"
248n/plot nt.pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx logy"
249n/plot nt.pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx logy"
250
251n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints red logx"
252n/plot nt.pk/pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx"
253n/plot nt.pk/pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx"
254addline -10 1 10 1
255
256#### Le spectre version Delta^2
257set D2 k*k*k*pk/(2*M_PI*M_PI)
258n/plot nt.$D2%k k>0 ! "nsta crossmarker3 connectpoints logx"
259
260#### Test des transferts dans Histo et TVector
261zone 1 2
262n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"
263disp hpkz "same red"
264
265disp vpkz
266c++exec cout<<hpkz(0)<<" =?= "<<vpkz(0)<<endl;
267
268zone
269disp hpkz
270disp halea "same red"
271
272#### Le spectre tabule
273n/plot nt.pk%k pk>0&&k>0 ! "nsta connectpoints logx logy"
274n/plot nt.pktab%k pktab>0&&k>0 ! "nsta connectpoints same red logx logy"
275
276n/plot nt.(pktab-pk)%k k>0&&pk>0 ! "nsta connectpoints"
277n/plot nt.(pktab-pk)/pk%k k>0&&pk>0 ! "nsta connectpoints"
278
279n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints logx"
280n/plot nt.pktab/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx"
281
282*/
Note: See TracBrowser for help on using the repository browser.