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

/*  ------------------------ 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 InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
---------------------------------------------------------------  */

#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 ;
//--- Declaration des fonctions 
TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
                            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 [PixNoiseLevel] " << 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];

    double pixsignoise = 0.;
    bool fgaddnoise=false;
    if (narg>6) {
      pixsignoise=atof(arg[6]);
      fgaddnoise=true;
    }


    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;
    }

    if (fgaddnoise) { 
      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
      BeamEffect::AddNoise(inmaplss, pixsignoise);
    }

    tm.Split(" After read ");
    TArray<r_8> synctemp, specidx; 
    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
    double freq0 = Freq0MHz;
    double dfreq = FreqSizeMHz/(double)NFreq;
    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, freq0, dfreq, 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, double freq0, double dfreq, 
                            TArray<r_8>& synctemp, TArray<r_8>& specidx)
// Inputs : maplss, mapsyc, freq0, dfreq
// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
// Return_Array : foreground subtracted LSS signal 
{
  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 : Frequence premier index en k (MHz)
  // double dfreq :   // 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;
}