Ajout programme tfm2dt.cc extraction DataTable Auto-cross et MAJ Makefile, Reza 20/12/2018

makefile

@@ -10,13 +10,13 @@ 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  $(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 visamm trkacfit 
+all : rdvisip4  visi2ntac  visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm trkacfit tfm2dt 
 
 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)/visamm $(EXE)/trkacfit
+        $(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra  $(EXE)/p4vdblist  $(EXE)/rdthermrfi $(EXE)/visamm $(EXE)/trkacfit $(EXE)/tfm2dt
 	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)/visamm.o $(OBJ)/trkacfit
+        $(OBJ)/filt_blind.o  $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o $(OBJ)/rdthermrfi.o $(OBJ)/visamm.o $(OBJ)/trkacfit.o $(OBJ)/tfm2dt.o 
 	rm -f $(MYOLISTHERE)
 
 depend : 
@@ -93,6 +93,16 @@ $(EXE)/trkacfit : $(OBJ)/trkacfit.o $(MYOLISTHERE)
 $(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE) acbeam.h gacfit.h 
 	$(CXXCOMPILE)  -o $(OBJ)/trkacfit.o trkacfit.cc
 
+## programme pour extraire un DataTable a partir des TFM (fichier PPF), input programme trkacfit
+tfm2dt : $(EXE)/tfm2dt
+	echo '---tfm2dt made'
+
+$(EXE)/tfm2dt : $(OBJ)/tfm2dt.o $(MYOLISTHERE)
+	$(CXXLINK) -o $(EXE)/tfm2dt $(OBJ)/tfm2dt.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST) 
+
+$(OBJ)/tfm2dt.o : tfm2dt.cc $(MYINCLISTHERE) 
+	$(CXXCOMPILE)  -o $(OBJ)/tfm2dt.o tfm2dt.cc
+
 
 ######
 ## programme de lecture et cartes TFM pour voies thermometre et RFI  

tfm2dt.cc

