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

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

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

#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 Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName " 
         << "             [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] \n" << endl;
    return 1;
  }

  // decodage arguments 
  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 tbeamDoL=135;
  if (narg>4) {
    tbeamDoL=atof(arg[4]);
    if (tbeamDoL>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;
  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, 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);
      fresp_p=&resptbl;
    }
    else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    BeamEffect beam(*fresp_p);

    if (fgcorrbeam) {
      double DoL = tbeamDoL;
      double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
      int typcb = 2;
      //      if (fgresptbl) typcb=22;
      Four2DResponse tbeam(typcb, DoL, DoL );
      cout << "applobe[3]: calling Correct2RefLobe() with target beam D/Lambda=" << DoL  
           << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << 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;
}