source: Sophya/trunk/Cosmo/RadioBeam/pknoise.cc@ 4049

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de calcul du spectre de puissance de bruit pour 
  un interferometre (spectre moyenne sur 3D -> P_noise(k) )

    R. Ansari , C. Magneville - Juin 2010 

  Usage:  pknoise pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName
          -parname : -noise , -renmax, -scut, -ngen, -z, -bsize, -cszmpc, -nbin -prt 
---------------------------------------------------------------  */

#include "machdefs.h"
#include "sopnamsp.h"
#include <iostream>
#include <string>
#include <math.h>

#include <typeinfo>

#include "mdish.h"
#include "specpk.h"
#include "interfconfigs.h"
#include "radutil.h"

#include "histinit.h"
// #include "fiosinit.h"
// #include "fitsioserver.h"

#include "randr48.h"      

#include "timing.h"
#include "ctimer.h"
#include "slininterp.h"

typedef DR48RandGen RandomGenerator ;

// ---------------------------------------------------------------------
// Test main 
// R. Ansari  - Avril-Dec 2010 
// ---------------------------------------------------------------------

void Usage() 
{
  cout << " Usage:  pknoise [-parname value] Diameter/Four2DRespTableFile OutPPFName \n" 
       << "    -noise NoiseLevel (default=1.) \n"
       << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
       << "    -ngen NGen (default=0) number of noise fourier amp generations \n"
       << "       NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
       << "       NGen=-1 -> Call ComputeNoisePk(), dont change Four3DPk cell size \n"
       << "    -bsize sx,sy,sz : 3D real space box size (default=512x512x256) \n"   
       << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"   
       << "    -z redshift : redshift  (default=1.0) \n"
       << "    -pnf : Apply frequency dependent noise factor=ComovDa^2*(1+z)^2/Hz (default= NON) \n"
       << "    -scut SCutValue (default= -100.) \n"
       << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
       << "    -prt PrtLev,PrtModulo (default=0,10) \n"
       << "    -nbin PkNBin : nomber of histogram bins for Pnoise(k) (default= 256) \n"
       << "  NOTE : this program needs redshift_da.txt andredshift_Hz.txt files" << endl;
  return;
}

