source: Sophya/trunk/Cosmo/RadioBeam/applobe.cc@ 3789

/*  ------------------------ Projet BAORadio -------------------- 
  Programme de convolution avec le lobe d'un cube 3D (angles,freq) 
    R. Ansari , C. Magneville - Juin 2010 

  Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamArcmin]
---------------------------------------------------------------  */

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

#include "array.h"
#include "histats.h"

#include "swfitsdtable.h"
#include "fitshdtable.h"

#include "randr48.h"
#include "vector3d.h"

// #include "xastropack.h"     -- Pour faire les conversions de coordonnees celestes

#include "radutil.h"
#include "lobe.h"

// Pour l'initialisation des modules 
#include "tarrinit.h"
#include "histinit.h"
#include "fiosinit.h"     

#include "timing.h"
#include "ctimer.h"

#include "cubedef.h"

//----------------------------------------------------------------------------
//----------------------------------------------------------------------------   
int main(int narg, char* arg[])
{
  // Sophya modules initialization
  TArrayInitiator  _inia;
  HiStatsInitiator  _inih;
  FitsIOServerInitiator  _inif;
  //------- AU LIEU DE ------>  SophyaInit();  

  InitTim();   // Initializing the CPU timer
  Timer tm("applobe");

  if (narg < 3) {
    cout << "Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamArcmin]\n" << endl;
    return 1;
  }

  // decodage arguments 
  string outname = arg[2];
  string inname = arg[1];
  bool fgresample=false;
  double fresamp=1.;
  if (narg>3) {
    fresamp=atof(arg[3]);
    if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2))  fgresample=true; 
  }
  bool fgcorrbeam=false;
  double tbeamarcmin=15.;
  if (narg>4) {
    tbeamarcmin=atof(arg[4]);
    if (tbeamarcmin>1.e-3)  fgcorrbeam=true;
  }

  int rc = 91;

  cout << " ====== applobe :   Input NVSS catalog name= " << inname << " OutName=" << outname;
  bool fginmap=true;
  try {
    TArray<r_4> incube;
    cout << "applobe[1]: reading input 3D map (cube) from file " << inname << endl;
    {
      PInPersist pin(inname);
      pin >> incube;
    }
    incube.Show(); 

    double dxdeg = ThetaSizeDegre/(double)NTheta;
    double dydeg = PhiSizeDegre/(double)NPhi;
    double dx = DegreeToRadian(dxdeg);
    double dy = DegreeToRadian(dydeg);
    double dfreq = FreqSizeMHz/(double)NFreq;

    cout << " X,Y map size in degrees , X/Phi=" << PhiSizeDegre << " Y/Theta=" << ThetaSizeDegre
         << " \n dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx << " dy_rad=" << dy << ")" 
         << " FreqSize=" << FreqSizeMHz << " dfreq=dz= " << dfreq << " MHz" << endl;

    double mean, sigma;
    MeanSigma(incube, mean, sigma);
    cout << " InCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;

    cout << "applobe[2]: creating Four2DResponse and BeamEffect objects..." << endl;
    H21Conversions conv;
    conv.setFrequency(Freq0MHz);
    double lambda = conv.getLambda();
    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda, lambda);    
    BeamEffect beam(fresp);

    if (fgcorrbeam) {
      double DoL = 1.22/ArcminToRadian(tbeamarcmin);
      Four2DResponse tbeam(2, DoL, DoL );
      cout << "applobe[3]: calling Correct2RefLobe() with target beam:" << tbeamarcmin 
           << " arcmin -> D/Lambda=" << DoL << endl;
      beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq);
    }
    else {
      cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
      beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
    }
    TArray< r_4 > outcube;
    if (fgresample) {
      cout << "applobe[4]: calling ReSample(incube," << fresamp << "," << ",1.) ... " << endl;
      outcube.Share(beam.ReSample(incube, fresamp, fresamp, 1.));
    }
    else outcube.Share(incube);

    outcube.Show(); 
    MeanSigma(outcube, mean, sigma);
    cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;

    // On sauve le cube de sortie
    {
      cout << " applobe[5]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << outcube;
    }

    rc = 0;
  }
  catch (PThrowable& exc) {
    cerr << " applobe.cc catched SOPHYA Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n applobe.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " applobe.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

  cout << ">>>> applobe[9] ------- FIN ----------- Rc=" << rc << endl;
  return rc;
}