source: Sophya/trunk/Cosmo/RadioBeam/calcpk.cc@ 3825

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de calcul du spectre de puissance (3D) a partir d'un 
  cube de donnees (delta rho/rho ou delta T/T)
    R. Ansari , C. Magneville - Juin 2010 

  Usage: calcpk In3DMap OutPk  [InRenormFactor] 
---------------------------------------------------------------  */

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

#include <typeinfo>

#include "specpk.h"
#include "histats.h"

#include "qhist.h" 
#include "lobe.h" 
#include "cubedef.h" 

#include "histinit.h"
#include "fftwserver.h"
#include "randr48.h"      

#include "ctimer.h"

typedef ThSDR48RandGen RandomGenerator ;

//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
int main(int narg, const char* arg[])
{
  if (narg<3) {
    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
    return 1;
  }
  Timer tm("calcpk");
  int rc = 0;
  try {
    string inppfname = arg[1];
    string outname = arg[2];
    double renfact=1.;
    bool fgrenfact=false;
    if (narg>3) {
      renfact=atof(arg[3]);
      fgrenfact=true;
    }
    double pixsignoise = 0.;
    bool fgaddnoise=false;
    if (narg>4) {
      pixsignoise=atof(arg[4]);
      fgaddnoise=true;
    }

    cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
    TArray<r_4> inmap;
    {
      PInPersist pin(inppfname);
      pin >> inmap;
    }
    if (fgrenfact) { 
      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
      inmap *= renfact;
    }
    double mean, sigma;
    MeanSigma(inmap, mean, sigma);
    cout << " InMap sizes " << inmap.InfoString() << endl;
    inmap.Show(); 
    cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
    if (fgaddnoise) { 
      BeamEffect::AddNoise(inmap, pixsignoise);
    }

    tm.Split(" After read ");
    
    cout << "calcpk[2] : computing 3D Fourier coefficients ... " << endl;
    FFTWServer ffts;
    TArray< complex<r_4> > four3d;
    ffts.FFTForward(inmap, four3d);
    tm.Split(" After FFTForward ");

    cout << "calcpk[3] : computing power spectrum ... " << endl;
    RandomGenerator rg;
    Four3DPk pkc(four3d, rg);
    double dkxmpc = DeuxPI/(double)inmap.SizeX()/XCellSizeMpc;
    double dkympc = DeuxPI/(double)inmap.SizeY()/YCellSizeMpc;
    double dkzmpc = DeuxPI/(double)inmap.SizeZ()/ZCellSizeMpc;
    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);

    HProf hp = pkc.ComputePk(0.,256);
    {
      cout << "calcpk[4] : writing profile P(k) to PPF file  " << outname << endl;
      POutPersist po(outname);
      po << hp;
      outname = "f3d_" + outname;
      cout << "calcpk[4.b] : writing four3d to PPF file  " << outname << endl;
      POutPersist pof(outname);
      pof << four3d;
    }
    tm.Split(" Done CalcP(k) ");    
  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " calcpk.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 << " calcpk.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " calcpk.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
  return rc;    
}