source: Sophya/trunk/Cosmo/RadioBeam/mdish.h@ 3771

//  Classes to compute Multi-Dish or CRT-like radio interferometer response
//  R. Ansari - Avril-Mai 2010 

#ifndef MDISH_SEEN
#define MDISH_SEEN

#include "machdefs.h"      // SOPHYA .h  
#include "sopnamsp.h"      // SOPHYA .h  
#include <math.h>
#include <iostream>
#include <vector>
#include <string> 

#include "genericfunc.h"   // SOPHYA .h  
#include "array.h"         // SOPHYA .h  
#include "histats.h"       // SOPHYA .h 

#define DeuxPI 2.*M_PI 

// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (angles theta,phi)
// typ=1 : Reponse gaussienne parabole diametre D exp[ - 0.5  (lambda  k_g / D )^2 ]
// typ=2 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda   
// typ=3 : Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda) 
class Four2DResponse  {
public:
  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
  Four2DResponse(int typ, double dx, double dy);

  Four2DResponse(Four2DResponse const& a)
  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; }
  Four2DResponse& operator=(Four2DResponse const& a)
  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }

  // Return the 2D response for wave vector (kx,ky)
  virtual double Value(double kx, double ky); 
  inline  double operator()(double kx, double ky) 
  { return Value(kx, ky); } 
  virtual Histo2D GetResponse(int nx=256, int ny=256);  
  inline double D() { return dx_; } ; 
  inline double Dx() { return dx_; } ; 
  inline double Dy() { return dy_; } ; 

  int typ_;
  double dx_, dy_; 
};

// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
class Four2DRespTable : public  Four2DResponse {
public:
  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
  Four2DRespTable(Histo2D const & hrep, double d);
  // Return the 2D response for wave vector (kx,ky)
  virtual double Value(double kx, double ky); 
  Histo2D hrep_;
};

// Classe toute simple pour representer un element de reception de type dish
class Dish {
public:
  // Circular dish 
  Dish(int id, double x, double y, double diam)
    :  id_(id), X(x), Y(y), D(diam), Dx(D), Dy(D), fgcirc_(true)   {   } 
  // Receiver with rectangular type answer in kx,ky plane 
  Dish(int id, double x, double y, double dx, double dy)
    :  id_(id), X(x), Y(y), D(sqrt(dx*dy)), Dx(dx), Dy(dy), fgcirc_(false)   {   } 
  Dish(Dish const& a) 
    :  id_(a.id_), X(a.X), Y(a.Y), D(a.D), Dx(a.Dx), Dy(a.Dy), fgcirc_(a.fgcirc_)     {   } 
  inline bool isCircular() { return fgcirc_; }
  inline int ReflectorId() { return id_; }
  inline double Diameter() { return D; }
  inline double DiameterX() { return Dx; }
  inline double DiameterY() { return Dx; }

  int id_;   // numero de reflecteur
  double D,X,Y;
  double Dx, Dy;
  bool fgcirc_;  // false -> rectangular dish 
};

// -----------------------------------
// -- Classe ressemblant a un histo 2D 
class QHis2D {
public:
  QHis2D();
  QHis2D(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
  void Define(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
  void SetZeroBin(r_8 x=0, r_8 y=0);
  inline double WBinX() { return dxb; }
  inline double WBinY() { return dyb; }
  Histo2D Convert();
  
  r_8 xmin,xmax,ymin,ymax;
  r_8 dxb,dyb;
  sa_size_t nx,ny;
  sa_size_t ixb0, jyb0;
  TArray<r_8> aw;
  double sumw0;
};

// -------------------------------------------------------------------
// -- Pour calculer la reponse ds le plan kx,ky d'un system MultiDish 
class MultiDish {
public:
  MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto=false);

  inline void SetThetaPhiRange(double thetamax=0., int ntet=1, double phimax=0., int nphi=1)
  { thetamax_=thetamax; ntet_=ntet; phimax_=phimax; nphi_=nphi; }

  inline void SetRespHisNBins(int nx=128, int ny=128)
  { nx_=nx; ny_=ny; }
  Histo2D GetResponse();

  double CumulResp(Four2DResponse& rd, double theta=0., double phi=0.);
  inline size_t NbDishes() { return dishes_.size(); }
  inline Dish&  operator[](size_t k) { return dishes_[k]; }

  virtual Histo2D PosDist(int nx=30, int ny=30, double dmax=0.);  

protected:
  double AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh);

  double lambda_, dmax_;
  vector<Dish> dishes_;
  bool fgnoauto_;
  double thetamax_, phimax_; 
  int ntet_,nphi_; 
  int nx_, ny_;
  //   Histo2D h2w_, h2cnt_;
  QHis2D h2w_;
  int mcnt_;
};


#endif