Commit 7e1bcb90 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

changes for pypico-3.3

parent f933b6fe
......@@ -43,6 +43,9 @@ MnPicoEngine::MnPicoEngine(const Variables& _upar,const int& lmax, const string&
_lmax=lmax;
planck_assert(file_present(picodata),string("missing pico data file: ")+picodata);
pico_init(true);
cout<<"loading pico datafile: " << picodata << endl;
pico=pico_load(const_cast<char*>(picodata.c_str()));
pico_set_verbose(pico,verbose);
......@@ -77,6 +80,9 @@ MnPicoEngine::MnPicoEngine(const Variables& _upar,const int& lmax, const string&
//add default values
PyDict_SetItemString(dict,"force",Py_True);
PyDict_SetItemString(dict,"scalar_nrun(1)",PyFloat_FromDouble(0));
//kpivot=0.05
PyDict_SetItemString(dict,"pivot_scalar",PyFloat_FromDouble(0.05));
//output keys
outputs=PyList_New(6);
......@@ -105,18 +111,14 @@ MnPicoEngine::~MnPicoEngine()
bool
MnPicoEngine::updateParValues(const std::vector<double>& vals){
double ns=vals[ins];
double logA=vals[iAs];
double mnu=vals[inu];
const double kout=0.05;
const double kin=0.002;
for (size_t i=0;i<_index.size();i++){
if (string(picoNames[i])=="scalar_amp(1)") {
double As=exp(logA)*1e-10;
double Asout=As/pow(kout/kin,ns-1);
PyDict_SetItemString(dict,picoNames[i],PyFloat_FromDouble(Asout));
PyDict_SetItemString(dict,picoNames[i],PyFloat_FromDouble(As));
}
else if (string(picoNames[i])=="omnuh2"){
PyDict_SetItemString(dict,picoNames[i],PyFloat_FromDouble(mnu/93.1));
......
......@@ -6,27 +6,7 @@
extern "C" {
#endif
PyObject* Py_Check(PyObject *result);
PyObject* pico_load(char *file);
PyObject* pico_compute_result(PyObject *pPico, int nparams, char *names[], double values[]);
PyObject* pico_compute_result_dict(PyObject *pPico, PyObject *pParams, PyObject *pOutputs);
void pico_get_output_len(PyObject *pResult, char *key, int *nresult);
void pico_read_output(PyObject *pResult, char *key, double **result, int *istart, int *iend);
bool pico_has_output(PyObject *pPico, char* output);
void pico_set_verbose(PyObject *pPico, bool verbose);
bool pico_is_verbose(PyObject *pPico);
void pico_free_result(PyObject *pResult);
#include"pico.h" // should be in ypur include path
#ifdef __cplusplus
}
......
......@@ -11,51 +11,34 @@ using namespace std;
int main(int argc,char** argv){
pico_init(true);
cout << "loading " << PICODATA << endl;
PyObject* pico=pico_load(PICODATA);
Py_Check(pico);
pico_set_verbose(pico,true);
char* name[]={ "ombh2","omch2","theta","re_optical_depth","scalar_spectral_index(1)","scalar_nrun(1)","scalar_amp(1)","omnuh2","massive_neutrinos","helium_fraction"};
const int nparams=sizeof(name)/sizeof(char*);
double* val=new double[nparams];
val[0]=0.2206934E-01;
val[1]=0.1202510;
val[2]=0.1041305E+01/100.;
val[3]=0.9273149E-01;
const double ns=0.9581583;
val[4]=ns;
val[5]=0;
const double logA=0.3095895E+01;
const double kout=0.05;
const double kin=0.002;
double As=exp(logA)*1e-10;
val[6]=As/pow(kout/kin,ns-1);
val[7]=0.06/93.1;
val[8]=3.046;
val[9]=0.248;
//for (size_t i=0;i<10;i++) cout << name[i] <<"=" << val[i] << endl;
//with dictionnary
//inputs
PyObject* params=PyDict_New();
PyDict_SetItemString(params,"force",Py_True);
for (size_t i=0;i<nparams;i++) PyDict_SetItemString(params,name[i],PyFloat_FromDouble(val[i]));
//test l'option force
//wanted outputs
char* outnames[]={"lensed_TT","lensed_EE","lensed_TE","lensed_BB"};
//add pivot scale
PyDict_SetItemString(params,"pivot_scalar",PyFloat_FromDouble(0.05));
PyDict_SetItemString(params,"ombh2",PyFloat_FromDouble(0.2206934E-01));
PyDict_SetItemString(params,"omch2",PyFloat_FromDouble(0.1202510E+00));
PyDict_SetItemString(params,"theta",PyFloat_FromDouble(0.1041305E+01/100.));
PyDict_SetItemString(params,"re_optical_depth",PyFloat_FromDouble(0.07));
PyDict_SetItemString(params,"scalar_spectral_index(1)",PyFloat_FromDouble(0.96));
PyDict_SetItemString(params,"scalar_nrun(1)",PyFloat_FromDouble(0.));
const double logA=3.02;
PyDict_SetItemString(params,"scalar_amp(1)",PyFloat_FromDouble(exp(logA)*1e-10));
const double mnu=0.06;
PyDict_SetItemString(params,"omnuh2",PyFloat_FromDouble(mnu/93));
PyDict_SetItemString(params,"massive_neutrinos",PyFloat_FromDouble(3.046));
PyDict_SetItemString(params,"helium_fraction",PyFloat_FromDouble(0.24));
PyDict_SetItemString(params,"Alens",PyFloat_FromDouble(1));
//outputs
char* outnames[]={"cl_TT","cl_EE","cl_TE","cl_BB"};
const int NOUT=sizeof(outnames)/sizeof(char*);
PyObject* output=PyList_New(NOUT);
......@@ -74,7 +57,7 @@ int main(int argc,char** argv){
double* ClBB=NULL;
int lmin=2;
int lmax=5000;
int lmax=10;
//check lower than nel
int N;
pico_get_output_len(pResult,outnames[0],&N);
......@@ -88,16 +71,12 @@ int main(int argc,char** argv){
pico_read_output(pResult, outnames[3],&ClBB,&zero,&lmax);
string filename("planck_wp_hlGrid.pico");
ofstream of(filename.c_str());
cout <<"writing " << filename << endl;
of << scientific << setprecision(12) ;
cout << scientific;
for (size_t l=lmin;l<=lmax;l++){
double fac=l*(l+1)/(2*M_PI);
of << l << "\t" << ClTT[l]/fac << "\t" << ClTE[l]/fac << "\t" << ClEE[l]/fac << "\t" << ClBB[l]/fac << endl;
cout << l << "\t" << ClTT[l]/fac << "\t" << ClTE[l]/fac << "\t" << ClEE[l]/fac << "\t" << ClBB[l]/fac << endl;
}
//clean
delete [] ClTT;
delete [] ClEE;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment