Commit 4ebb56ba authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Ajout programme visamm.cc de calcul de tableaux de visibilites (Time/RA selon...

Ajout programme visamm.cc de calcul de tableaux de visibilites (Time/RA selon x, Auto,CxCorr selon Y) , par bande de frequence, adapte pour l'etape de reconstruction des cartes, Reza 15/08/2018
parent 9d293aa5
......@@ -6,17 +6,17 @@ 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 p4phcor.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h p4phcor.h p4freqselmgr.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o $(OBJ)/p4freqselmgr.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra $(EXE)/p4vdblist $(EXE)/rdthermrfi
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra $(EXE)/p4vdblist $(EXE)/rdthermrfi $(EXE)/visamm
rm -f $(OBJ)/rdvisip4.o $(OBJ)/visi2ntac.o $(OBJ)/visi2dtacx.o $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o $(OBJ)/visiavg.o \
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o $(OBJ)/rdthermrfi.o
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o $(OBJ)/rdthermrfi.o $(OBJ)/visamm.o
rm -f $(MYOLISTHERE)
depend :
......@@ -46,6 +46,9 @@ $(OBJ)/p4gvcor.o : p4gvcor.cc $(MYINCLISTHERE)
######
$(OBJ)/p4phcor.o : p4phcor.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4phcor.o p4phcor.cc
######
$(OBJ)/p4freqselmgr.o : p4freqselmgr.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4freqselmgr.o p4freqselmgr.cc
###############################################################
###### Compilation et link des executables
......@@ -69,6 +72,16 @@ $(EXE)/vis2ra : $(OBJ)/vis2ra.o $(MYOLISTHERE)
$(OBJ)/vis2ra.o : vis2ra.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/vis2ra.o vis2ra.cc
## Calcul de Tableaux de visibilites - un tableau pour chaque frequence, adapte pour le map making
visamm : $(EXE)/visamm
echo '---visamm made'
$(EXE)/visamm : $(OBJ)/visamm.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/visamm $(OBJ)/visamm.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/visamm.o : visamm.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visamm.o visamm.cc
######
## programme de lecture et cartes TFM pour voies thermometre et RFI
rdthermrfi : $(EXE)/rdthermrfi
......
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU
visamm: programme de lecture des fichiers matrices de visibilites
de PAON4, calcul de tableau de visibilites moyennees/filtrees
adaptees pour ingestion par l'etape de map-making
R.Ansari, O. Perdereau - LAL , Aout 2018
---------------------------------------------------------- */
// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <string>
#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h"
#include "datatable.h"
#include "histinit.h"
#include "matharr.h"
#include "timestamp.h"
#include "utilarr.h"
#include "fitsioserver.h"
#include "fiosinit.h"
// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"
// include lecteur de fichiers visibilites
#include "visp4winreader.h"
#include "p4freqselmgr.h"
int Usage(void);
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"
<< " Note that Frequency bands should be specified \n" << endl;
P4AnaParams::UsageOptions();
cout<< endl;
return 1;
}
//----------------------------------------------------
int main(int narg, const char* arg[])
{
// --- Decoding parameters
if( (narg<2) || ((narg>1)&&(strcmp(arg[1],"-h")==0) ) ) return Usage();
FitsIOServerInit();
P4AnaParams params;
params.DecodeArgs(narg, arg);
string outfile = params.outfile_;
if (outfile.length()<1) outfile = "visamm.ppf";
string fitsoutfile = params.fitsoutfile_;
if (fitsoutfile.length()>=1) {
fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
}
if (params.fbands_.size() < 1) {
cout << " visamm/ERROR : No frequency band specified -> STOP \n"<<endl;
return 1;
}
int prtlev = params.prtlev_;
params.Print(cout);
cout <<"visamm/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
// ---
int rc = 0;
try {
HiStatsInitiator _inia;
ResourceUsage resu;
P4RAMapUtil ram(params);
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++) {
cout << "Band["<<i<<"] "<<params.fbands_[i]<<endl;
freqbands.addBand(params.fbands_[i]);
}
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;
}
cout << " List of HH X-cor rows:"<<endl;
for(size_t k=0; k<KVCXHH.size(); k++) {
cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<" ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
}
*/
VisiP4WindowReader wreader(params);
long Imin = wreader.getReader().getSerialFirst();
long Imax = wreader.getReader().getSerialLast();
long Istep = wreader.getReader().getSerialStep();
cout << "visamm/Info: processing visibility matrix serial/sequence number range "
<<Imin<<" <= seq <= " << Imax << " with step="<<Istep<<endl;
cout << " WindowSize="<<wreader.getWindowSize()<<" -> TotalNbWindows="<<wreader.getTotalNbWindows()<<endl;
bool fgok=true;
// Vecteurs des temps et des ra
TVector< double > timevec(wreader.getTotalNbWindows());
TVector< double > ravec(wreader.getTotalNbWindows());
TimeStamp dateobs, cfdate,datestart;
TimeStamp dateorg(2015,1,1,12,0,0.); // Date origine 1 jan 2015
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;
while (fgok) {
//reads next visibility matrix window
fgok = wreader.Shift();
if (!fgok) break;
TMatrix< complex<r_4> > vismtx = wreader.getAverageVisMtx(cfdate);
if (cnt==0) {
// 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;
}
// 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];
}
}
// Getting the 6 H-H cross-correlations
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];
}
}
dateobs=cfdate;
timevec(kt) = dateobs.TimeDifferenceSeconds(dateobs, datestart);
ravec(kt) = P4Coords::RAFromTimeTU(dateobs);
kt++; cnt++;
if ( (cnt>0) && (cnt%50==0) ) {
cout<<"visamm/Info: read/fill VisMtxCount="<<cnt
<<" /Max="<<wreader.getTotalNbWindows()<<" DateObs="<<dateobs<<endl;
}
} // Fin de la boucle de lecture
cout<<"vis2ra/Info: count="<<cnt<<" visimtx read "<<endl;
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();
char buff[32];
sprintf(buff,"vamm_%d",(int)k);
po << PPFNameTag(buff) << arr;
}
po << PPFNameTag("TimeVec") << timevec ;
po << PPFNameTag("RAVec") << ravec ;
cout << resu; // Update est fait lors du print
}
catch (PException& exc) {
cerr << " vis2ra.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
cerr << "\n vis2ra.cc std::exception :"
<< (string)typeid(sex).name() << "\n msg= "
<< sex.what() << endl;
rc = 78;
}
catch (...) {
cerr << " vis2ra.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> vis2ra.cc ------- END ----------- RC=" << rc << endl;
return rc;
}
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