source: Sophya/trunk/Cosmo/SimLSS/cmvtransf.cc@ 3892

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "ntuple.h"

#include "constcosmo.h"
#include "pkspectrum.h"

void usage(void);
void usage(void) {
cout
<<"cmvtransf [options]"<<endl
<<"  -k npt,k1,k2 : k range in Mpc^-1"<<endl
<<"  -U h100,Om0,Ob0,tcmb : cosmology"<<endl
<<"  -F filename : read also CAMB file"<<endl
<<endl;
}

int main(int narg,char *arg[])
{
 double k1 = 1e-6, k2 = 10.; int_4 npt = 100;
 double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
 string fcmbfile = "";

 char c;
  while((c = getopt(narg,arg,"hk:U:F:")) != -1) {
  switch (c) {
  case 'k' :
    sscanf(optarg,"%d,%lf,%lf",&npt,&k1,&k2);
    break;
  case 'U' :
    sscanf(optarg,"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
    break;
  case 'F' :
    fcmbfile = optarg;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }

 if(k1<=0.) {cout<<"k1 must be >0"<<endl; return -1;}
 if(npt<=0) npt = 1000;
 if(k2<=k1) {k1=1.e-4; k2=1.e+2;}
 cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
 cout<<"fcmbfile="<<fcmbfile<<endl;

 cout<<"TransfertEisenstein with baryon"<<endl;
 TransfertEisenstein tf(h100,Om0-Ob0,Ob0,tcmb);
 cout<<"KPeak = "<<tf.KPeak()<<"  1/Mpc"<<endl;

 cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.1"<<endl;
 TransfertEisenstein tfnosc1(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.1
 tfnosc1.SetNoOscEnv(1);

 cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.2"<<endl;
 TransfertEisenstein tfnosc2(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.2
 tfnosc2.SetNoOscEnv(2);

 cout<<endl<<"TransfertEisenstein without baryon"<<endl;
 TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons

 TransfertTabulate tf_cmbfile;
 bool cmbfileOK = false;
 if(fcmbfile.size()>0) {
   cout<<endl<<"TransfertTabulate for CAMB"<<endl;
   tf_cmbfile.ReadCAMB(fcmbfile,h100);
   tf_cmbfile.SetInterpTyp(1);
   if(tf_cmbfile.NPoints()>0) cmbfileOK = true;
 }

 const int n = 6;
 const char *vname[n] = {"k","t","tnosc1","tnosc2","tnob","tcmbf"};
 NTuple nt(n,vname);
 double xnt[n];
 for(int i=0;i<n;i++) xnt[i]=0.;

 double lnk1 = log10(k1), lnk2=log10(k2), dlnk=(lnk2-lnk1)/npt;
 for(double lnk=lnk1;lnk<lnk2+dlnk/2.;lnk+=dlnk) {
   double k = pow(10.,lnk);
   xnt[0]  = k;
   xnt[1]  = tf(k);
   xnt[2]  = tfnosc1(k);
   xnt[3]  = tfnosc2(k);
   xnt[4]  = tfnob(k);
   if(cmbfileOK) xnt[5]  = tf_cmbfile(k);
   nt.Fill(xnt);
 }

 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvtransf.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);

 return 0;
}

/*
openppf cmvtransf.ppf

# Eisenstein
n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
n/plot nt.tnob%k tnob>0&&k>0 ! "nsta red same connectpoints logx logy"
n/plot nt.tnosc2%k tnosc2>0&&k>0 ! "nsta blue same connectpoints logx logy"
n/plot nt.tnosc1%k tnosc1>0&&k>0 ! "nsta green same connectpoints logx logy"

n/plot nt.t/tnob%k   tnob>0   ! "nsta red connectpoints"
n/plot nt.t/tnosc2%k tnosc2>0 ! "nsta blue same connectpoints"
n/plot nt.t/tnosc1%k tnosc1>0 ! "nsta green same connectpoints"

# CAMB
n/plot nt.t%k t>0&&k>0 ! "nsta connectpoints logx logy"
n/plot nt.tcmbf%k tcmbf>0&&k>0 ! "nsta red same connectpoints logx logy"

n/plot nt.(tcmbf-t)%k tcmbf>0. ! "nsta connectpoints"
n/plot nt.(tcmbf-t)/tcmbf%k tcmbf>0. ! "nsta connectpoints"

n/plot nt.t/tnosc2%k tnosc2>0&&t>0&&k>0 ! "nsta connectpoints logx logy"
n/plot nt.tcmbf/tnosc2%k tnosc2>0&&tcmbf>0&&k>0 ! "nsta connectpoints same red logx logy"
*/