Skip to content

GitLab

Commit 24d20cf1 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Ajout classe beam cross-cor, modifs programme trkacfit pour extraire les... 

Ajout classe beam cross-cor, modifs programme trkacfit pour extraire les cross-cor en vue du fit , Reza 20/12/2018
parent e37109b7
Hide whitespace changes
Inline Side-by-side
Showing with 106 additions and 12 deletions
acbeam.h
View file @ 24d20cf1
......@@ -79,6 +79,17 @@ public:
{
return NormFac_;
}
// return the dish size
inline double getDdish() const
{
return D_;
}
// return the associated beam wavelength
inline double getLambda() const
{
return lambda_;
}
// 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)
{
......
cxbeam.h 0 → 100644
View file @ 24d20cf1
/* PAON4 analysis software
Class representing an cross correlation beam (adapted from JSkyMap code)
R. Ansari, December 2018 */
#ifndef CXBEAM_H_SEEN
#define CXBEAM_H_SEEN
#include "acbeam.h"
#include "unitvector.h"
using namespace std;
using namespace SOPHYA;
class CxBeam { // Auto-correlation beam
public:
CxBeam(ACBeam & a1, ACBeam & a2, Vector3d & baseline)
: a1_(a1), a2_(a2), baseline_(baseline)
{
if (fabs(a1.getLambda()-a2.getLambda())>1e-9)
throw ParmError("CxBeam::CxBeam() - not same lambda in the two different beams !");
facphase_=2.*M_PI/a1.getLambda();
}
CxBeam(CxBeam const & c)
: a1_(c.a1_), a2_(c.a2_), baseline_(c.baseline_), facphase_(c.facphase_)
{
}
// copy operator
CxBeam& operator = (CxBeam const& a)
{
a1_=c.a1_; a2_=c.a2_; baseline_=c.baseline_; facphase_=c.facphase_;
return (*this);
}
//--- return the Beam value for a given direction
inline complex<double> Value(LongitudeLatitude const& ll)
{
return Value( UnitVector(ll) );
}
inline complex<double> Value(double theta, double phi)
{
return Value( UnitVector(LongitudeLatitude(theta, phi)) );
}
//--- Beam value for a given direction
complex<double> Value(UnitVector const& uvo)
{
double amp = sqrt(a1_.Value(uvo)*a2_.Value(uvo));
double phase=facphase_*(baseline_.Psc(uvo));
return (amp*complex<r_8>(cos(phase), sin(phase)) );
}
ACBeam a1_, a2_;
Vector3d baseline_;
double facphase_;
};
trkacfit.cc
View file @ 24d20cf1
......@@ -60,7 +60,9 @@ static string outcheckfilename; // output filename for fit parameters
static int prtlevel=1; // print level
//---- input data definition
static size_t NB_ANTENNES = 4; // Number of antenna ( 4 for PAON4
static size_t NB_ANTENNES = 4; // Number of antenna ( 4 for PAON4 )
static size_t NB_CXCORS = 6; // Number of Cx correlations ( 6 for PAON4 )
static size_t NTRK = 1; // Number of datasets = number of track files
//----- Visibility data sets
static string input_base_path; // base path (directory) for track and data files
......@@ -74,9 +76,11 @@ static double zenang=0.; // zenith angle = shift with declination, u
//-------- Declaration of functions in this file - code after the main, below
int decode_args(int narg, char* arg[]);
int read_p4_autocors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data);
int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data,
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr);
int read_srctracks(vector< vector<double> > & v_time_sat, vector< vector<double> > & v_sat_elev,
vector< vector<double> > & v_sat_azim, vector< SLinInterp1D > & v_interp_elev,
vector< SLinInterp1D > & v_interp_sazim);
......@@ -108,8 +112,10 @@ int main (int narg, char* arg[])
vector< vector< vector<double> > > vv_err(NTRK);
vector< vector<double> > v_min_data(NTRK);
vector< vector<double> > v_max_data(NTRK);
vector< vector< vector< complex<double> > > > vv_cxdata(NTRK);
vector< vector< vector<double> > > vv_cxerr(NTRK);
int tot_npoints = read_p4_autocors(v_time_data, vv_data, vv_err, v_min_data, v_max_data);
int tot_npoints = read_p4_autocxcors(v_time_data, vv_data, vv_err, v_min_data, v_max_data, vv_cxdata, vv_cxerr);
vector< vector<double> > v_time_sat(NTRK);
vector< vector<double> > v_sat_elev(NTRK);
......@@ -176,13 +182,16 @@ int main (int narg, char* arg[])
//-----------------------------------------------------------------------------------------
/* --Fonction-- */
int read_p4_autocors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data)
int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data,
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr)
{
int tot_npoints = 0; // total number of points for fit
cout << "---- trkacfit/read_p4_autocors() ; reading 4 PAON4 auto-correlation signals/DataTables ..."<<endl;
cout << "---- trkacfit/read_p4_autocxcors() ; reading 4 PAON4 auto-correlation & 6 Cross-cor signals/DataTables ..."<<endl;
const char * acname[4]={"V11","V22","V33","V44"};
const char * cxname[6]={"V12","V13","V14","V23","V24","V34"};
for(size_t j=0; j<dataflnm.size(); j++) {
string flnm = input_base_path+dataflnm[j]+".ppf";
cout << "1."<<j+1<<" Extracting data from data file DataTable: " << flnm<<endl
......@@ -194,8 +203,8 @@ int read_p4_autocors(vector< vector<double> > & v_time_data, vector< vector< vec
size_t ktime = dt_data.IndexNom("timesec");
vector<double> vtm;
dt_data.GetColumn(ktime, vtm);
vector< vector<double> > v_vac(4);
for(size_t ii=0; ii<4; ii++) { // 4 auto-correlations
vector< vector<double> > v_vac(NB_ANTENNES);
for(size_t ii=0; ii<NB_ANTENNES; ii++) { // 4 auto-correlations
size_t kac = dt_data.IndexNom(acname[ii]);
dt_data.GetColumn(kac, v_vac[ii]);
vector<double> vtmp, vetmp;
......@@ -204,19 +213,38 @@ int read_p4_autocors(vector< vector<double> > & v_time_data, vector< vector< vec
v_min_data[j].push_back(9.e19);
v_max_data[j].push_back(-9.e19);
}
vector< vector <complex<double> > > v_vcx(NB_CXCORS);
for(size_t ii=0; ii<NB_CXCORS; ii++) { // 6 cross-correlations
size_t kac = dt_data.IndexNom(cxname[ii]);
dt_data.GetColumn(kac, v_vcx[ii]);
vector< complex<double> > vtmp;
vector<double> vetmp;
vv_cxdata[j].push_back(vtmp);
vv_cxerr[j].push_back(vetmp);
}
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
vector< vector< complex<double> > > & v_cxdata = vv_cxdata[j];
vector< vector<double> > & v_cxerr = vv_cxerr[j];
for(size_t k=0; k<vtm.size(); k++) {
if ((vtm[k]<tstart[j])||(vtm[k]>tend[j])) continue;
v_time_data[j].push_back(vtm[k]);
for(size_t ii=0; ii<4; ii++) {
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
vector<double> & vac = v_vac[ii];
v_data[ii].push_back(vac[k]);
v_err[ii].push_back(0.1*sqrt(fabs(vac[k]))); // calcul d'erreur, a affiner
if (vac[k]<v_min_data[j][ii]) v_min_data[j][ii]=vac[k];
if (vac[k]>v_max_data[j][ii]) v_max_data[j][ii]=vac[k];
}
for(size_t ii=0; ii<NB_CXCORS; ii++) { // 6 cross-correlations
vector< complex<double> > & vcx = v_vcx[ii];
v_cxdata[ii].push_back(vcx[k]);
v_cxerr[ii].push_back(0.1*std::abs(vcx[k]));
}
}
tot_npoints += v_time_data[j].size(); // total number of points for fit
cout << " ... Done for " << j+1 << " data size="<<v_time_data[j].size()<<endl;
cout << " Data-AutoCor Min[A1...A4]=";
......
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