//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
int main(int narg, const char* arg[])
{
  if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) )  {
    Usage();
    return 1;
  }
  cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
  // make sure SOPHYA modules are initialized 
  SophyaInit();  
  //  FitsIOServerInit();
  InitTim();
  //--- decoding command line arguments 
  string tits="pknoise";
  char cbuff[64];
  bool fgresptbl=false;
  double DIAMETRE=100.;
  string resptblname;
  double NoiseLevel=1.;
  
  bool fgrenorm=false;
  double rmax=1.;
  int NMAX = 0;
  double SCut=-100.;
  bool fgautoscut=true;
  double FacSCut=-SCut;
  double z_Redshift=1.0 ;  // z=0.7 : 21 cm at z=0.7 -> 0.357 m  
  double comov_dA_z=3330.;  // Comoving radial distance = (1+z)dA
  double H_z=120.5;  // Hubble_param(z) 
  int box3dsz[3]={512,512,256};
  double cellsz[3]={3.,3.,3.};
  bool fgnoisefreq=false;
  int NBINPK=256;
  int prtlev=0;
  int prtmod=10;

  int ka=1;
  while (ka<(narg-1)) {
    if (strcmp(arg[ka],"-noise")==0) {
      NoiseLevel=atof(arg[ka+1]);
      ka+=2;
    }
    else if (strcmp(arg[ka],"-renmax")==0) {
      rmax=atof(arg[ka+1]);  fgrenorm=true;   ka+=2;
    }
    else if (strcmp(arg[ka],"-scut")==0) {
      SCut=atof(arg[ka+1]);    ka+=2; 
      if (SCut<0.) { FacSCut=-SCut;  fgautoscut=true; }
    }
    else if (strcmp(arg[ka],"-ngen")==0) {
      NMAX=atoi(arg[ka+1]);    ka+=2;
    }
    else if (strcmp(arg[ka],"-z")==0) {
      z_Redshift=atof(arg[ka+1]);  ka+=2;
    }
    else if (strcmp(arg[ka],"-bsize")==0) {
      sscanf(arg[ka+1],"%d,%d,%d",box3dsz,box3dsz+1,box3dsz+2);   ka+=2;
    }
    else if (strcmp(arg[ka],"-cszmpc")==0) {
      sscanf(arg[ka+1],"%lg,%lg,%lg",cellsz,cellsz+1,cellsz+2);   ka+=2;
    }
    else if (strcmp(arg[ka],"-nbin")==0) {
      NBINPK=atoi(arg[ka+1]);    ka+=2;
    }
    else if (strcmp(arg[ka],"-pnf")==0) {
      fgnoisefreq=true;   ka++;
    }
    else if (strcmp(arg[ka],"-prt")==0) {
      sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod);   ka+=2;
    }
    else break; 
  }

  if ((ka+1)>=narg) {
    cout << " pknoise / Argument error " << endl;
    Usage();
    return 2;
  }
  if (isdigit(*arg[ka])) {
    fgresptbl=false;
    DIAMETRE=atof(arg[ka]);
    sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
  }
  else { 
    resptblname=arg[ka];
    fgresptbl=true;
    sprintf(cbuff,"pknoise_RespTblName=%s", arg[ka]);
  }
  tits=cbuff;
  string outfile = arg[ka+1];  
  if (outfile==".")  outfile = "pknoise.ppf";
 //-- end command line arguments
  
  int rc = 1;  
  try {  // exception handling try bloc at top level
    cout << " pknoise[0] : Executing, output file= " << outfile << endl;  
    POutPersist po(outfile);
    
    cout << " pknoise[0.b] : Creating interpolating objects from redshift_da.txt and redshift_Hz.txt " << endl;
    SLinInterp1D redshift2da;
    string interpfilename="redshift_da.txt";
    redshift2da.ReadXYFromFile(interpfilename);
    SLinInterp1D redshift2hz;
    interpfilename="redshift_Hz.txt";
    redshift2hz.ReadXYFromFile(interpfilename);
    comov_dA_z=redshift2da(z_Redshift);
    H_z=redshift2hz(z_Redshift);
    cout << " Redshift= " << z_Redshift << " -> ComovdA= " << comov_dA_z
         << " Mpc  H(z)= " << H_z << " km/s/Mpc " << endl;

    H21Conversions conv;
    conv.setRedshift(z_Redshift);
    double lambda = conv.getLambda();

    double Dol = DIAMETRE/lambda;

    Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
    Four2DResponse* arep_p=&arep;
    Four2DRespTable resptbl;
    if (fgresptbl) {
      cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
      resptbl.readFromPPF(resptblname);
      cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda(" 
           << lambda << ")" << endl;
      resptbl.setLambda(lambda);
      arep_p=&resptbl;
      if (fgrenorm) {
        cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")"; 
        double omax=resptbl.renormalize(rmax);
        cout << " ... Old Max=" << omax << endl;
      }
    }
    else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    Histo2D h2drep = arep_p->GetResponse();
    Histo2D h2drepuv = arep_p->GetResponseUV();
    double repmax= h2drep.VMax();
    if (fgautoscut) {
      SCut = FacSCut/repmax;
      cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " --> SCut=" << FacSCut << "/repmax=" 
           << SCut << endl;
    }
    else cout << " pknoise[1.b]: Four2DResponse.RepMax=" << repmax << " , SCut=" << SCut << endl;
    
    cout << " pknoise[2]: Instanciating object type Four3DPk  " << endl;
    RandomGenerator rg;

    double dfreq=H_z*cellsz[2]/(1+z_Redshift)/lambda/1000.;
    double freq0=conv.getFrequency()-dfreq*box3dsz[2]/2;

    double dkxmpc=2.*DeuxPI; 
    double dkympc=2.*DeuxPI; 
    double dkzmpc=2.*DeuxPI; 
    double angscale=1.;
    Vector angscales(box3dsz[2]);
    double pnoise0=comov_dA_z*comov_dA_z*(1.+z_Redshift)*(1.+z_Redshift)/H_z;
    Vector pnoisevec(box3dsz[2]);
    pnoisevec=1.;
    angscales=angscale;
    if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[2]=128; }
    else {
      angscale=comov_dA_z;
      for(int kz=0; kz<box3dsz[2]; kz++) angscales(kz)=redshift2da(1420.4/(freq0+(double)kz*dfreq)-1.);
      dkxmpc = DeuxPI/(double)box3dsz[0]/cellsz[0];
      dkympc = DeuxPI/(double)box3dsz[1]/cellsz[1];
      dkzmpc = DeuxPI/(double)box3dsz[2]/cellsz[2];
    }
    cout << " pknoise[2.b]: ComovDist=" << comov_dA_z << " Mpc H(z)=" << H_z << " Mpc/km/s -> angscale=" << angscale << endl; 
    cout << " pknoise[2.c]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")=" 
         << conv.getFrequency() << " zmin=" << 1420.4/freq0-1. << " zmax=" << 1420.4/(freq0+box3dsz[2]*dfreq)-1. << endl;
    if(fgnoisefreq) {
      for(int kz=0; kz<box3dsz[2]; kz++) {
        double zreds=1420.4/(freq0+(double)kz*dfreq)-1.;
        double comda=redshift2da(zreds);
        pnoisevec(kz)=comda*comda*(1.+zreds)*(1.+zreds)/redshift2hz(zreds)/pnoise0;
      }
      cout << " pknoise[2.d]: PnoiseFactor:" << pnoisevec(0) << " ... " << pnoisevec(pnoisevec.Size()-1) << endl;
    }
    cout << " pknoise[2.e]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
         << box3dsz[2] << ")" << endl;
    Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
    cout << " pknoise[2.f]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
         << ") cell size (Mpc) : " << cellsz[0] << "x" << cellsz[1] << "x" << cellsz[2] << endl;
    m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
    m3d.SetPrtLevel(prtlev,prtmod);

    if (NMAX>0) {
      cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
      PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
      pkn.SetFreqRange(freq0, dfreq);
      pkn.SetAngScaleConversion(angscales);
      pkn.SetPNoiseFactor(pnoisevec);
      pkn.SetPrtLevel(prtlev,prtmod);
      HProf hpn = pkn.Compute(NBINPK);
    }
    else {
      cout << " pknoise[3]: Computing Noise P(k) m3d.ComputeNoisePk(...) " << endl;
      HProf hpn = m3d.ComputeNoisePk(*(arep_p),angscale,SCut,NBINPK);
    }
        
    HProf hpnoise=m3d.GetPk();
    DataTable dtnoise;
    Histo fracmodok=m3d.FillPkDataTable(dtnoise);
    HProf h1dnoise=arep_p->GetProjNoiseLevel();
    HProf h1drep=arep_p->GetProjResponse();
    cout << " pknoise[3.b]: writing dtnoise,hpn,h2rep... with tags to " << outfile << endl;
    po << PPFNameTag("dtnoise") << dtnoise;
    po << PPFNameTag("hpnoise") << hpnoise;
    po << PPFNameTag("fracmodok") << fracmodok;
    po << PPFNameTag("h1dnoise") << h1dnoise;
    po << PPFNameTag("h1drep") << h1drep;
    po << PPFNameTag("h2drep") << h2drep;
    po << PPFNameTag("h2drepuv") << h2drepuv;

    rc = 0;
  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << "\n...exc.Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {  // catching standard C++ exceptions
    cerr << " pknoise.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  PrtTim("End-pknoise");
  cout << " ==== End of pknoise.cc program  Rc= " << rc << endl;
  return rc;    
}