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

//  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 "qhist.h" 

#ifndef DeuxPI
#define DeuxPI 2.*M_PI 
#endif

// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (frequences angulaires 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) 
// typ=22 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda  avec un trou au centre  

class Four2DResponse  {
public:
  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
  Four2DResponse(int typ, double dx, double dy, double lambda=1.);

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

  inline void setLambdaRef(double lambda=1.)
  { lambdaref_ = lambda; }
  inline void setLambda(double lambda=1.)
  { lambda_ = lambda;   lambda_ratio_ = lambda_/lambdaref_; }   
  
  // 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=255, int ny=255);  
  inline double D() { return dx_; } ; 
  inline double Dx() { return dx_; } ; 
  inline double Dy() { return dy_; } ; 

  int typ_;
  double dx_, dy_; 
  double lambdaref_, lambda_;
  double lambda_ratio_;     // lambdaref_/lambda_;

};


// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
class Four2DRespTable : public  Four2DResponse {
public:
  // Constructeur sans argument, utilise pour lire depuis un fichier 
  Four2DRespTable();
  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
  Four2DRespTable(Histo2D const & hrep, double d, double lambda=1.);
  // Apres renormalisaton Value(kx,ky) <= max
  double renormalize(double max=1.);
  // Return the 2D response for wave vector (kx,ky)
  virtual double Value(double kx, double ky); 

  void writeToPPF(string flnm);
  void readFromPPF(string flnm);

  Histo2D hrep_;
};


// -- Four2DRespRatio: Retourne le rapport de la reponse entre deux objets Four2DResponse
class Four2DRespRatio : public  Four2DResponse {
public:
  Four2DRespRatio(Four2DResponse& a, Four2DResponse& b, double divzthr=5.e-2);
  // Return the ratio a.Value(kx,ky) / b.Value(kx, ky) - with protection against divide by zero 
  virtual double Value(double kx, double ky); 
  Four2DResponse& a_;
  Four2DResponse& b_;
  double divzthr_;
};

// 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 Dy; }

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


// -------------------------------------------------------------------
// -- 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);

  // Pour phi, les angles phi, -phi, phi+pi, -(phi+pi) sont prises en compte
  inline void SetThetaPhiRange(double thetamax=0., int ntet=1, double phimax=0., int nphi=1)
  { thetamax_=thetamax; ntet_=ntet; phimax_=phimax; nphi_=nphi; }

  inline int SetPrtLevel(int lev=0, int prtmod=10) 
  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }

  inline void SetRespHisNBins(int nx=128, int ny=128)
  { nx_=nx; ny_=ny; }
  inline void SetBeamNSamples(int nx=128, int ny=128)
  { beamnx_=nx; beamny_=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_;   // nb de bins de l'histo de reponse 
  int beamnx_, beamny_; // nb de points d'echantillonnage du beam

  //   Histo2D h2w_, h2cnt_;
  QHis2D h2w_;
  int mcnt_;
  int prtlev_,prtmodulo_;
};


#endif