source: Sophya/trunk/Cosmo/RadioBeam/calcpk2.cc@ 3784

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de calcul du spectre de puissance (3D) a partir d'un 
  cube de  delta T/T LSS, d'un cube  delta T/T LSS synchrotron
  ou radio-sources, apres ajustement / soustraction d'une loi de 
  puissance en frequence (l'axe Z du tableau doit etre en frequence)

    R. Ansari , C. Magneville - Juin 2010 

Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk 
---------------------------------------------------------------  */

#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 "histinit.h"
#include "fftwserver.h"
#include "randr48.h"      

#include "ctimer.h"

typedef ThSDR48RandGen RandomGenerator ;
//--- Declaration des fonctions 
TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx);

//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
/* --Fonction-- */
int main(int narg, const char* arg[])
{
  if (narg<6) {
    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk  " << endl;
    return 1;
  }
  Timer tm("calcpk2");
  int rc = 0;
  try {
    string inppflss = arg[1];
    r_4 rfaclss = atof(arg[2]);
    string inppfsync = arg[3];
    r_4 rfacsync = atof(arg[4]);
    string outname = arg[5];
    TArray<r_4> inmaplss, inmapsync;
    const char * tits[2]={"LSS", "Sync"};
    for(int ks=0; ks<2; ks++)  {
      string& ppfname=inppflss;
      r_4 rfac=rfaclss;
      TArray<r_4>* inmap=&inmaplss; 
      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&inmapsync; }
      cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname 
           << "  RenormFactor=" << rfac << endl;
      PInPersist pin(ppfname);
      pin >> (*inmap);
      (*inmap) *= rfac;
      double mean, sigma;
      MeanSigma(*inmap, mean, sigma);
      cout << " ...InMap sizes " << inmap->InfoString() << endl;
      inmap->Show(); 
      cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
    }
    tm.Split(" After read ");
    TArray<r_8> synctemp, specidx; 
    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, synctemp, specidx);
    double mean, sigma;
    MeanSigma(inmap, mean, sigma);
    cout << " After cleaning: Mean=" << mean << " Sigma=" << sigma << endl;      
    tm.Split(" After CleanForeground");

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

    cout << "calcpk2[5] : computing power spectrum ... " << endl;
    RandomGenerator rg;
    Four3DPk pkc(four3d, rg);
    
    HProf hp = pkc.ComputePk(0.,256);

    tm.Split(" Done ComputePk ");    
    
    cout << "calcpk2[4] : writing profile P(k) and foreground maps to  " << outname << endl;
    POutPersist po(outname);
    po << PPFNameTag("Pk") << hp;
    po << PPFNameTag("Tsync") << synctemp;
    po << PPFNameTag("async") << specidx;

  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " calcpk2.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 << " calcpk2.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " calcpk2.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
  return rc;    
}


/* --Fonction-- */
TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx)
{
  bool smo;
  if (!maplss.CompareSizes(mapsync,smo) ) {
    cout << " CleanForeground/ERROR sizes " << endl;
    maplss.Show();  mapsync.Show(); 
    throw ParmError("CleanForeground- maplss , mapsync not the same size !");
  }
  sa_size_t sz[5];   sz[0]=maplss.SizeX();  sz[1]=maplss.SizeY();
  synctemp.SetSize(2, sz); 
  specidx.SetSize(2, sz);
  TArray<r_4>& omap=maplss;
  Vector vlnf(maplss.SizeZ());
  int nprt = 0;
  double freq0=840.;  // Frequence premier index en k (MHz)
  double dfreq=1.;   // largeur en frequence de chaque plan (Mhz)  
  for(sa_size_t i=0; i<maplss.SizeX(); i++) 
    for(sa_size_t j=0; j<maplss.SizeY(); j++)  {
      r_8 s1, sx, sx2, sy, sxy;
      s1=sx=sx2=sy=sxy=0.;
      for(sa_size_t k=0; k<maplss.SizeZ(); k++)  {
        double lnf=log((double)k*dfreq+freq0);
        vlnf(k)=lnf;
        double ttot=(r_8)(mapsync(i,j,k)+maplss(i,j,k));
        if (ttot < 1.e-5) continue;
        double lntt=log(ttot);
        s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
        sy+=lntt;    sxy+=(lnf*lntt);
      }
      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
      double T0 = exp(beta+alpha*vlnf(0));
      if ((i%6==0)&&(j%7==0)) 
        cout << "CleanForeground[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
             << " (mapsync=" << mapsync(i,j,0) << " ... " << mapsync(i,j,125) << ")" << endl; 
      synctemp(i,j) = T0;
      specidx(i,j) = alpha;
      for(sa_size_t k=0; k<maplss.SizeZ(); k++) {
        r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
        omap(i,j,k) = mapsync(i,j,k)+maplss(i,j,k)-fittedtemp;
      }
    }
  return omap;
}