adds the possibility to read template data as fits

Python/viewtfmap.py

@@ -203,43 +203,6 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
             else :
                 plt.savefig(filename+"_"+ names[num]+"_TFM"+ext+".png")
             
-        axcolor = 'lightgoldenrodyellow'
-        axmin = fig.add_axes([0.075, 0.05, 0.35, 0.03], facecolor = axcolor )
-        axmax = fig.add_axes([0.075, 0.02, 0.35, 0.03], facecolor = axcolor)
-#       axrange = fig.add_axes([0.6, 0.05, 0.35, 0.03], facecolor = axcolor)
-#       axcent = fig.add_axes([0.6, 0.02, 0.35, 0.03], facecolor = axcolor)
-
-        smin = Slider(axmin, 'Min', vmin,vmax, valinit=vmin)
-        smax = Slider(axmax, 'Max', vmin,vmax, valinit=vmax)
-        
-        def update(val):
-            print "min=" + str(smin.val)
-            print "max=" + str(smax.val)
-            cbar.norm.vmin=smin.val
-            cbar.norm.vmax=smax.val
-            start = time.time()
-            cbar.draw_all()
-            #elp = time.time()
-            #print str(elp-start)  + " seconds" 
-            im.set_norm(cbar.norm)
-            #elp = time.time()
-            #print str(elp-start)  + " seconds" 
-            cbar.patch.figure.canvas.draw()
-            #elp = time.time()
-            #print "end->"+str(elp-start)  + " seconds" 
-            #cen.set_val( (smax.val + smin.val)/2.)
-            #rng.set_val( (smax.val - smin.val)   )
-            
-            #m.set_clim([smin.val,smax.val])
-            #ig.canvas.draw()
-        
-        smin.on_changed(update)
-        smax.on_changed(update)
-
-        #        class Index(object):
-        #            ind =  cm_cycle.index(cbar.get_cmap().name)
-        #            cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
-        
 
         axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
         bnext = Button(axnext,cm_cycle[cm_index+1] )
@@ -350,7 +313,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
     elif    cmplx:
         return times/3600.,img,freqs,ras
     else :
-        return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev,smin,smax,mskdat  
+        return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev,mskdat  
 
 
  

p4gvcor.cc

@@ -21,6 +21,9 @@
 
 #include "nbtrixx.h"
 #include "cspline.h"
+#include "fitsioserver.h"
+#include "fiosinit.h"
+
 using namespace SOPHYA ;
 using namespace std;
 
@@ -31,7 +34,7 @@ P4gvCor::P4gvCor()
 }
 
 P4gvCor::P4gvCor(string const & param_filename, string const & filt_blind_name)
-: Gt_corr_param_defined_(false), glob_normalisation_defined_(false)
+  : Gt_corr_param_defined_(false), glob_normalisation_defined_(false)
 {
   InitNormalisation();
   setGtCorrection(param_filename, filt_blind_name);
@@ -99,21 +102,48 @@ void P4gvCor::setGtCorrection(string const & param_filename, string const & filt
    throw IOExc("P4gvCor::setGtCorrection() : NO DATA FILE   ");
  }
  
- cout << "P4gvCor::setGtCorrection() - reading PPF fileterd blind channel signal from "<<filt_blind_name<<endl;
- PInPersist pinp(filt_blind_name);
-
- pinp >>  PPFNameTag("TVFilt_Times") >> filt_vtimes;
- pinp >>  PPFNameTag("TVFilt_Values") >> filt_values;
- 
  vector <double> x;
  vector <double> y;
- for (sa_size_t k =0 ; k<filt_vtimes.Size() ; k++){
-   x.push_back(filt_vtimes(k));
-   y.push_back(filt_values(k));
 
+ string ppfstr(".ppf");
+ string ftsstr(".fits");
+ if (filt_blind_name.find(ppfstr  ) != std::string::npos) {
+   std::cout << " ppf found!" << '\n';
+   // lecture de template ds un ppf
+   PInPersist pinp(filt_blind_name);
+   cout << "P4gvCor::setGtCorrection() - reading PPF fileterd blind channel signal from "<<filt_blind_name<<endl;
+   pinp >>  PPFNameTag("TVFilt_Times") >> filt_vtimes;
+   pinp >>  PPFNameTag("TVFilt_Values") >> filt_values;
+   
+   for (sa_size_t k =0 ; k<filt_vtimes.Size() ; k++){
+     x.push_back(filt_vtimes(k));
+     y.push_back(filt_values(k));
+     
+   }
+ } else {
+   if (filt_blind_name.find(ftsstr  ) != std::string::npos) {
+     cout <<" reading in fits file for calib templates "<<endl;
+     FitsInOutFile fior(filt_blind_name,FitsInOutFile::Fits_RO); 
+     fior.MoveAbsToHDU (2) ; // move to hdu 2
+     cout << fior << endl;  
+
+     fior >> filt_vtimes;
+
+     fior >> filt_values;
+
+     for (sa_size_t k =0 ; k<filt_vtimes.Size() ; k++){
+       x.push_back(filt_vtimes(k));
+       y.push_back(filt_values(k));
+     
+     }
+
+   }else{
+     throw ParmError(" P4gvCor::setGtCorrection ERROR : filt_blind_nam -> file type cant be identified  !!" )  ;
+
+   }
  }
  my_interp_.DefinePoints (x,y);
- Gt_corr_param_defined_ = true;
+ Gt_corr_param_defined_ = true; 
 }
 
 
@@ -161,17 +191,17 @@ void P4gvCor::applyGain(TMatrix <complex <float> > & mtx,  TimeStamp const & tms
   vector< double > visi_gains(venc.getTotalNbVisibilities());
 
   //--- computing gains 
-  double vbl = my_interp_(tmstp.ToDays());  
-  
+  double vbl = my_interp_(tmstp.ToDays());  // time in days
+  //cout << " P4gvCor::applyGain at "<<tmstp.ToDays()<<" get "<<vbl<< endl;
 
   for(size_t k=0; k<visi_gains.size(); k++) {
     std::pair<sa_size_t , sa_size_t> pij = venc.getFeedPair(k);
     sa_size_t i=pij.first;
     sa_size_t j=pij.second;
-    double ggi = slopes_[i]*vbl+offsets_[i];
+    double ggi = offsets_[i] + slopes_[i]*vbl;
     double gg=ggi;
     if (i!=j) { 
-      double ggj = slopes_[j]*vbl+offsets_[j];
+      double ggj = offsets_[j] + slopes_[j]*vbl;
       gg=sqrt(ggi*ggj);
     }
     visi_gains[k]=1./gg;
@@ -179,7 +209,7 @@ void P4gvCor::applyGain(TMatrix <complex <float> > & mtx,  TimeStamp const & tms
     if (gg==0.) 
       visi_gains[k]=1. ;
     if (fgdonorm)  visi_gains[k]*= visi_norm_[k];
-    
+    //cout <<" chan "<<k<<" gain " << visi_gains[k]<<endl;
   }
 
   //--- applying gains ...

p4gvcor.h

@@ -55,7 +55,7 @@ protected:
   std::string filt_blind_filename_;
   vector <double> slopes_ ;
   vector <double> offsets_ ;
-  
+  bool read_from_ppf;
   TVector<double> filt_vtimes;
   TVector<double> filt_values;
   SLinInterp1D my_interp_;