+/*--- 
+   PAON4 analysis : 
+   Extracting auto & cross correlation DataTable from TFM maps (input PPF file) 
+
+   Example command to run : 
+  csh> ../AnaPaon4/Objs/trkacfit -D 4.5 -sdec 10 -out A21_fac.txt dt_A21,trk_A21_43057,81500,84500,1278.5  dt_A21,trk_A21_43055,90750,93250,1278.5
+
+ ----------------------------------------------------------*/
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <iostream>
+#include <fstream>
+
+#include <cmath>
+#include <complex>
+
+#include "array.h"
+#include "randr48.h"
+#include "skymap.h"
+#include "samba.h"
+#include "fftwserver.h"
+#include "strutilxx.h"
+
+#include "ctimer.h"
+#include "tarrinit.h"
+#include "skymapinit.h"
+#include "fitsioserver.h"
+
+using namespace std; 
+using namespace SOPHYA; 
+
+
+//--- Global parameters 
+static size_t NB_ANTENNES = 4;     // Number of antenna ( 4 for PAON4 
+
+static string infilename;        // Input PPF Time-Frequancy maps filename 
+static string outfilename;       // output filename for auto/cross datatable (PPF or FITS) 
+static bool fgoutfits=false;     // true -> write output as FITS file 
+
+static double FREQ_MIN, FREQ_MAX;
+static bool fgautonorm = false;   // if true, normalise with the maximum value of each autocorrelation ...
+static bool fgdoVV = false;   // if true, extract V autocor and VxV cross-cor
+
+static int prtlevel=0;  // print level 
+
+static TArray<float> TFM_1, TFM_2, TFM_3, TFM_4;
+static TArray< complex<float> >  TFM_12, TFM_13, TFM_14, TFM_23, TFM_24, TFM_34;
+static TVector<double> TimeVec;
+static TVector<double> RAVec;
+
+
+//-------- Declaration of functions in this file - code after the main, below 
+int decode_args(int narg, char* arg[]);
+int read_p4_tfm(bool fgvv=false);
+DataTable Extract_DT();
+
+//--------------------------------------------------------------
+//---------  Main program -------------
+//--------------------------------------------------------------
+int main (int narg, char* arg[])
+{
+  cout << " ------ tfm2dt.cc : extracting auto/cross visibilities datatable from TFM ppf file  ------- " << endl;
+  int rc = 0;
+  try {
+    SophyaInit();
+    int rcda=decode_args(narg, arg);
+    if (rcda)  return rcda;
+    read_p4_tfm(fgdoVV);
+    DataTable dtvis = Extract_DT();
+
+    size_t ll = outfilename.length();
+    size_t pp = outfilename.rfind('.');
+    if ((pp<ll)&&(outfilename.substr(pp)==".fits"))  {
+      cout << "----- tfm2dt/Info Saving  DataTable to FITS file "<<outfilename<<endl;
+      FitsInOutFile fos(outfilename,FitsInOutFile::Fits_Create);
+      fos << dtvis; 
+    }
+    else {
+      cout << "----- tfm2dt/Info Saving  DataTable to PPF file "<<outfilename<<endl;
+      POutPersist pof(outfilename);
+      pof << dtvis;
+    }
+
+  }  // End of try bloc 
+  catch (PThrowable & exc) {
+    cerr << " tfm2dt.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
+         << " - Msg= " << exc.Msg() << endl;
+    rc = 99;
+  }
+  catch (std::exception & e) {
+    cerr << " tfm2dt.cc: Catched std::exception "  
+         << " - what()= " << e.what() << endl;
+    rc = 98;
+  }
+  catch (...) {
+    ;    cerr << " tfm2dt.cc: some other exception (...) was caught ! " << endl;
+    rc = 97;
+  }
+  cout << " ------------ END OF tfm2dt.cc  (Rc= " << rc << ") ------------ " << endl;
+  return rc;
+}
+
+/* --Fonction-- */
+int decode_args(int narg, char* arg[])
+//! Decode program arguments 
+{
+  bool fglonghelp=false;
+  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  fglonghelp=true;
+  if ((narg<2)||fglonghelp) {
+    cout << "  tfm2dt : Extract visibility datatable from TFM PPF file for a given frequency band  \n"
+	 << "  Usage: tfm2dt [-doVV] [-autonorm] [-prt PrintLevel] InFileName FreqRange OutFileName \n"
+	 << "     -doVV :  extract V autoCor VxV crossCor (default H, HxH) \n" 
+	 << "     -autonorm : Normalise using max value over time of each AutoCorrelation (for the freq. band) \n" 
+	 << "     -prt  PrintLevel : define print level \n"
+	 << "   o InFileName : input PPF file name, produced by visiavg.cc \n" 
+  	 << "   o FreqRange : Frequency range to be averaged, in the form FreqMin,FreqMax (in MHz) \n" 
+	 << "   o OutFileName : Output file name to save the DataTable (.ppf or .fits) \n" << endl;
+    return 1;
+  }
+
+
+  vector<string> lastargs;
+  while (narg>1) {
+    string fbo = arg[1];
+    if (fbo=="-doVV") {  // Extract V auto , VxV cross
+      fgdoVV=true;   arg++;  narg--;  lastargs.clear(); 
+    }
+    else if (fbo=="-autonorm") {  // Normalise using max of auto-cor 
+      fgautonorm=true;   arg++;  narg--;  lastargs.clear(); 
+    }
+    else if (fbo=="-prt") {  // print level 
+      if (narg<2) { cout << "tfm2dt/decode_args missing/bad argument, -h for help " << endl;  return -1; }
+      prtlevel=atoi(arg[2]);   arg+=2; narg-=2;  lastargs.clear(); 
+    }
+    else { arg++; narg--; lastargs.push_back(fbo); }
+  }
+
+  if (lastargs.size() < 3) {
+    cout << "  tfm2dt/decode_args/ERROR : missing input/output specification arguments (tfm2dt -h for usage)  \n";
+    return 9;
+  }
+  infilename = lastargs[0];
+  sscanf(lastargs[1].c_str(),"%lg,%lg",&FREQ_MIN,&FREQ_MAX);
+  outfilename = lastargs[2];
+  cout << " tfm2dt/Info : InFileName= "<<infilename<< " OutFileName= "<<outfilename<<" FreqMin,Max="
+       <<FREQ_MIN<<" , "<<FREQ_MAX<<endl;
+  return 0;
+}
+
+/* --Fonction-- */
+int read_p4_tfm(bool fgvv)
+{
+  cout << "---- tfm2dt/read_p4_tfm() ; reading PAON4 TFM arrays from file "<<infilename<<endl;
+  const char * acnamesH[4]={"TFM_1H","TFM_2H","TFM_3H","TFM_4H"};
+  const char * acnamesV[4]={"TFM_1V","TFM_2V","TFM_3V","TFM_4V"};
+
+  const char * cxnamesH[6]={"TFM_1H2H","TFM_1H3H","TFM_1H4H","TFM_2H3H","TFM_2H4H","TFM_3H4H"};
+  const char * cxnamesV[6]={"TFM_1V2V","TFM_1V3V","TFM_1V4V","TFM_2V3V","TFM_2V4V","TFM_3V4V"};
+
+  const char ** acnames = (fgvv?acnamesV:acnamesH);
+  const char ** cxnames = (fgvv?cxnamesV:cxnamesH); 
+
+  PInPersist pin(infilename);
+
+  cout << "... Reading AutoCor TFM : ";
+  for(size_t i=0; i<4; i++) cout << acnames[i] << " ";  cout << endl;
+  pin >> PPFNameTag(acnames[0]) >> TFM_1;
+  pin >> PPFNameTag(acnames[1]) >> TFM_2;
+  pin >> PPFNameTag(acnames[2]) >> TFM_3;
+  pin >> PPFNameTag(acnames[3]) >> TFM_4;
+
+  cout << "... Reading CrossCor TFM : ";
+  for(size_t i=0; i<6; i++) cout << cxnames[i] << " ";  cout << endl;
+  pin >> PPFNameTag(cxnames[0]) >> TFM_12;
+  pin >> PPFNameTag(cxnames[1]) >> TFM_13;
+  pin >> PPFNameTag(cxnames[2]) >> TFM_14;
+  pin >> PPFNameTag(cxnames[3]) >> TFM_23;
+  pin >> PPFNameTag(cxnames[4]) >> TFM_24;
+  pin >> PPFNameTag(cxnames[5]) >> TFM_34;
+
+  cout << "... Reading TimeVec RAVec : ";
+  pin >> PPFNameTag("TimeVec") >> TimeVec;
+  pin >> PPFNameTag("RAVec") >> RAVec;
+
+  return 0;
+}
+
+/* --Fonction-- */
+DataTable Extract_DT()
+{
+
+  double delta_freq=250./(double)TFM_1.SizeY(); 
+  sa_size_t JMIN = (FREQ_MIN-1250.+0.5)/delta_freq;
+  sa_size_t JMAX = (FREQ_MAX-1250.+0.5)/delta_freq;
+  cout << "---- tfm2dt/Extract_DT/Info: FREQ_MIN,MAX="<<FREQ_MIN<<","<<FREQ_MAX<<" -> JMIN,MAX="<<JMIN<<","<<JMAX<<endl;
+
+  sa_size_t SZX=TFM_1.SizeX(); 
+  Vector V11(SZX); V11=0.; Vector V22(SZX); V22=0.;
+  Vector V33(SZX); V33=0.; Vector V44(SZX); V44=0.; 
+  
+  TVector< complex<double> > V12(SZX); V12= complex<double>(0.,0.);
+  TVector< complex<double> > V13(SZX); V13= complex<double>(0.,0.);
+  TVector< complex<double> > V14(SZX); V14= complex<double>(0.,0.);
+  TVector< complex<double> > V23(SZX); V23= complex<double>(0.,0.);
+  TVector< complex<double> > V24(SZX); V24= complex<double>(0.,0.);
+  TVector< complex<double> > V34(SZX); V34= complex<double>(0.,0.);
+  
+  
+  for(sa_size_t i=0; i<V11.Size(); i++) {
+    for(sa_size_t j=JMIN; j<=JMAX; j++) {
+      V11(i)+=TFM_1(i,j); V22(i)+=TFM_2(i,j); 
+      V33(i)+=TFM_3(i,j); V44(i)+=TFM_4(i,j);
+      V12(i)+=TFM_12(i,j);
+      V13(i)+=TFM_13(i,j);
+      V14(i)+=TFM_14(i,j);
+      V23(i)+=TFM_23(i,j);
+      V24(i)+=TFM_24(i,j);
+      V34(i)+=TFM_34(i,j);
+    }
+  }
+  
+  double fnorm=1./(double)(JMAX-JMIN+1);
+  complex<double> znorm=complex<double>(fnorm,0.);
+  V11 *= fnorm;  V22 *= fnorm;
+  V33 *= fnorm;  V44 *= fnorm;
+  V12 *= znorm;  V13 *= fnorm;  V14 *= znorm; 
+  V23 *= znorm;  V24 *= fnorm;  V34 *= znorm; 
+  
+  double min,max;
+  vector<double> vfmax;
+  
+  V11.MinMax(min,max);   vfmax.push_back(max); 
+  cout << " V11->Max="<<max<<endl;
+  V22.MinMax(min,max);   vfmax.push_back(max); 
+  cout << " V22->Max="<<max<<endl;  
+  V33.MinMax(min,max);   vfmax.push_back(max); 
+  cout << " V33->Max="<<max<<endl;    
+  V44.MinMax(min,max);   vfmax.push_back(max); 
+  cout << " V44->Max="<<max<<endl;  
+  
+  vector<double>  vfnorm(NB_ANTENNES);  
+  for(size_t ii=0; ii<NB_ANTENNES; ii++)  vfnorm[ii]=1.;
+  if (fgautonorm) {
+    cout << "---- tfm2dt/Extract_DT/Info: Normalising using max AutoCor values ..."<<endl;
+    for(size_t ii=0; ii<NB_ANTENNES; ii++)  vfnorm[ii]=vfmax[ii];
+    V11 /= vfnorm[0];  V22 /= vfnorm[1];
+    V33 /= vfnorm[2];  V44 /= vfnorm[3];
+    
+    V12 /= complex<double>(sqrt(vfnorm[0]*vfnorm[1]));
+    V13 /= complex<double>(sqrt(vfnorm[0]*vfnorm[2]));
+    V14 /= complex<double>(sqrt(vfnorm[0]*vfnorm[3]));
+    V23 /= complex<double>(sqrt(vfnorm[1]*vfnorm[2]));
+    V24 /= complex<double>(sqrt(vfnorm[1]*vfnorm[3]));
+    V34 /= complex<double>(sqrt(vfnorm[2]*vfnorm[3]));
+  }
+  
+  DataTable dtvis(128);
+  
+  dtvis.AddDoubleColumn("timesec");
+  dtvis.AddDoubleColumn("ra_hour");
+  dtvis.AddDoubleColumn("V11");
+  dtvis.AddDoubleColumn("V22");
+  dtvis.AddDoubleColumn("V33");
+  dtvis.AddDoubleColumn("V44");
+  dtvis.AddDoubleComplexColumn("V12");
+  dtvis.AddDoubleComplexColumn("V13");
+  dtvis.AddDoubleComplexColumn("V14");
+  dtvis.AddDoubleComplexColumn("V23");
+  dtvis.AddDoubleComplexColumn("V24");
+  dtvis.AddDoubleComplexColumn("V34");
+  
+  DataTableRow  row = dtvis.EmptyRow();
+  for(size_t i=0; i<(size_t)SZX; i++)  {
+    row[0] = TimeVec(i);
+    row[1] = RAVec(i);
+    row[2] = V11(i);   row[3] = V22(i);
+    row[4] = V33(i);   row[5] = V44(i);   
+    row[6] = V12(i);   row[7] = V13(i);
+    row[8] = V14(i);   row[9] = V23(i);
+    row[10] = V24(i);   row[11] = V34(i);
+    dtvis.AddRow(row); 
+  }
+  
+  dtvis.SetShowMinMaxFlag(true);
+  cout<<dtvis;
+  return dtvis;
+}