Commit 3ea34d2f authored by perdereau's avatar perdereau
Browse files

change type for min/max freq (was int...) + upgrade vis2ra to handle

vv,hv`
parent fe54e5d6
......@@ -78,7 +78,7 @@ P4gnuGain::P4gnuGain(string const & gainfilename, double gthr)
P4FreqBand myp4fre;
int maxg=gains_.NCols()/2;
double maxg=gains_.NCols()/2.;
for (sa_size_t c = gains_.NCols()/2 ; c < gains_.NCols() ; c++) {
int ok=0;
for(sa_size_t r=0;r<gains_.NRows();r++) {
......@@ -87,23 +87,25 @@ P4gnuGain::P4gnuGain(string const & gainfilename, double gthr)
}
if (ok==1) break;
maxg=c;
maxg= (double) c;
}
int ming=gains_.NCols()/2;
double ming=gains_.NCols()/2;
for (sa_size_t c = gains_.NCols()/2 ; c >=0 ; c--) {
int ok=0;
for(sa_size_t r=0;r<gains_.NRows();r++) {
if (gains_(r,c) ==0 ) { ok=1; }
if (gains_(r,c) == 0 ) { ok=1; }
}
if (ok==1) break;
ming=c;
ming=(double) c;
}
minGoodF_ = myp4fre.freqstart_ + ming* myp4fre.getP4FreqResolution();
maxGoodF_ = myp4fre.freqstart_ + maxg* myp4fre.getP4FreqResolution();
cout << "Frequences limites "<< minGoodF_ << " <> " << maxGoodF_ << endl;
for(sa_size_t r=0; r<vmtxgains_.NRows(); r++) {
std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(r);
sa_size_t i=pij.first;
......
......@@ -43,15 +43,15 @@ public:
//! return the multiplicative complex gain matrix, applicable to the full visibility matrix.
inline SOPHYA::TMatrix< complex<r_4> > const & getVisiGainMatrix() const { return vmtxgains_; }
inline int minGoodF(){return minGoodF_;}
inline int maxGoodF(){return maxGoodF_;}
inline double minGoodF(){return minGoodF_;}
inline double maxGoodF(){return maxGoodF_;}
protected:
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
static SOPHYA::Matrix computeGain2(SOPHYA::Matrix & acs, int SZW=96, float frackeep=0.75);
// ----- fonction pour normaliser les gains
static SOPHYA::Vector normalizeGain2(SOPHYA::Matrix & gain);
int minGoodF_,maxGoodF_;
double minGoodF_,maxGoodF_;
std::string gainfilename_;
SOPHYA::TMatrix<r_8> gains_;
SOPHYA::TVector<r_8> gnorm_;
......
......@@ -89,6 +89,8 @@ int main(int narg, const char* arg[])
P4AVisiNumEncoder visiencod;
vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
vector<sa_size_t> KVCXVV = visiencod.getAllVCrossCor();
vector<sa_size_t> KVCXHV = visiencod.getAllHVCrossCor();
cout << " List of AutoCorrelation rows:"<<endl;
for(size_t k=0; k<KVAC.size(); k++) {
cout << "KVAC["<<k<<"]="<<KVAC[k]<<" ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
......@@ -129,12 +131,25 @@ int main(int narg, const char* arg[])
int cnt=0, cntnt=0, pcntnt=0;
//----- Count map for projection into (RA,Freq) maps
TArray< r_4 > racntmap;
//----- 6 H-H cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vram;
//----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vram_rp_sq;
vector< TArray< r_4 > > vram_ip_sq;
//----- 6 V-V cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vram_vv;
//----- 6 V-V cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vram_vv_rp_sq;
vector< TArray< r_4 > > vram_vv_ip_sq;
//----- 16 H-V cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vram_hv;
//----- 16 H-V cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vram_hv_rp_sq;
vector< TArray< r_4 > > vram_hv_ip_sq;
//----- 8 auto-corr TimeFrequency maps
vector< TArray< r_4 > > vramac;
//----- 8 auto-corr TimeFrequency variance maps
......@@ -165,6 +180,28 @@ 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) );
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;
for(int k=0; k<6; k++) vram_vv.push_back( TArray< complex<r_4> >(ramSX, ramSY) );
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) );
}
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) );
}
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
......@@ -179,7 +216,7 @@ int main(int narg, const char* arg[])
for(sa_size_t jy=0; jy<ramSY; jy++) racntmap(idxra,jy) += 1.;
//---- 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(k));
TVector<r_4> vac = real(vismtx.Row(KVAC[k]));
TVector<r_4> vacsq = vac.SquareElt();
TArray< r_4 > & ramap = vramac[k];
TArray< r_4 > & ramap_sq = vramac_sq[k];
......@@ -191,7 +228,7 @@ int main(int narg, const char* arg[])
//---- 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
TVector< complex<r_4> > vcx = vismtx.Row(k);
TVector< complex<r_4> > vcx = vismtx.Row(KVCXHH[k]);
TVector<r_4> vcxprsq = (real(vcx)).SquareElt();
TVector<r_4> vcxpisq = (imag(vcx)).SquareElt();
TArray< complex<r_4> > & ramap = vram[k];
......@@ -204,6 +241,41 @@ int main(int narg, const char* arg[])
}
} //----- end of loop over the 6 Xcor
if (params.doVV_){
//---- On s'occupe des 6 cross-correlations V-V si demande
for(size_t k=0; k<KVCXVV.size(); k++) { // loop over the 6 Xcor
TVector< complex<r_4> > vcx = vismtx.Row(KVCXVV[k]);
TVector<r_4> vcxprsq = (real(vcx)).SquareElt();
TVector<r_4> vcxpisq = (imag(vcx)).SquareElt();
TArray< complex<r_4> > & ramap = vram_vv[k];
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();
}
} //----- end of loop over the 6 VV Xcor
}// end doVV_
if (params.doHV_){
//---- On s'occupe des 16 cross-correlations H-V si demande
for(size_t k=0; k<KVCXHV.size(); k++) { // loop over the 16 Xcor
TVector< complex<r_4> > vcx = vismtx.Row(KVCXHV[k]);
TVector<r_4> vcxprsq = (real(vcx)).SquareElt();
TVector<r_4> vcxpisq = (imag(vcx)).SquareElt();
TArray< complex<r_4> > & ramap = vram_hv[k];
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();
}
} //----- end of loop over the 16 HV Xcor
}// end doHV_
// double tdif = cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
// timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart); // centre du bin
// ... done
......@@ -222,6 +294,7 @@ int main(int narg, const char* arg[])
cout<<"vis2ra/Info: Normalising RA-maps ... " << endl;
for(size_t k=0; k<KVAC.size(); k++) { // loop over 8 Autocor
TArray< r_4 > & ramap = vramac[k];
TArray< r_4 > & ramap_sq = vramac_sq[k];
......@@ -234,6 +307,7 @@ int main(int narg, const char* arg[])
}
}
} // end of loop over 8 Autocor
for(size_t k=0; k<KVCXHH.size(); k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & ramap = vram[k];
TArray< r_4 > & ramapsqpr = vram_rp_sq[k];
......@@ -250,6 +324,51 @@ int main(int narg, const char* arg[])
}
} // End of loop over the 6 Xcor
// normalisation of VV is asked
if (params.doVV_){
cout<<"vis2ra/Info: Normalising RA-maps (VV) ... " << endl;
for(size_t k=0; k<KVCXVV.size(); k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & ramap = vram_vv[k];
TArray< r_4 > & ramapsqpr = vram_vv_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_vv_ip_sq[k];
for(sa_size_t ix=0; ix<ramSX; ix++) {
for(sa_size_t jy=0; jy<ramSY; jy++) {
if (racntmap(ix,jy) < 1.e-9) continue;
r_4 fnorm = 1./racntmap(ix,jy);
complex<r_4> znorm = complex<r_4>(fnorm,0.);
ramap(ix, jy) *= znorm;
ramapsqpr(ix, jy) *= fnorm;
ramapsqpi(ix, jy) *= fnorm;
}
}
} // End of loop over the 6 Xcor
}//end of doVV_
// normalisation of HV is asked
if (params.doHV_){
cout<<"vis2ra/Info: Normalising RA-maps (HV) ... " << endl;
for(size_t k=0; k<KVCXHV.size(); k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & ramap = vram_hv[k];
TArray< r_4 > & ramapsqpr = vram_hv_rp_sq[k];
TArray< r_4 > & ramapsqpi = vram_hv_ip_sq[k];
for(sa_size_t ix=0; ix<ramSX; ix++) {
for(sa_size_t jy=0; jy<ramSY; jy++) {
if (racntmap(ix,jy) < 1.e-9) continue;
r_4 fnorm = 1./racntmap(ix,jy);
complex<r_4> znorm = complex<r_4>(fnorm,0.);
ramap(ix, jy) *= znorm;
ramapsqpr(ix, jy) *= fnorm;
ramapsqpi(ix, jy) *= fnorm;
}
}
} // End of loop over the Xcor
}// end of if(doHV_)
if (fitsoutfile.length()>=1){
cout << " fitsoutfile" <<fitsoutfile<< endl;
fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
......@@ -313,9 +432,92 @@ int main(int narg, const char* arg[])
ext_names.push_back(viramCC_names[k]);
(*fos)<< ramap_sqpi;
}
}
if (params.doVV_){ // Option VV
cout<<" vis2ra/Info: Saving 6 V-V cross-corr time-frequency maps to PPF file "<<outfile<<endl;
const char* ramVV_names[6]={"RAM_1V2V", "RAM_1V3V", "RAM_1V4V", "RAM_2V3V", "RAM_2V4V", "RAM_3V4V"};
const char* vrramVV_names[6]={"RVARRAM_1V2V", "RVARRAM_1V3V", "RVARRAM_1V4V", "RVARRAM_2V3V", "RVARRAM_2V4V", "RVARRAM_3V4V"};
const char* viramVV_names[6]={"IVARRAM_1V2V", "IVARRAM_1V3V", "IVARRAM_1V4V", "IVARRAM_2V3V", "IVARRAM_2V4V", "IVARRAM_3V4V"};
for(int k=0; k<6; k++) { // loop over the 6 VV xCor
TArray< complex<r_4> > & ramap = vram_vv[k];
TArray< r_4 > & ramap_sqpr = vram_vv_rp_sq[k];
TArray< r_4 > & ramap_sqpi = vram_vv_ip_sq[k];
poram << PPFNameTag(ramVV_names[k]) << ramap;
if (fos != NULL) {
ext_names.push_back(string(ramVV_names[k])+"_real");
(*fos)<< real(ramap);
ext_names.push_back(string(ramVV_names[k])+"_imag");
(*fos)<< imag(ramap);
}
TArray< r_4 > ramapr = real(ramap);
TArray< r_4 > ramapi = imag(ramap);
ramap_sqpr -= ramapr.SquareElt() ;
ramapi = ramapi.MulElt(ramapi,ramapi) ;
ramap_sqpi -= ramapi.SquareElt() ;
poram << PPFNameTag(vrramVV_names[k]) << ramap_sqpr;
poram << PPFNameTag(viramVV_names[k]) << ramap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrramVV_names[k]);
(*fos)<< ramap_sqpr;
ext_names.push_back(viramVV_names[k]);
(*fos)<< ramap_sqpi;
}
}// end loop on 6 VV xCor
}// end of option VV
if (params.doHV_){ // Option HV
const char* ramHV_names[16]={"RAM_1H1V", "RAM_1H2V", "RAM_1H3V", "RAM_1H4V",
"RAM_2H1V", "RAM_2H2V", "RAM_2H3V", "RAM_2H4V",
"RAM_3H1V", "RAM_3H2V", "RAM_3H3V", "RAM_3H4V",
"RAM_4H1V", "RAM_4H2V", "RAM_4H3V", "RAM_4H4V"
};
const char* vrramHV_names[16]={"RVARRAM_1H1V", "RVARRAM_1H2V", "RVARRAM_1H3V", "RVARRAM_1H4V",
"RVARRAM_2H1V", "RVARRAM_2H2V", "RVARRAM_2H3V", "RVARRAM_2H4V",
"RVARRAM_3H1V", "RVARRAM_3H2V", "RVARRAM_3H3V", "RVARRAM_3H4V",
"RVARRAM_4H1V", "RVARRAM_4H2V", "RVARRAM_4H3V", "RVARRAM_4H4V"
};
const char* viramHV_names[16]={"IVARRAM_1H1V", "IVARRAM_1H2V", "IVARRAM_1H3V", "IVARRAM_1H4V",
"IVARRAM_2H1V", "IVARRAM_2H2V", "IVARRAM_2H3V", "IVARRAM_2H4V",
"IVARRAM_3H1V", "IVARRAM_3H2V", "IVARRAM_3H3V", "IVARRAM_3H4V",
"IVARRAM_4H1V", "IVARRAM_4H2V", "IVARRAM_4H3V", "IVARRAM_4H4V"
};
for(int k=0; k<16; k++) { // loop over the 16 HV xCor
TArray< complex<r_4> > & ramap = vram_hv[k];
TArray< r_4 > & ramap_sqpr = vram_hv_rp_sq[k];
TArray< r_4 > & ramap_sqpi = vram_hv_ip_sq[k];
poram << PPFNameTag(ramHV_names[k]) << ramap;
if (fos != NULL) {
ext_names.push_back(string(ramHV_names[k])+"_real");
(*fos)<< real(ramap);
ext_names.push_back(string(ramHV_names[k])+"_imag");
(*fos)<< imag(ramap);
}
TArray< r_4 > ramapr = real(ramap);
TArray< r_4 > ramapi = imag(ramap);
ramap_sqpr -= ramapr.SquareElt() ;
ramapi = ramapi.MulElt(ramapi,ramapi) ;
ramap_sqpi -= ramapi.SquareElt() ;
poram << PPFNameTag(vrramHV_names[k]) << ramap_sqpr;
poram << PPFNameTag(viramHV_names[k]) << ramap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrramHV_names[k]);
(*fos)<< ramap_sqpr;
ext_names.push_back(viramHV_names[k]);
(*fos)<< ramap_sqpi;
}
}
}// end option HV
// --- FIN sauvegarde cartes temps-frequence
// resu.Update();
P4FreqBand myp4fre;
......
......@@ -638,6 +638,8 @@ int main(int narg, const char* arg[])
lim_freq(0) = myp4fre.freqstart_;
lim_freq(1) = myp4fre.freqend_;
}
cout << " frequences limites "<< lim_freq(0) <<" ; "<< lim_freq(1)<<endl;
potfm << PPFNameTag("FreqLims") << lim_freq ;
potfm << PPFNameTag("TimeVec") << timevec ;
potfm << PPFNameTag("RAVec") << ravec ;
......
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