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

/*  ------------------------ 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 Diameter/Four2DRespTableFile OutPPFName 
                 [NGen=1] [S2Cut=100.] [z_redshift=0.7]
---------------------------------------------------------------  */

#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 
// ---------------------------------------------------------------------



//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
int main(int narg, const char* arg[])
{
  if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) )  {
    cout << " Usage:  pknoise Diameter/Four2DRespTableFile OutPPFName [RenormalizeMax]\n" 
         << "         [NGen=1] [S2Cut=100.] [z_redshift=0.7] [PrtLev,PrtModulo=0,10]" << endl;
    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;
  if (isdigit(*arg[1])) {
    fgresptbl=false;
    DIAMETRE=atof(arg[1]);
    sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
  }
  else { 
    resptblname=arg[1];
    fgresptbl=true;
    sprintf(cbuff,"pknoise_RespTblName=%s", arg[1]);
  }
  tits=cbuff;
  string outfile = arg[2];  
  if (outfile==".")  outfile = "pknoise.ppf";
  bool fgrenorm=false;
  double rmax=1.;
  if (narg>3) {
    rmax=atof(arg[3]);
    fgrenorm=true;
  }
  int NMAX = 1;
  if (narg>4) NMAX = atoi(arg[4]);
  double SCut=0.;
  if (narg>5) SCut = atof(arg[5]);
  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
  if (narg>6) z_Redshift = atof(arg[6]);
  int prtlev=0;
  int prtmod=10;
  if (narg>7) sscanf(arg[7],"%d,%d",&prtlev,&prtmod);
  //-- 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);
      arep_p=&resptbl;
      if (fgrenorm) {
        cout << " pknoise[1.b] 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;

    
    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;
    PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
    pkn.SetPrtLevel(prtlev,prtmod);
    HProf hpn = pkn.Compute();
    po << hpn;
 
    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;    
}