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

  Usage: applobe In3DPPFName Out3DPPFName
---------------------------------------------------------------  */

#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 \n" << endl;
    return 1;
  }

  // decodage arguments 
  string outname = arg[2];
  string inname = arg[1];
  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;
    }
    cout << incube;

    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;
    // Dish de 64 m de diametre
    double Diametre = 64.;
    H21Conversions conv;
    conv.setFrequency(Freq0MHz);
    double lambda = conv.getLambda();
    Four2DResponse fresp(2, Diametre/lambda, Diametre/lambda);
    BeamEffect beam(fresp);

    cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
    beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);

    cout << "applobe[4]: calling ReSample(incube, 0.25, 0.25, 1.) ... " << endl;
    TArray< r_4 > outcube = beam.ReSample(incube, 0.25, 0.25, 1.);
    
    cout << outcube;
    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;
}