Commit f0c53ba6 authored by perdereau's avatar perdereau
Browse files

few cosmetic changes ; add handling code for new polynomial phase

correction  ; correct bugs in vis2ra
parent 6fae0b71
......@@ -640,9 +640,11 @@ P4RAMapUtil::P4RAMapUtil(double RAcenter, double RAwidth, double RAresol_mn)
sa_size_t P4RAMapUtil::getMapIndex(double ra)
{
double delra = ra-RAcenter_;
sa_size_t idx = (sa_size_t)(delra*invres_hours_)+size_/2;
//DBG cout << " delra="<<delra<<" idx="<<idx<<endl;
//DBG
//cout << " delra="<<delra<<" idx="<<idx<<" center "<< RAcenter_ <<" size "<< size_<< endl;
if (idx<0) { // on ajoute 24 heures a ra et on recommece
idx += size_24hours_;
}
......
......@@ -250,6 +250,7 @@ public:
P4RAMapUtil(P4AnaParams const& param);
P4RAMapUtil(double RAcenter, double RAwidth=24., double RAresol_mn=1.);
inline sa_size_t getMapSize() { return size_; }
//! return the index in map for ra (in decimal hours) - return negative if outside map
sa_size_t getMapIndex(double ra);
......
......@@ -72,7 +72,7 @@ void P4gvCor::InitNormalisation()
if there are more that 3 parameters on each line (row) , the 4th parameter is considered as the overall normalisation parameter
\param filt_blind_name : filetered blind channel signal PPF filename, signal as a function of time
*/
void P4gvCor::setGtCorrection(string const & param_filename, string const & filt_blind_filename)
void P4gvCor::setGtCorrection(string const & param_filename, string const & filt_blind_name)
{
TArray <double> arr ;
......@@ -99,8 +99,8 @@ void P4gvCor::setGtCorrection(string const & param_filename, string const & filt
throw IOExc("P4gvCor::setGtCorrection() : NO DATA FILE ");
}
cout << "P4gvCor::setGtCorrection() - reading PPF fileterd blind channel signal from "<<filt_blind_filename<<endl;
PInPersist pinp(filt_blind_filename);
cout << "P4gvCor::setGtCorrection() - reading PPF fileterd blind channel signal from "<<filt_blind_name<<endl;
PInPersist pinp(filt_blind_name);
pinp >> PPFNameTag("TVFilt_Times") >> filt_vtimes;
pinp >> PPFNameTag("TVFilt_Values") >> filt_values;
......@@ -110,7 +110,7 @@ void P4gvCor::setGtCorrection(string const & param_filename, string const & filt
for (sa_size_t k =0 ; k<filt_vtimes.Size() ; k++){
x.push_back(filt_vtimes(k));
y.push_back(filt_values(k));
cout << " decoding template : t="<<filt_vtimes(k) << " val = "<< filt_values(k) << endl;
}
my_interp_.DefinePoints (x,y);
Gt_corr_param_defined_ = true;
......@@ -174,8 +174,8 @@ void P4gvCor::applyGain(TMatrix <complex <float> > & mtx, TimeStamp const & tms
double ggj = slopes_[j]*vbl+offsets_[j];
gg=sqrt(ggi*ggj);
}
if (fgdonorm) gg *= visi_norm_[k];
visi_gains[k]=1./gg;
if (fgdonorm) visi_gains[k]*= visi_norm_[k];
}
......
......@@ -68,6 +68,18 @@ double P4PhaseCor::jec_phase_poly(double freq, double * p){
}
return(ret);
}
double P4PhaseCor::phase_pypoly(double freq, double * p){
double ret = 0.;
double del = 1.;
int deg = p[0];
double delta = freq-p[1];
for (int k=0 ; k <= deg ; k ++){
ret += p[k+2] *pow(delta,deg-k);
}
return(ret);
}
void P4PhaseCor::readFitParamFile(string const & filename)
{
......@@ -85,14 +97,17 @@ void P4PhaseCor::readFitParamFile(string const & filename)
sscanf(line.c_str(),"%s %s %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", mot, ftype, p, p+1, p+2,p+3,p+4,p+5,p+6,p+7,p+8,p+9,p+10,p+11,p+12);
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;
if (strcmp(key,"Phi")!=0) {
cout << " Phi not in "<<key<<endl;
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")&&(sftyp!="Polynomial")) {
if ((sftyp!="Sawtooth")&&(sftyp!="Linear")&&(sftyp!="Polynomial")&&(sftyp!="Polynome")) {
cout<<"P4PhaseCor::readFitParamFile()/ERROR bad fitted function type= "<<sftyp<<endl;
continue;
}
......@@ -101,7 +116,10 @@ void P4PhaseCor::readFitParamFile(string const & filename)
double freq = P4FreqBand::getP4Frequency((int)k);
if (sftyp=="Linear") phases(k) = jec_phase_linear(freq, p);
else if (sftyp=="Sawtooth") phases(k) = jec_phase_sawtooth(freq, p);
else phases(k) = jec_phase_poly(freq, p);
else if (sftyp=="Polynomial") phases(k) = jec_phase_poly(freq, p);
else phases(k) = phase_pypoly(freq, p);
}
cout << " P4PhaseCor::readFitParamFile()/Info: Setting phases (fit-type="<<sftyp<<") for feed="<<feed<<" -> fid="<<fid<<endl;
setFeedPhases(fid, phases);
......@@ -142,6 +160,7 @@ void P4PhaseCor::setFeedPhases(size_t numfeed, TVector<double> vphase)
if (! vphase.CompareSizes(v_phases_[numfeed], smo) )
throw SzMismatchError("P4PhaseCor::setFeedPhases()/ERROR Bad size for vphase");
v_phases_[numfeed] = vphase;
fgcorphasok_=false;
}
......@@ -177,6 +196,7 @@ void P4PhaseCor::computePhaseCorrectionMatrix()
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) );
}
}
}
......
......@@ -39,6 +39,13 @@ public:
Phi3H Sawtooth 0.05565275701111126877 -1435.2843925781251073
Phi4H Sawtooth 0.070963965382242408242 -1393.1309238281251055
\endverbatim
or
\verbatim
2.0 Polynom 1373.0877047317845 6.0 6.285802318653258e-13 4.9695307019287343e-11 -1.0776922683612119e-08 -1.0260776080700936e-06 3.6209213468801214e-05 0.0031258250785635264 -2.474941671890262
3.0 Polynom 1373.3238707042551 7.0 -6.906887825590439e-15 8.688203385849563e-13 1.6463244340003824e-10 -1.6684660244478277e-08 -1.1777761116566551e-06 5.40549458809292e-05 -0.005545857771127432 -0.19635476731191293
4.0 Polynom 1376.207728968306 10.0 -1.930997090869567e-20 -1.2452810626934814e-18 6.380842950410297e-16 2.336938329098397e-14 -7.48472627594635e-12 3.316167366128346e-11 3.403078006610292e-08 -2.1592536781344065e-06 -2.787174415656901e-05 0.004667981257677257 0.22993894526629288
\endverbatim
*/
void readFitParamFile(string const & filename);
......@@ -80,6 +87,11 @@ public:
p[2]+p[3]*(freq-p[0])+... */
static double jec_phase_poly(double freq, double * p);
/*! \brief return the phase value as a polynom function of frequency, using Jean-Eric Sawtooth fit
parameters p[0,1,...] : p[0] center freq p[1]= degree of poly, p[2].... p[n] coeffs
in decreasing power order (python-like convention)*/
static double phase_pypoly(double freq, double * p);
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_;
......
......@@ -124,11 +124,11 @@ int main(int narg, const char* arg[])
// y avait cnt++ & continue ?
}// end of resizing
else{
// Accumulating
visum += vismtx;
meanvismtx += vismtx;
}
if (I==0) dateobs=cfdate; // start filling a new time bin
// Accumulating
I++; // incrementing DeltaTime counter
......
......@@ -83,6 +83,7 @@ int main(int narg, const char* arg[])
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<" DeltaIAvg="<<deltaIavg<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
P4RAMapUtil ram(params);
......@@ -122,7 +123,13 @@ int main(int narg, const char* arg[])
// vecteur de noms
vector <string> ext_names;
// un vecteur avec les temps
// TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg);
TVector< double > RAVec(ram.getMapSize());
double rast=ram.RAcenter_-ram.RAwidth_/2.+.5*ram.RAresol_mn_/60.;
for (int kk=0 ; kk< ram.getMapSize() ; kk++){
RAVec[kk]=rast;
rast += ram.RAresol_mn_/60.;
}
TMatrix< complex<r_4> > vismtx;
TimeStamp dateobs, cfdate,datestart;
......@@ -169,7 +176,7 @@ int main(int narg, const char* arg[])
ramSY=vismtx.NCols()/TFMfbin;
racntmap.SetSize2D(ramSX, ramSY, true); // true pour forcer la mise a zero
//allocating 8 Auto-Corr RA-frequency maps
cout<<"vis2ra/Info: allocating 8 AutoCor rightascension-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
cout<<"vis2ra/Info: allocating 8 AutoCor RA-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
for(int k=0; k<8; k++) vramac.push_back( TArray< r_4 >(ramSX, ramSY) );
cout <<" and 8 for the the variance maps "<<endl;
for(int k=0; k<8; k++) vramac_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
......@@ -180,7 +187,7 @@ int main(int narg, const char* arg[])
cout << "and the 2x6 for squares of rela & imag parts " << endl;
for(int k=0; k<6; k++) vram_rp_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
for(int k=0; k<6; k++) vram_ip_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
cout << " HH xcor done "<<endl;
if (params.doVV_){
//allocating 6 Cross-Corr H-H time-frequency maps
cout<<"vis2ra/Info: allocating V-V cross-cor RA-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
......@@ -188,32 +195,29 @@ int main(int narg, const char* arg[])
cout << "and the 2x6 for squares of rela & imag parts " << endl;
for(int k=0; k<6; k++) vram_vv_rp_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
for(int k=0; k<6; k++) vram_vv_ip_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
} // endif doVV
}
if (params.doHV_){
//allocating 16 Cross-Corr H-H time-frequency maps
cout<<"vis2ra/Info: allocating H-V cross-cor RA-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
for(int k=0; k<16; k++) vram_hv.push_back( TArray< complex<r_4> >(ramSX, ramSY) );
cout << "and the 2x6 for squares of rela & imag parts " << endl;
for(int k=0; k<16; k++) vram_hv_rp_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
for(int k=0; k<16; k++) vram_hv_ip_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
}
for(int k=0; k<16; k++) vram_hv_ip_sq.push_back( TArray< r_4 >(ramSX, ramSY) );
} // endif doHV
// recupere le jour de depart @ 0h
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
cout << "vis2ra/Info: alloc done cnt=" <<cnt<<endl;
}
cout << "vis2ra/Info: alloc done cnt=" <<cnt<<endl;
} // end if cnt==00
dateobs=cfdate;
double racourant = P4Coords::RAFromTimeTU(dateobs);
sa_size_t idxra = ram.getMapIndex(racourant);
//cout << idxra <<" return by ram.getMapIndex for "<< racourant << endl;
cnt++;
if (idxra < 0) continue;
// On est dans la carte en RA , on remplit d'abord la carte de count
for(sa_size_t jy=0; jy<ramSY; jy++) racntmap(idxra,jy) += 1.;
int nbval=0;
//---- On s'occupe d'abord des autocorrelations P1 ... P8
for(size_t k=0; k<KVAC.size(); k++) { // Loop over the 8 auto-correlations
TVector<r_4> vac = real(vismtx.Row(KVAC[k]));
......@@ -221,10 +225,13 @@ int main(int narg, const char* arg[])
TArray< r_4 > & ramap = vramac[k];
TArray< r_4 > & ramap_sq = vramac_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) { // frequency binning
nbval = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Size();
ramap(idxra, jy) += vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramap_sq(idxra, jy) += vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 8 AutoCor
// On est dans la carte en RA , on remplit aussi la carte de count
for(sa_size_t jy=0; jy<ramSY; jy++) racntmap(idxra,jy) += nbval;
//---- On s'occupe des 6 cross-correlations 1H-2H ... 3H-4H
for(size_t k=0; k<KVCXHH.size(); k++) { // loop over the 6 Xcor
......@@ -235,9 +242,9 @@ int main(int narg, const char* arg[])
TArray< r_4 > & ramapsqpr = vram_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_ip_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) {
ramap(idxra, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramap(idxra, jy) += vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) += vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) += vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 6 Xcor
......@@ -251,9 +258,9 @@ int main(int narg, const char* arg[])
TArray< r_4 > & ramapsqpr = vram_vv_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_vv_ip_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) {
ramap(idxra, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramap(idxra, jy) += vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) += vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) += vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 6 VV Xcor
}// end doVV_
......@@ -268,9 +275,9 @@ int main(int narg, const char* arg[])
TArray< r_4 > & ramapsqpr = vram_hv_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_hv_ip_sq[k];
for(sa_size_t jy=0; jy<ramSY; jy++) {
ramap(idxra, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramap(idxra, jy) += vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpr(idxra, jy) += vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
ramapsqpi(idxra, jy) += vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
} //----- end of loop over the 16 HV Xcor
}// end doHV_
......@@ -280,7 +287,7 @@ int main(int narg, const char* arg[])
// timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart); // centre du bin
// ... done
cnt++;
if ( (cnt>0) && (cnt%100==0) ) {
cout << " vis2ra/Info: RA-map fill, VisMtxCount= "<<cnt<<endl;
// cout << " vis2ra/Info: RA-map fill, VisMtxCount="<<cnt<< " /Max="<<Imax<<" DateObs="<<dateobs<<endl;
......@@ -514,7 +521,6 @@ int main(int narg, const char* arg[])
ext_names.push_back(viramHV_names[k]);
(*fos)<< ramap_sqpi;
}
}
}// end option HV
......@@ -527,9 +533,12 @@ int main(int narg, const char* arg[])
avg_freqs(kf) = frbase ;
poram << PPFNameTag("FreqVec") << avg_freqs;
poram << PPFNameTag("RAVec") << RAVec;
if (fos != NULL) {
ext_names.push_back("Frequences");
(*fos)<< avg_freqs ;
ext_names.push_back("RAs");
(*fos)<<RAVec ;
cout << " number of objs in fits "<< ext_names.size() << endl;
cout << ext_names << endl;
delete(fos);
......
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