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

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

  Usage: applobe [-g -t -fib]  Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName 
                 [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] 
   o -t : triangle beam shape in k-space
   o -g : gaussian beam shape in k-space 
   o -fib: application d'un lobe fixe (independant de la frequence) 
   o Diameter/Four2DRespTableFile : diametre dish ou nom de fichier PPF reponse 2D
   o In3DPPFName Out3DPPFName : Noms des fichiers (TArray<r_4> 3D) entree-sortie 
   o TargetBeamDiam : Correction à un beam fiduciel, independant de la frequence 
     avec specification de la valeur du diametre pour la frequence la plus basse
   o ResmapleFactor : reechantillonnage du cube (2 -> surechantillonnage, 0.5 : 1 sur 2)
---------------------------------------------------------------  */

#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

  if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
    cout << "Usage: applobe [-t -g -fib -kzf/-nzf -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n" 
         << "               [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
      cout << "   -t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
           << "   -fib :  Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n" 
           << "   -kzf -nzf :  Keep (default) or Not zero space frequency when applying lobes (BeamEffect class) \n" 
           << "   -mxr val: Max beam correction factor (default=10.) \n"
           << "   Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n" 
           << "   In3DPPFName Out3DPPFName: Input/Output PPF file names  (TArray<r_4> 3D) \n" 
           << "   TargetBeamDiam: Corrected beam target diameter (D/Lambda at lowest frequency) \n" 
           << "      DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm)) \n"
           << "   ResmapleFactor : Resampling (2-> oversampling, 0.5 : 1/2 undersampling) \n" << endl;    
    }  
    return 1;
  }

  Timer tm("applobe");

  // decodage arguments 
  bool fggaussian=true;  // true -> gaussian beam 
  bool fixedbeam=false;  // true -> apply freq. independent beam 
  double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe 
  bool preservefreq0=true; // true -> keep zero frequency when appyling lobe 

  // decodage argument optionnel 
  bool fgoptarg=true;
  while (fgoptarg) {
    string fbo = arg[1];
    if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
    else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
    else if (fbo=="-fib")  { fixedbeam=true; arg++; narg--; }
    else if (fbo=="-kzf")  { preservefreq0=true; arg++; narg--; }
    else if (fbo=="-nzf")  { preservefreq0=false; arg++; narg--; }
    else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
    else fgoptarg=false;
  }
  if (narg < 4) {
    cout << " applobe/error arguments , applobe -h for help " << endl;
    return 2;
  }

  bool fgresptbl=true;
  double DIAMETRE=100.;
  string resptblname;
  if (isdigit(*arg[1])) {
    fgresptbl=false;
    DIAMETRE=atof(arg[1]);
  }
  else  resptblname=arg[1];
  
  string inname = arg[2];
  string outname = arg[3];
  bool fgcorrbeam=false;
  double tbeamDiam=50.;
  if (narg>4) {
    tbeamDiam=atof(arg[4]);
    if (tbeamDiam>1.)  fgcorrbeam=true;
  }
  bool fgresample=false;
  double fresamp=1.;
  if (narg>5) {
    fresamp=atof(arg[5]);
    if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2))  fgresample=true; 
  }

  int rc = 91;

  cout << " ====== applobe :   Input skycube  name= " << inname << " OutName=" << outname << endl;
  cout << ((fggaussian)?" Gaussian ":" Triangular") << ((fixedbeam)?" Fixed (freq.independent)":"") << " beams" << endl;
  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();
    
    int typbeam=(fggaussian)?1:2;
    Four2DResponse fresp(typbeam, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
    Four2DResponse* fresp_p=&fresp;
    Four2DRespTable resptbl;
    if (fgresptbl) {
      cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
      resptbl.readFromPPF(resptblname);
      resptbl.renormalize(1.);
      fresp_p=&resptbl;
    }
    else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    BeamEffect beam(*fresp_p,preservefreq0);

    if (fgcorrbeam) {
      double DoL = tbeamDiam/lambda;
      double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
      int typcb = (fggaussian)?1:2;
      //      if (fgresptbl) typcb=22;
      Four2DResponse tbeam(typcb, DoL, DoL );
      cout << "applobe[3]: calling Correct2RefLobe() with target beam Diameter=" << tbeamDiam 
           << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb 
           << " MaxRatio=" << maxratio << endl;
      beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq, maxratio);
    }
    else if (fixedbeam) {
      cout << "applobe[3]: calling ApplyLobe() (freq.independent beam) ... " << endl;
      beam.ApplyLobe(incube, dx, dy, Freq0MHz);
    }
    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;
}