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

/*  ------------------------ 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 , -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"

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"
       << "    -z redshift (default=0.7) \n"
       << "    -scut SCutValue (default= -100.) \n"
       << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
       << "    -prt PrtLev,PrtModulo (default=0,10) " << 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=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
  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],"-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);

    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();
    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;
    Four3DPk m3d(rg);
    m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI); 
    m3d.SetPrtLevel(prtlev,prtmod);

    cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    if (NMAX>0) {
      PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
      pkn.SetPrtLevel(prtlev,prtmod);
      HProf hpn = pkn.Compute();
      cout << " pknoise[3.b]: writing hpn noise profile to " << outfile << endl;
      po << hpn;
    }
    else {
      Histo fracmodok;
      DataTable dtnoise;
      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,SCut);
      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") << hpn;
      po << PPFNameTag("fracmodok") << fracmodok;
      po << PPFNameTag("h1dnoise") << h1dnoise;
      po << PPFNameTag("h1drep") << h1drep;
      po << PPFNameTag("h2drep") << h2drep;
    }
    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;    
}