Commit d5491102 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Ajout fichier acbeam.h pour fit angle de visee des antennes, Reza 12/12/2018

parent 4dce16a0
/* 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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment