1/ Ecriture du programme de projection selon ra et adaptation de la classe de decodage des

arguments , adaptation makefile et classe utilitaire P4RAMapUtil ajoutee
    Reza + Olivier, 29/01/2018

anap4_exemple.d

@@ -26,6 +26,9 @@ phcor   PhaseCorrParamFileName
 ####  Time - frequency maps
 ## @tfm   TFM_TimeBin  TFM_FreqBin
 @tfm   10 4
+####  Right Ascension maps 
+## @ra RA_center_hours  RA_width_hours RA_resol_minutes 
+@ra 12. 24. 1. 
 ####  PrintLevel
 #   @prtlev   PrintLevel  PrintModulo
 @prtlev  1 20 

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 
 
 #  Define our target list 
-all : rdvisip4  visi2ntac  visi2dtacx visi2tmfreq p4conv2fits msvis2dt visiavg filt_blind ckconvphases
+all : rdvisip4  visi2ntac  visi2dtacx visi2tmfreq p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra 
 
 clean :
 	rm -f $(EXE)/rdvisip4  $(EXE)/visi2ntac  $(EXE)/visi2dtacx  $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
-        $(EXE)/filt_blind $(EXE)/ckconvphases 
+        $(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra
 	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)/filt_blind.o  $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o 
 	rm -f $(MYOLISTHERE)
 
 depend : 
@@ -59,6 +59,17 @@ $(EXE)/visiavg : $(OBJ)/visiavg.o $(MYOLISTHERE)
 $(OBJ)/visiavg.o : visiavg.cc $(MYINCLISTHERE) 
 	$(CXXCOMPILE)  -o $(OBJ)/visiavg.o visiavg.cc
 
+## Calcul de cartes en ra (Right Ascension) a partir de donnees de visibilites
+vis2ra : $(EXE)/vis2ra
+	echo '---vis2ra made'
+
+$(EXE)/vis2ra : $(OBJ)/vis2ra.o $(MYOLISTHERE)
+	$(CXXLINK) -o $(EXE)/vis2ra $(OBJ)/vis2ra.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST) 
+
+$(OBJ)/vis2ra.o : vis2ra.cc $(MYINCLISTHERE) 
+	$(CXXCOMPILE)  -o $(OBJ)/vis2ra.o vis2ra.cc
+
+
 ## filtrage median dun canal
 filt_blind : $(EXE)/filt_blind
 	echo '---filt_blind made'

p4autils.cc

@@ -200,7 +200,8 @@ P4AnaParams::P4AnaParams()
     fgserall_(false),  fgtmsel_(false), inwsz_(1), fgreorderfreq_(false), 
     gain_gnu_file_(""), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""), 
     outfile_("p4a.ppf"), fitsoutfile_(""), prtlev_(0),  prtmodulo_(1), 
-    Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false)
+    Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false),
+    fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.)
 {
   
 }
@@ -213,7 +214,9 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
     phase_corr_param_file_(a.phase_corr_param_file_),  
     outfile_(a.outfile_), fitsoutfile_(a.fitsoutfile_), prtlev_(a.prtlev_), prtmodulo_(a.prtmodulo_), 
     Nmean_(a.Nmean_), fgTFM_(a.fgTFM_), TFMtimebin_(a.TFMtimebin_), TFMfreqbin_(a.TFMfreqbin_), 
-    fgdtable_(a.fgdtable_), fbands_(a.fbands_), lastargs_(a.lastargs_)
+    fgdtable_(a.fgdtable_), fbands_(a.fbands_),
+    fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),
+    lastargs_(a.lastargs_)
 {
 }
 
@@ -228,7 +231,9 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
   outfile_=a.outfile_;  fitsoutfile_=a.fitsoutfile_;
   prtlev_=a.prtlev_;  prtmodulo_=a.prtmodulo_;
   Nmean_=a.Nmean_;  fgTFM_=a.fgTFM_;  TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
-  fgdtable_=a.fgdtable_;   fbands_=a.fbands_;  lastargs_=a.lastargs_;
+  fgdtable_=a.fgdtable_;   fbands_=a.fbands_;
+  fgRA_=a.fgRA_;  RAcenter_=a.RAcenter_;  RAwidth_=a.RAwidth_;   RAresol_mn_=a.RAresol_mn_;
+  lastargs_=a.lastargs_;
   return *this;
 }
 
@@ -308,6 +313,12 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
       sscanf(arg[2],"%lg,%lg",&f0,&df);   fbands_.push_back(P4FreqBand(f0,df));    fgdtable_=true;
       arg+=2; narg-=2;   lastargs_.clear(); 
     }
+    else if (fbo=="-ra") {
+      if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
+      int ja,jb;
+      sscanf(arg[2],"%lg,%lg,%lg",&RAcenter_,&RAwidth_,&RAresol_mn_);  fgRA_=true;
+      arg+=2; narg-=2;  lastargs_.clear(); 
+    }
     else if (fbo=="-prt") {
       if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl;  return -1; }
       sscanf(arg[2],"%d,%d",&prtlev_,&prtmodulo_);
@@ -399,7 +410,11 @@ int P4AnaParams::ReadDatacardFile(const char * dcfilename)
     cnt++;
   }
   if (dc.HasKey("tfm"))  {  // @tfm   TFM_TimeBin  TFM_FreqBin
-    TFMtimebin_=dc.IParam("tfm",0,1);   TFMfreqbin_=dc.IParam("tfm",1,1); 
+    TFMtimebin_=dc.IParam("tfm",0,1);   TFMfreqbin_=dc.IParam("tfm",1,1);  fgTFM_=true;
+    cnt++;
+  }
+  if (dc.HasKey("ra"))  {  // @ra RA_center_hours  RA_width_hours RA_resol_minutes 
+    RAcenter_=dc.DParam("ra",0,12.);   RAwidth_=dc.DParam("ra",1,24.);   RAresol_mn_=dc.DParam("ra",2,1.);  fgRA_=true;
     cnt++;
   }
   if (dc.HasKey("prtlev"))  {  // @prtlev   PrintLevel  PrintModulo 
@@ -418,7 +433,7 @@ int P4AnaParams::UsageOptions()
       << "    [-dcf DataCardFileName ]  [-in5 InputPath_BAO5]   [-in6 InputPath_BAO6]  \n"
       << "    [-inrange Imin,Imax[,Istep] ]  [-inall] [-inseltm YYYY-MM-DDThh:mm:ss duration_minutes] \n"
       << "    [-o outfilename] [-gnu gain_filename] [-nmean N] [-reorderfreq] \n"
-      << "    [-tfm timebin,freqbin]  [-freq f0,df] [-prt level] \n"
+      << "    [-tfm timebin,freqbin] [-ra RAcenterHours,RAwidthHours,RAresolMinutes] [-freq f0,df] [-prt level] \n"
       << " -dcf DatacardFileName : read and decode parameters from datacard file  \n"
       << " -in5/6 InputPath_BAO5/6 : input file path for BAO5 and BAO6 \n"
       << " -inrange Imin,Imax[,Istep] : define range of visibilite matrices to be processed (Imin<=I<=Imax with step Istep) \n"
@@ -435,6 +450,7 @@ int P4AnaParams::UsageOptions()
       << " -nmean N: compute average (and sigma...) using N visibility matrices \n"
       << " -reorderfreq: activate frequency reordering in VisiP4Reader, suitable for FFT firmware \n"
       << " -tfm timebin,freqbin : time-frequency maps, averaging timebin=nb of visi matrices, freqbin frequency components \n"
+      << " -ra RAcenterHours,RAwidthHours,RAresolMinutes: define RA (right ascension) map center,total width (hours), resolution (minutes) \n" 
       << " -freq f0,df : definition of a frequency band centered on f0, with total width df (in MHz) \n"
       << "     several frequency bands can be specified by mutiple -freq options \n"
       << " -prt level,modulo: define print level and print modulo count \n"
@@ -468,6 +484,8 @@ ostream& P4AnaParams::Print(ostream& os) const
   os<<" Gain variation  : parameters in "<<  gain_var_file_<< "\n filtered vectors in "<< gain_blind_filt_file_ <<endl;
   os<<" Phase correction  : parameters in "<<  phase_corr_param_file_ <<endl;
   os<<" Nmean="<<Nmean_<<" TFM: TimeBin="<<TFMtimebin_<<" FreqBin="<<TFMfreqbin_<<endl;
+  os<<" Right Ascension: center_="<<RAcenter_<<" Width="<<RAwidth_<<" (hours) Resolution="<<RAresol_mn_
+    <<"("<<RAcenter_-RAwidth_/2<<"<ra<"<<RAcenter_-RAwidth_/2<<")"<<endl;
   if (fbands_.size()==0) os<<" No frequency band defined ..."<<endl;
   else for(size_t i=0; i<fbands_.size(); i++)  os<<"FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
   return os;
@@ -493,3 +511,43 @@ double P4Coords::RAFromTimeTU(TimeStamp const & tu){
   return (ra);
 
 }
+
+
+//------------------------------------------------------
+//----------------- class P4RAMapUtil ------------------
+//------------------------------------------------------
+P4RAMapUtil::P4RAMapUtil(P4AnaParams const& param)
+  : RAcenter_(param.RAcenter_), RAwidth_(param.RAwidth_), RAresol_mn_(param.RAresol_mn_)
+{
+  if ((RAcenter_<0.)||(RAcenter_>=24.)||(RAwidth_<=1.e-19)||(RAwidth_>24.)||(RAresol_mn_<1.e-19)||(RAresol_mn_>360.))
+    throw ParmError("P4RAMapUtil::P4RAMapUtil(P4AnaParams const&) Bad RAcenter_ or RAwidth_ or RAresol_mn_ parameters");
+  invres_hours_ = 60./RAresol_mn_;
+  size_ = (sa_size_t)(RAwidth_*invres_hours_);
+  size_24hours_ = (sa_size_t)(24.*invres_hours_);
+  cout << " P4RAMapUtil:  Resol(min)="<<RAresol_mn_<<" size="<< size_<<" size for 24 hours="<<size_24hours_<<endl; 
+}
+P4RAMapUtil::P4RAMapUtil(double RAcenter, double RAwidth, double RAresol_mn)
+  : RAcenter_(RAcenter), RAwidth_(RAwidth), RAresol_mn_(RAresol_mn)
+{
+  if ((RAcenter_<0.)||(RAcenter_>=24.)||(RAwidth_<=1.e-19)||(RAwidth_>24.)||(RAresol_mn_<1.e-19)||(RAresol_mn_>360.))
+    throw ParmError("P4RAMapUtil::P4RAMapUtil(RAcenter,...) Bad RAcenter_ or RAwidth_ or RAresol_mn_ parameters");
+  invres_hours_ = 60./RAresol_mn_;
+  size_ = (sa_size_t)(RAwidth_*invres_hours_);
+  size_24hours_ = (sa_size_t)(24.*invres_hours_);
+  cout << " P4RAMapUtil:  Resol(min)="<<RAresol_mn_<<" size="<< size_<<" size for 24 hours="<<size_24hours_<<endl; 
+}
+
+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;
+  if (idx<0) {  // on ajoute 24 heures a ra et on recommece
+    idx += size_24hours_;
+  }
+  else if (idx>=size_) {  // on soustrait 24 heures a ra et on retente
+    idx -= size_24hours_;
+  }
+  if ((idx>=0)&&(idx<size_))  return idx;
+  return -1;
+}

