Commit 7dcd115e authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Amelioration programme visamm.cc : projection en fonction de ra, argument...

Amelioration programme visamm.cc : projection en fonction de ra, argument -elev et autres ameliorations... Reza 18/07/2018
parent baa2cc43
......@@ -38,7 +38,7 @@ size_t P4FreqSelectorFilterMgr::readFreqBandDefinitioFile(string const & filenam
size_t P4FreqSelectorFilterMgr::addBand(P4FreqBand& fb)
{
v_bands_.push_back( pair<sa_size_t , sa_size_t>(fb.getFirstFreqChannel() , fb.getLastFreqChannel() ) );
weight_bands_.push_back(complex<float>((float)(1./(double)(fb.getLastFreqChannel()-fb.getFirstFreqChannel())), 1.));
weight_bands_.push_back(complex<float>((float)(1./(double)(fb.getLastFreqChannel()-fb.getFirstFreqChannel()+1)), 0.));
v_freqbands_.push_back(fb);
return v_bands_.size();
}
......
......@@ -48,7 +48,7 @@ int Usage(void)
{
cout << " --- visamm.cc : Reads PPF files produced by mfacq and produces visibility \n"
<< " arrays suitable as input for map making step \n"
<< " Usage: visamm [-arguments] \n"
<< " Usage: visamm [-arguments] [-elev elevationAngleDegree] \n"
<< " Note that Frequency bands should be specified \n" << endl;
P4AnaParams::UsageOptions();
cout<< endl;
......@@ -64,6 +64,14 @@ int main(int narg, const char* arg[])
P4AnaParams params;
params.DecodeArgs(narg, arg);
// -- if elevation angle set/defined
bool fgelevset=false;
double elevation=0.;
if ((params.lastargs_.size()>1)&&(params.lastargs_[0]=="-elev")) {
elevation=atof(params.lastargs_[1].c_str());
fgelevset=true;
}
string outfile = params.outfile_;
if (outfile.length()<1) outfile = "visamm.ppf";
......@@ -82,7 +90,7 @@ int main(int narg, const char* arg[])
cout <<"visamm/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
if (fgelevset) cout << " ... elevation set to "<<elevation<<endl;
// ---
int rc = 0;
......@@ -91,7 +99,7 @@ int main(int narg, const char* arg[])
HiStatsInitiator _inia;
ResourceUsage resu;
P4RAMapUtil ram(params);
// Frequency band definition
P4FreqSelectorFilterMgr freqbands;
cout << "visamm/Info : Computing Visibility array for " << params.fbands_.size() << " Frequency bands"<<endl;
for(size_t i=0; i<params.fbands_.size(); i++) {
......@@ -105,6 +113,9 @@ int main(int narg, const char* arg[])
vector<sa_size_t> KVCXVV = visiencod.getAllVCrossCor();
vector<sa_size_t> KVCXHV = visiencod.getAllHVCrossCor();
// ra-map definition
P4RAMapUtil ram(params);
VisiP4WindowReader wreader(params);
long Imin = wreader.getReader().getSerialFirst();
......@@ -125,12 +136,22 @@ int main(int narg, const char* arg[])
double mttag;
int cnt=0, cntnt=0, pcntnt=0;
//----- Array of visibilities for HH , one array for each frequency
vector< TArray< complex<r_4> > > vamm;
for(size_t i=0; i<freqbands.size(); i++)
vamm.push_back( TArray< complex<r_4> >((sa_size_t)wreader.getTotalNbWindows(), (sa_size_t)4+KVCXHH.size()) );
sa_size_t kt=0;
vector< TArray<r_4> > vammac;
for(size_t i=0; i<freqbands.size(); i++) {
// Arrays of 4 auto-correlations
vammac.push_back( TArray< r_4 >((sa_size_t)ram.getMapSize(), 4 ) );
// Arrays of 1 Auto-correlation + 6 cross-correlation
vamm.push_back( TArray< complex<r_4> >((sa_size_t)ram.getMapSize(), (sa_size_t)(KVCXHH.size()+1) ) );
}
// vamm.push_back( TArray< complex<r_4> >((sa_size_t)wreader.getTotalNbWindows(), (sa_size_t)1+KVCXHH.size()) );
//----- Fill count for ra-maps
TArray< r_4 > racntmap((sa_size_t)ram.getMapSize());
sa_size_t kt=0; //-- time index
while (fgok) {
//reads next visibility matrix window
......@@ -142,32 +163,43 @@ int main(int narg, const char* arg[])
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
}
if (kt >= wreader.getTotalNbWindows()) { // Ca ne devrait pas arriver
cout << "visamm/BUG: kt >= wreader.getTotalNbWindows() -> getting out of the read loop ..." << endl;
break;
}
dateobs=cfdate;
double racourant = P4Coords::RAFromTimeTU(dateobs);
sa_size_t idxra = ram.getMapIndex(racourant);
timevec(kt) = dateobs.TimeDifferenceSeconds(dateobs, datestart);
ravec(kt) = racourant;
if (idxra < 0) continue; // Out of ra-map
// updating count map
racntmap(idxra) += (r_4)1.;
// To compute auto-correlation average in each frequency band
vector< complex<float> > avgac(freqbands.size());
for(size_t i=0; i<avgac.size(); i++) avgac[i]=complex<float>(0.,0.);
// Getting the 4 H auto-correlations
for(size_t i=0; i<4; i++) {
TVector< complex<r_4> > vir = vismtx.Row(KVAC[i]);
vector< complex<float> > vacm = freqbands.getAverage(vir);
for(size_t j=0; j<vacm.size(); j++) {
TArray< complex<r_4> > & arr = vamm[j];
arr(kt, (sa_size_t)i) = vacm[j];
TArray<r_4> & arr = vammac[j];
arr(idxra, (sa_size_t)i) += vacm[j].real();
avgac[j] += vacm[j]; // to compute average auto-correlation
}
}
// Getting the 6 H-H cross-correlations
// Updating the auto-correlation component in the output visibility array
for(size_t j=0; j<freqbands.size(); j++) {
avgac[j] *= complex<r_4>(0.25 , 0.); // to compute the average of the 4 auto-corr
TArray< complex<r_4> > & arr = vamm[j];
arr(idxra, (sa_size_t)0) += avgac[j];
}
// Getting the 6 H-H cross-correlations and updating the corresponding component in the output visibility array
for(size_t i=0; i<KVCXHH.size(); i++) {
TVector< complex<r_4> > vir = vismtx.Row(KVCXHH[i]);
vector< complex<float> > vcxm = freqbands.getAverage(vir);
for(size_t j=0; j<vcxm.size(); j++) {
TArray< complex<r_4> > & arr = vamm[j];
arr(kt, (sa_size_t)i+4) = vcxm[j];
TArray< complex<r_4> > & arrZ = vamm[j];
arrZ(idxra, (sa_size_t)i+4) += vcxm[j];
}
}
dateobs=cfdate;
timevec(kt) = dateobs.TimeDifferenceSeconds(dateobs, datestart);
ravec(kt) = P4Coords::RAFromTimeTU(dateobs);
kt++; cnt++;
......@@ -178,22 +210,59 @@ int main(int narg, const char* arg[])
} // Fin de la boucle de lecture
cout<<"vis2ra/Info: count="<<cnt<<" visimtx read "<<endl;
r_4 cmin, cmax, cmean;
racntmap.MinMax(cmin,cmax);
cmean=racntmap.Sum()/(r_4)racntmap.Size();
cout<<" visamm/Info: Normalising using count map - Count/ Min="<<cmin<<" Max="<<cmax<<" Mean="<<cmean<<endl;
TVector<r_4> norm((sa_size_t)racntmap.Size());
TVector< complex<r_4> > normZ((sa_size_t)racntmap.Size());
norm = (r_4)0.;
normZ = complex<r_4>(0.,0.);
for(sa_size_t i=0; i<racntmap.Size(); i++) {
if (racntmap(i)<1e-9) continue;
r_4 nfac=1./racntmap(i);
norm(i) = nfac; normZ(i) = complex<r_4>(nfac,0.);
}
for(size_t l=0; l<freqbands.size(); l++) {
TArray<r_4> & arr = vammac[l];
TArray< complex<r_4> > & arrZ = vamm[l];
for(sa_size_t j=0; j<arr.SizeY(); j++) {
arr(Range::all(),Range(j),Range::all()).Mul(norm);
}
for(sa_size_t j=0; j<arrZ.SizeY(); j++) {
arrZ(Range::all(),Range(j),Range::all()).Mul(normZ);
}
DVList visaI;
visaI["DATADESC"]="PAON4 AutoCorrelations (1H,2H,3H,4H)";
visaI["DATASET"]=params.inpath5_;
if (fgelevset) visaI["ELEVATION"]=elevation;
visaI["FREQCENTER"]=freqbands.getBand(l).getCentralFreq();
visaI["FREQBAND"]=freqbands.getBand(l).getBandwidth();
arr.Info()=visaI;
visaI["DATADESC"]="PAON4 Auto/Cx HH visibilities AC,1-2,1-3,1-4,2-3,2-4,3-4";
arrZ.Info()=visaI;
}
POutPersist po(outfile);
// --- saving AutoCorr time-frequency maps
cout<<" visamm/Info: Saving " << vamm.size() << " Single freq Visibility array to PPF file "<<outfile<<endl;
for(size_t k=0; k<vamm.size(); k++) { // loop over all the specified frequency bands
TArray< complex<r_4> > & arr = vamm[k];
arr.Info()["FREQCENTER"]=freqbands.getBand(k).getCentralFreq();
arr.Info()["FREQBAND"]=freqbands.getBand(k).getBandwidth();
TArray< complex<r_4> > & arrZ = vamm[k];
TArray<r_4> & arr = vammac[k];
char buff[32];
sprintf(buff,"vamm_%d",(int)k);
cout << " Writing array for FreqBand="<<freqbands.getBand(k)<< " with name="<<buff<<endl;
cout << " Writing Visibility array for FreqBand="<<freqbands.getBand(k)<< " with name="<<buff<<endl;
po << PPFNameTag(buff) << arrZ;
sprintf(buff,"aca_%d",(int)k);
cout << " Writing AutoCorrelation array for FreqBand="<<freqbands.getBand(k)<< " with name="<<buff<<endl;
po << PPFNameTag(buff) << arr;
}
cout << " Writing TimeVec & RAVec arrays ..."<<endl;
po << PPFNameTag("TimeVec") << timevec ;
po << PPFNameTag("RAVec") << ravec ;
po << PPFNameTag("RACount") << racntmap ;
cout << resu; // Update est fait lors du print
......
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