Commit edfcf9e8 authored by Reza Ansari's avatar Reza Ansari
Browse files

Implementation presque complete de la classe P4PhaseCor (non testee) et maj de...

Implementation presque complete de la classe P4PhaseCor (non testee) et maj de makefile, Reza 13/01/2018
parent 3c1d7776
......@@ -6,8 +6,8 @@ OBJ = ./Objs/
EXE = ./Objs/
# List of include files of this package, and .o files to handle dependencies
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h p4phcor.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq p4conv2fits msvis2dt visiavg filt_blind
......@@ -41,6 +41,9 @@ $(OBJ)/p4gnugain.o : p4gnugain.cc $(MYINCLISTHERE)
######
$(OBJ)/p4gvcor.o : p4gvcor.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4gvcor.o p4gvcor.cc
######
$(OBJ)/p4phcor.o : p4phcor.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4phcor.o p4phcor.cc
###############################################################
###### Compilation et link des executables
......
......@@ -38,23 +38,122 @@ See AnaFringes/postFringeV2.cc
r_8 val = 2.*atan(1./tan(M_PI*(xshift+0.5)));
*/
P4PhaseCor::P4PhaseCor(string const & filename,string const & filt_blind_name)
: parmsfilename_(filename),filt_blind_filename_(filt_blind_name)
P4PhaseCor::P4PhaseCor()
: fgcorphasok_(false)
{
//--- on initialise les 8 vecteurs des phases - phases nulle pour tout les feeds
for(size_t i=0; i<P4AVisiNumEncoder::getTotalNbFeeds(); i++) {
TVector< double > vph((sa_size_t)P4FreqBand::getP4NbFreqChannels);
v_phases_.push_back(vph);
}
}
P4PhaseCor::P4PhaseCor(string const & filename)
: fgcorphasok_(false)
{
//--- on initialise les 8 vecteurs des phases - phases nulle pour tout les feeds
for(size_t i=0; i<P4AVisiNumEncoder::getTotalNbFeeds(); i++) {
TVector< double > vph((sa_size_t)P4FreqBand::getP4NbFreqChannels);
v_phases_.push_back(vph);
}
readFitParamFile(filename);
}
double P4PhaseCor::getPhase(size_t numch, double freq )
void P4PhaseCor::readFitParamFile(string const & filename)
{
cout << "--- P4PhaseCor::readFitParamFile( filename= "<<filename<<" )"<<endl;
std::ifstream is(filename, std::ifstream::in);
string line;
size_t iv, nl;
char mot[32], ftype[32], key[4], feed[3];
double p[3];
while (!is.eof()) {
line = "";
getline(is, line);
if ((is.good() || is.eof()) && (line.length()>0) && (line[0]!='#')) {
p[0]=p[1]=p[2]= -9.e-99;
sscanf(line.c_str(),"%s %s %lg %lg %lg", mot, ftype, p, p+1, p+2);
key[0]=mot[0]; key[1]=mot[1]; key[2]=mot[2]; key[3]='\0';
feed[0]=mot[3]; feed[1]=mot[4]; feed[2]='\0';
if (strcmp(key,"Phi")!=0) continue;
sa_size_t fid = P4AVisiNumEncoder::getFeedNum(feed);
string sftyp = ftype;
if (fid < 0) {
cout<<"P4PhaseCor::readFitParamFile()/ERROR bad feed name= "<<feed<<" -> fid="<<fid<<endl;
continue;
}
if ((sftyp!="Sawtooth")&&(sftyp!="Linear")) {
cout<<"P4PhaseCor::readFitParamFile()/ERROR bad fitted function type= "<<sftyp<<endl;
continue;
}
TVector<double> phases((sa_size_t)P4FreqBand::getP4NbFreqChannels);
for (sa_size_t k=0; k<phases.Size(); k++) {
double freq = P4FreqBand::getP4Frequency((int)k);
if (sftyp=="Linear") phases(k) = jec_phase_linear(freq, p);
else phases(k) = jec_phase_sawtooth(freq, p);
}
cout << " P4PhaseCor::readFitParamFile()/Info: Setting phases (fit-type="<<sftyp<<") for feed="<<feed<<" -> fid="<<fid<<endl;
setFeedPhases(fid, phases);
}
}
}
void P4PhaseCor::setFeedPhases(size_t numfeed, TVector<double> vphase)
{
if (numfeed >= v_phases_.size())
throw ParmError("P4PhaseCor::setFeedPhases()/ERROR Out-of-range numfeed");
bool smo;
if (! vphase.CompareSizes(v_phases_[numfeed], smo) )
throw SzMismatchError("P4PhaseCor::setFeedPhases()/ERROR Bad size for vphase");
v_phases_[numfeed] = vphase;
fgcorphasok_=false;
}
TVector<double> P4PhaseCor::getFeedPhases(size_t numfeed) const
{
if (numfeed >= v_phases_.size())
throw ParmError("P4PhaseCor::getFeedPhases()/ERROR Out-of-range numfeed");
TVector<double> rvec;
rvec = v_phases_[numfeed];
return rvec;
}
double P4PhaseCor::getPhase(size_t numfeed, double freq) const
{
if (numfeed >= v_phases_.size())
throw ParmError("P4PhaseCor::getPhase()/ERROR Out-of-range numfeed");
sa_size_t kf = P4FreqBand::getP4FrequencyChannel(freq);
if ((kf < 0)||(kf>=v_phases_[numfeed].Size()))
throw ParmError("P4PhaseCor::getPhase()/ERROR Bad frequency value");
return v_phases_[numfeed](kf);
}
void P4PhaseCor::computePhaseCorrectionMatrix()
{
P4AVisiNumEncoder venc;
corphasemtx_.ReSize(venc.getTotalNbVisibilities(), P4FreqBand::getP4NbFreqChannels());
for(sa_size_t r=0; r<corphasemtx_.NRows(); r++) {
std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(r);
size_t i=(size_t)pij.first;
size_t j=(size_t)pij.second;
for(sa_size_t c=0; c<corphasemtx_.NCols(); c++) {
if (i == j ) corphasemtx_(r,c) = complex< r_4 >( (r_4)1., (r_4)0.);
else {
double deltaphi = v_phases_[i](c)-v_phases_[j](c); // CHECK : faut-il un signe moins ici ?
corphasemtx_(r,c) = complex< r_4 >( cos(deltaphi), sin(deltaphi) );
}
}
}
fgcorphasok_=true;
}
void P4gvCor::applyPhase(TMatrix <complex <float> > & mtx)
void P4PhaseCor::applyPhase(TMatrix <complex <float> > & mtx)
{
if (!fgcorphasok_) computePhaseCorrectionMatrix();
mtx.Mul(corphasemtx_);
return;
}
......
......@@ -26,16 +26,55 @@ using namespace std;
//! \brief class for managing phase corrections for PAON4 visbility matrices
class P4PhaseCor {
public:
//! Default constructor, phases for all feeds initialised to zero
P4PhaseCor();
//! \brief Constructor with the specification of the name of text file containing phase(freq) fit parameters (see readFitParamFile())
P4PhaseCor(string const & filename);
/*! read the text file \b filename containing the fit paramaters for phase as a function of frequency
file format:
\verbatim
Phi2H Linear 1375.152587890625 -2.5804805755615234375 -0.0041318359374999998751
Phi3H Sawtooth 0.05565275701111126877 -1435.2843925781251073
Phi4H Sawtooth 0.070963965382242408242 -1393.1309238281251055
\endverbatim
*/
void readFitParamFile(string const & filename);
//! set/changes the values of the phase=func(freq) (vector \b vphase) for the feed \b numfeed (0...7 -> 1H,2H...,3V,4V)
void setFeedPhases(size_t numfeed, TVector<double> vphase);
//! return the vector of values of the phase=func(freq) (vector \b vphase) for the feed \b numfeed (0...7 -> 1H,2H...,3V,4V)
TVector<double> getFeedPhases(size_t numfeed) const;
//! return the phase correction (in radian) to be applied to feed indentified by numch and frequency freq in MHz
double getPhase(size_t numch, double freq );
double getPhase(size_t numfeed, double freq ) const;
//! Applies the phase corrections
void applyPhase(TMatrix <complex <float> > & mtx);
//! computes the phase correction matrix (this method is automatically called when needed)
void computePhaseCorrectionMatrix();
//---- See AnaFringes/postFringeV2.cc
/*! \brief return the phase value as a linear function of frequency, using Jean-Eric linear fit
parameters p[0,1,2] : p[1]+(freq-p[0])*p[2] */
static inline double jec_phase_linear(double freq, double * p)
{
return p[1]+p[2]*(freq-p[0]);
}
/*! \brief return the phase value as a sawtooth function of frequency, using Jean-Eric Sawtooth fit
parameters p[0,1] : 2 ArgTan[Cot(Pi{p[0](x+p[1])+0.5})] */
static inline double jec_phase_sawtooth(double freq, double * p)
{
return 2.*atan(1./tan(M_PI*((p[0]*(freq+p[1]))+0.5)));
}
protected:
// for each feed, there is a vector which contains the value of the phase as a function of frequency
std::vector< TVector< double > > v_phases;
// This matrix contains for each feed and frequancy exp (i phase)
TMatrix <complex <float> > cmplxphases_;
std::vector< TVector< double > > v_phases_;
bool fgcorphasok_; // true -> corphasemtx_ OK
// This matrix contains for each feed pair and frequancy exp (i phase)
TMatrix <complex <float> > corphasemtx_;
};
......
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