p4autils.h

@@ -206,6 +206,9 @@ public:
   //  definition des bandes de frequences pour les DataTables 
   bool fgdtable_;
   vector<P4FreqBand> fbands_;
+  // Definition de la carte en ra (Right Ascension) - heures decimales, sauf pour la resolution, en minutes 
+  bool fgRA_;   // True , on a defini les limites de la carte en RA 
+  double RAcenter_, RAwidth_, RAresol_mn_;
   // les arguments en fin de ligne de commande 
   vector<string> lastargs_;
 };
@@ -227,4 +230,19 @@ public :
   //! static method Converting a TU time into RA (Right Ascension sky coordinates) in decimal degrees 
   static double RAFromTimeTU(TimeStamp const & tu);
 };
+
+//! Class to facilitate projection into RA map
+class P4RAMapUtil {
+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);
+      
+  double RAcenter_, RAwidth_, RAresol_mn_;
+  double invres_hours_;
+  sa_size_t size_;
+  sa_size_t size_24hours_;
+};
 #endif

vis2ra.cc

+// Utilisation de SOPHYA pour faciliter les tests ...
+#include "sopnamsp.h"
+#include "machdefs.h"
+
+/* ---------------------------------------------------------- 
+   Projet BAORadio/PAON4 - (C) LAL/IRFU  
+
+   visiavg: programme de lecture des fichiers matrices de 
+   visibilites de PAON4, calcul cartes de visibilites en 
+   fonction de ra  (Right Ascension)
+   O. Perdereau, R.Ansari   -  LAL , Janvier 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 sophya mesure ressource CPU/memoire ...
+#include "resusage.h"
+#include "ctimer.h"
+#include "timing.h"
+
+// include lecteur de fichiers visibilites 
+#include "visp4winreader.h"
+
+#include "fitsioserver.h"
+#include "fiosinit.h"
+
+int Usage(void);
+int Usage(void)
+{
+  cout << " --- vis2ra.cc : Reads PPF files produced by mfacq time-frequency \n"
+       << "  and produces maps of visibilies as a function of ra (Right Ascension) \n"
+       << " Usage: vis2ra [-arguments] \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 = "visavg.ppf";
+
+  string fitsoutfile = params.fitsoutfile_;
+  if (fitsoutfile.length()>=1) {
+    fitsoutfile = "!"+fitsoutfile ; // adds '!' ?
+  }
+
+  int deltaIavg = params.TFMtimebin_;
+  sa_size_t TFMfbin = params.TFMfreqbin_;
+  int prtlev = params.prtlev_;
+  bool FgTFMAC = true;
+  bool FgTFMCX = true;
+  string desctfmap;
+  bool FgTFM = params.fgTFM_;   // true -> create time-frequency maps
+
+  params.Print(cout);
+  cout <<"vis2ra/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
+       <<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
+       <<" DeltaIAvg="<<deltaIavg<<"\n"
+       <<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
+  P4RAMapUtil ram(params);
+  
+
+  P4AVisiNumEncoder  visiencod;
+  vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
+  vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
+  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;
+  }
+
+
+  // ---
+  HiStatsInitiator _inia;
+  int rc = 0;
+
+  try {
+    ResourceUsage resu;
+
+
+    VisiP4WindowReader wreader(params);
+
+    long Imin = wreader.getReader().getSerialFirst();
+    long Imax =  wreader.getReader().getSerialLast();
+    long Istep =  wreader.getReader().getSerialStep();
+    cout << "vis2ra/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;
+    // vecteur de noms 
+    vector <string> ext_names;
+    // un vecteur avec les temps 
+    //    TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg); 
+    TMatrix< complex<r_4> > vismtx;
+
+    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;
+    //----- 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;
+
+    //----- 8 auto-corr TimeFrequency maps 
+    vector< TArray< r_4 > > vramac;
+    //----- 8 auto-corr TimeFrequency variance maps 
+    vector< TArray< r_4 > > vramac_sq;
+    
+    sa_size_t ramSX, ramSY;
+
+    while (fgok) {
+      //reads next visibility matrix window 
+      fgok = wreader.Shift();
+      if (!fgok)  break;
+      
+      vismtx = wreader.getAverageVisMtx(cfdate);
+      
+      if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 	
+	ramSX=ram.getMapSize();
+	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;
+	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) ); 
+
+	//allocating 6 Cross-Corr H-H time-frequency maps 
+	cout<<"vis2ra/Info: allocating H-H cross-cor RA-Frequency maps : RA->NX="<<ramSX<<" x Freq->NY="<<ramSY<<endl;
+	for(int k=0; k<6; k++) vram.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_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) );
+
+	// recupere le jour de depart @ 0h
+	datestart = TimeStamp(cfdate.DaysPart(),0.);
+ 
+      }
+      
+      dateobs=cfdate;
+      double racourant = P4Coords::RAFromTimeTU(dateobs);
+      sa_size_t idxra = ram.getMapIndex(racourant);
+      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.;
+      //---- 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> vacsq = vac.SquareElt(); 
+	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 
+	  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 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<r_4>  vcxprsq = (real(vcx)).SquareElt(); 
+	TVector<r_4>  vcxpisq = (imag(vcx)).SquareElt(); 
+	TArray< complex<r_4> > & ramap = vram[k];
+	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();
+	} 
+      }   //----- end of loop over the 6 Xcor 
+      
+      //	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
+      //	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
+      //  ... done
+      
+      cnt++;
+      if ( (cnt>0) && (cnt%50==0) ) {
+	//	cout << " vis2ra/Info: RA-map fill, VisMtxCount="<<cnt<< " /Max="<<Imax<<" DateObs="<<dateobs<<endl;
+      }
+    }   // Fin de la boucle de lecture 
+    cout<<"vis2ra/Info: count="<<cnt<<" visimtx read "<<endl;
+
+    // --- Sauvegarde cartes temps-frequence en fits 
+    //FitsABTWriter * fbtw = NULL;
+    FitsInOutFile  * fos = NULL ;
+
+
+    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];
+      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);
+	  ramap(ix, jy) *= fnorm;
+	  ramap_sq(ix, jy) *= fnorm;
+	}
+      }
+    }  // 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];
+      TArray< r_4 > & ramapsqpi = vram_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
+      
+    if (fitsoutfile.length()>=1){
+      cout << " fitsoutfile" <<fitsoutfile<< endl;
+      fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
+    }
+    
+
+    POutPersist poram(outfile);
+    char bufnam[10];
+    int numkey =0;
+    // --- saving AutoCorr time-frequency maps 
+    cout<<"  vis2ra/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
+    const char* ram_names[8]={"RAM_1H", "RAM_2H", "RAM_3H", "RAM_4H", "RAM_1V", "RAM_2V", "RAM_3V", "RAM_4V"};
+    const char* ramsq_names[8]={"VARRAM_1H", "VARRAM_2H", "VARRAM_3H", "VARRAM_4H", "VARRAM_1V", "VARRAM_2V", "VARRAM_3V", "VARRAM_4V"};
+    for(int k=0; k<8; k++)  {   // loop over the 8 AutoCorr 
+      TArray< r_4 > & ramap = vramac[k];
+      poram << PPFNameTag(ram_names[k]) << ramap;
+
+      if (fos != NULL) {
+	ext_names.push_back(ram_names[k]);
+	(*fos)<<  ramap;
+      }
+      TArray< r_4 > & ramapsq = vramac_sq[k];
+      ramapsq = ramapsq - ramap.SquareElt(); 
+      poram << PPFNameTag(ramsq_names[k]) << ramapsq;
+      if (fos != NULL) {
+	ext_names.push_back(ramsq_names[k]);
+	(*fos)<<  ramapsq;
+      }
+      
+    }
+
+    // --- renormalizing and saving H-H Cross-Corr time-frequency maps 
+    cout<<"  vis2ra/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
+
+    const char* ramCC_names[6]={"RAM_1H2H", "RAM_1H3H", "RAM_1H4H", "RAM_2H3H", "RAM_2H4H", "RAM_3H4H"};
+    const char* vrramCC_names[6]={"RVARRAM_1H2H", "RVARRAM_1H3H", "RVARRAM_1H4H", "RVARRAM_2H3H", "RVARRAM_2H4H", "RVARRAM_3H4H"};
+    const char* viramCC_names[6]={"IVARRAM_1H2H", "IVARRAM_1H3H", "IVARRAM_1H4H", "IVARRAM_2H3H", "IVARRAM_2H4H", "IVARRAM_3H4H"};
+
+    for(int k=0; k<6; k++)  {   // loop over the 6 Xcor 
+      TArray< complex<r_4> > & ramap = vram[k];
+      TArray< r_4 > & ramap_sqpr = vram_rp_sq[k];
+      TArray< r_4 > & ramap_sqpi = vram_ip_sq[k];
+      poram << PPFNameTag(ramCC_names[k]) << ramap;
+      if (fos != NULL) {
+	ext_names.push_back(string(ramCC_names[k])+"_real");
+	(*fos)<<  real(ramap);
+	ext_names.push_back(string(ramCC_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(vrramCC_names[k]) << ramap_sqpr;
+      poram << PPFNameTag(viramCC_names[k]) << ramap_sqpi;
+      if (fos != NULL) {
+	ext_names.push_back(vrramCC_names[k]);
+	(*fos)<<  ramap_sqpr;
+	ext_names.push_back(viramCC_names[k]);
+	(*fos)<<  ramap_sqpi;
+      }
+
+    }
+
+    // --- FIN sauvegarde cartes temps-frequence 
+    //    resu.Update();
+    P4FreqBand myp4fre;
+    TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
+    double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
+    for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
+      avg_freqs(kf) = frbase ;
+
+    poram << PPFNameTag("FreqVec") << avg_freqs;
+    if (fos != NULL) {
+      ext_names.push_back("Frequences");
+      (*fos)<< avg_freqs ;
+      cout << " number of objs in fits "<< ext_names.size() << endl;
+      cout << ext_names << endl;
+      delete(fos);
+    }
+
+    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;
+
+}
+