Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4

acbeam.h

+/*  PAON4 analysis software 
+    Class representing an autocorrelation beam (adapted from JSkyMap code) 
+      R. Ansari, December 2018                                             */
+
+#ifndef ACBEAM_H_SEEN
+#define ACBEAM_H_SEEN
+
+#include "skymap.h"
+#include "unitvector.h"
+
+using namespace std; 
+using namespace SOPHYA; 
+
+//-----------------------------------------------------------------
+//---------- BeamTP (Theta,Phi) class 
+// represents a single dish beam response 
+//-----------------------------------------------------------------
+class ACBeam  {  //  Auto-correlation beam 
+public:
+  // constructor: define a lobe  to a circular dish of diameter D pointing to the direction (theta_b,phi_b)
+  ACBeam(double D, double theta_b=0., double phi_b=0., double lambda=0.21)
+    : fggauss_(false), D_(D), lambda_(lambda), DoL_(D/lambda), 
+      bdir_(LongitudeLatitude(theta_b, phi_b)), bdir_theta_(theta_b), bdir_phi_(phi_b), NormFac_(1.)
+  {
+  }
+  // copy constructor
+  ACBeam(ACBeam const& a)
+    : fggauss_(a.fggauss_), D_(a.D_), lambda_(a.lambda_), DoL_(a.DoL_), 
+      bdir_(a.bdir_), bdir_theta_(a.bdir_theta_), bdir_phi_(a.bdir_phi_), NormFac_(a.NormFac_)
+  {
+  }
+
+  // copy operator 
+  ACBeam& operator = (ACBeam const& a)
+  {
+    fggauss_=a.fggauss_;   D_=a.D_;     lambda_=a.lambda_;   DoL_=a.DoL_;  
+    bdir_=a.bdir_;  bdir_theta_=a.bdir_theta_; bdir_phi_=a.bdir_phi_;  NormFac_=a.NormFac_;
+    return (*this);
+  }
+  // defining a gaussian beam 
+  inline void setGaussianLobe(bool fgg=true)
+  {  fggauss_=fgg; }
+  // To define a rectangular feed response (for cylinder feeds) 
+  // set/changes the beam axis direction 
+  inline void setBeamAxisDirection(double theta_b, double phi_b)
+  {
+    bdir_ = UnitVector(LongitudeLatitude(theta_b, phi_b));
+    bdir_theta_=theta_b; 
+    bdir_phi_=phi_b;
+  }
+  // set/changes the beam axis direction 
+  inline void setBeamAxisDirection(UnitVector const& beamdir)
+  {
+    bdir_ = beamdir;
+    bdir_theta_=bdir_.GetThetaPhi().Theta(); 
+    bdir_phi_=bdir_.GetThetaPhi().Phi(); 
+  }
+  // return the beam axis direction 
+  inline UnitVector getBeamAxisDirection() const 
+  {
+    return bdir_; 
+  }
+  // Define the normalization factor 
+  inline void setBeamNormFactor(double fac=1.)
+  {
+    NormFac_=fac;
+  }
+  // return the beam normalization factor 
+  inline double getBeamNormFactor()  const
+  {
+    return NormFac_;
+  }
+  // Compute the beam normalization factor such as the beam integral over the sphere is equal 1 (unity)
+  void AutoNormalizeBeam(int m=128, int prtlev=0)
+  {
+    //Timing    Timer tm("BeamTP::AutoNormalizeBeam()");  tm.Nop();
+    double sum=0.;
+    SphereHEALPix<double> sphb=getBeamSph(sum,m);
+    if (sum>1.e-12) {
+      NormFac_=((double)sphb.NbPixels())/sum/4./M_PI;
+      if (prtlev>0)
+	cout << " BeamTP::AutoNormalizeBeam() sum="<<sum<<" NPixels="<<sphb.NbPixels()<<" -> NormFac=N/(sum*4Pi)="<<NormFac_<<endl;
+    }
+    else {
+      NormFac_=1.;
+      cout << " BeamTP::AutoNormalizeBeam()/Warning sum="<<sum<<" -> NormFac=1"<<endl;      
+    }
+  }
+  //--- return the Beam value for a given direction 
+  inline double Value(LongitudeLatitude const& ll)
+  {
+    return Value( UnitVector(ll) ); 
+  }
+  inline double Value(double theta, double phi)
+  { 
+    return Value( UnitVector(LongitudeLatitude(theta, phi)) ); 
+  }
+  //--- Beam value for a given direction 
+  double Value(UnitVector const& uvo)
+  {
+ // circular beam response 
+    double alp=acos(bdir_.Psc(uvo))*DoL_;
+    if (fggauss_) {
+      return  ( exp(-2.*alp*alp)*NormFac_ );   // Jiao : Check the factor 2 in the exp( )
+    }
+    else {  
+      if (alp<1.e-19)  return NormFac_;
+      alp*=M_PI;  double xa=2.*j1(alp)/alp;
+      return (NormFac_*xa*xa); 
+    }
+  }
+
+  // Return the beam in the form of HEALPix spherical map 
+  inline SphereHEALPix<double> getBeamSph(int m=128)
+  {
+    double sum=0.;
+    return getBeamSph(sum, m);
+  }
+  // Return the beam in the form of HEALPix spherical map - as well as the sum over all map pixels
+  SphereHEALPix<double> getBeamSph(double & sum, int m=128)
+  {
+    SphereHEALPix<double> sph(m);
+    sum=0.;
+    double bval;
+    for(int i=0; i<sph.NbPixels(); i++) {
+      double theta,phi;
+      sph.PixThetaPhi(i, theta, phi);
+      bval=Value(theta,phi);   
+      sum += bval;
+      sph[i]=bval;
+    }
+    return sph;
+  }
+
+protected:
+  bool fggauss_;
+  double D_, lambda_;
+  double DoL_;
+  UnitVector bdir_;    //  beam-direction   
+  double  bdir_theta_, bdir_phi_;   // auxiliary information for beam direction - optimized access to angles Theta() and Phi()
+  double NormFac_;  // Normalization factor 
+};
+
+
+#endif