axes in viewtmap ; add fits check file in trkfit; check input

Python/viewtfmap.py

@@ -70,6 +70,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
     if timin ==0. :
         timin=times[0]/3600.
     timax = max(times)/3600.
+    print '  timin/max = '+str(timin)+', ' +str(timax)
     if ramod :
         timin = min(ras)
         timax= max(ras)
@@ -153,7 +154,8 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
         print vmax,np.asarray(img[300:900,]).max()
 
     if (not noplo) :
-        im = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
+        im = ax1.imshow(img,aspect='auto' , interpolation='none',vmin=vmin,vmax=vmax)
+        #m = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
         if (yrcm) :
             im.set_cmap('YlOrBr_r')
         cbar = plt.colorbar( im ) #,fraction=0.021,pad=.05)
@@ -180,13 +182,13 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
                 ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='left',verticalalignment='bottom',rotation=30)
         ax1.yaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
         ax1.xaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
-
+        
         ax2=ax1.twinx()
-        ax2.yaxis.set_view_interval(len(freqs),0.,ignore=True)
+        ax2.yaxis.set_view_interval(frmax,frmin  ,ignore=True)
         ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
         if not ramod :
             ax3=ax1.twiny()
-            ax3.xaxis.set_view_interval(0.,len(times))
+            ax3.xaxis.set_view_interval( timin,timax  ,ignore=True)            
             ax3.set_xlabel("time bin",x=0.85)
 
         if save:
@@ -254,7 +256,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False,   raplot=False,timeplot=F
         bnext.on_clicked(cm_next)
         bprev.on_clicked(cm_prev)
         
-        plt.tight_layout()
+        #plt.tight_layout()
         plt.subplots_adjust(bottom=0.15)
         plt.show()
 

Satellites/trk2dt.cc

@@ -89,6 +89,11 @@ int main(int narg, const char* arg[])
     DataTableRow  row = dt.EmptyRow();
 
     ifstream is(infile.c_str());
+    if (is.fail()){
+	cout<<" Problem opening file !! "<<endl;
+	throw IOExc(" trk2dt : bad file : "+infile+"  - theck it !!");
+      }
+
     string sline;
     while (!is.eof()) {
       sline = "";
@@ -141,7 +146,7 @@ int main(int narg, const char* arg[])
     }
   }
   catch (PException& exc) {
-    cerr << " rdvisip4.cc catched PException " << exc.Msg() << endl;
+    cerr << " trk2dt.cc catched PException " << exc.Msg() << endl;
     rc = 77;
   }  
   catch (std::exception& sex) {

trkfit.cc

@@ -7,6 +7,8 @@
 #include <iomanip>
 
 #include "pexceptions.h"
+
+
 #include "trkfit.h"
 #include "datacards.h"
 #include "array.h"
@@ -17,6 +19,9 @@
 #include "gcxfitbaseline.h"
 
 #include "p4autils.h"
+#include "fitsioserver.h" 
+using namespace std;
+using namespace SOPHYA;
 
 
 //------------------- Print Level for this file --------------------------
@@ -90,6 +95,10 @@ static int decode_trkcard(string const& key, string const& toks)
   char flnmdata[256], flnmtrk[256];
   char sflags[64];
   double ts,te,freq;
+  flnmdata[0]='\0';
+  flnmtrk[0] ='\0';
+  sflags[0] ='\0';
+  ts=te=freq=-9999.;
   sscanf(toks.c_str(),"%s %lg,%lg %lg %s %s",flnmdata,&ts,&te,&freq,flnmtrk,sflags);
 
   dataflnm_p_->push_back(flnmdata);
@@ -307,7 +316,7 @@ size_t AcxDataSet::ReadData(TrkInputDataSet & tkds)
 	vector< complex<double> > & vcx = v_vcx[ii];
 	v_cxdata[ii].push_back(vcx[k]);
 	double acx=std::abs(vcx[k]);
-	v_cxerr[ii].push_back(0.1*sqrt(acx));
+	v_cxerr[ii].push_back(0.1*sqrt(acx)); // calcul d'erreur, a affiner 
 	if (acx<v_min_cxdata[j][ii])  v_min_cxdata[j][ii]=acx;
 	if (acx>v_max_cxdata[j][ii])  v_max_cxdata[j][ii]=acx;
       }
@@ -649,10 +658,21 @@ int ACxSetFitter::doACfit(string outfilename)
 
 int ACxSetFitter::saveExpectedAC(string outcheckfilename)
 {
+
   if (outcheckfilename.length()<1)  return 1;
   cout << "-----ACxSetFitter::saveExpectedAC() : computing expected signal for fitted params , will be saved to file "
        <<outcheckfilename<<endl;
-  POutPersist pos(outcheckfilename);
+  vector <string> data_names;
+  FitsInOutFile *fos=NULL;
+  POutPersist *pos=NULL;
+  size_t ll = outcheckfilename.length();
+  size_t pp = outcheckfilename.rfind('.');
+  if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits"))  {
+    cout << "----- ACxSetFitter::saveExpectedAC : Saving  to FITS file "<<outcheckfilename<<endl;
+    fos=new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
+  } else {
+    pos = new POutPersist (outcheckfilename);
+  }
   size_t NB_ANTENNES = acxd_.getNbAutoCor();
   size_t NTRK = acxd_.NbTrk();
 
@@ -671,19 +691,56 @@ int ACxSetFitter::saveExpectedAC(string outcheckfilename)
     for(size_t j=0; j<NTRK; j++)  {
       double A = v_A[ii][j];
       double B = v_B[ii][j];
-      Vector signal = macs.getDataSignal(j);
-      sprintf(oname,"ac_%d_%d",(int)ii+1,(int)j+1);
-      pos << PPFNameTag(oname)<<signal;
-      Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
-      sprintf(oname,"simac_%d_%d",(int)ii+1,(int)j+1);      
-      pos << PPFNameTag(oname)<<expsignal;
+
       if (ii==0)  {
 	Vector tmvec = macs.getTimeVec(j);
 	sprintf(oname,"tim_%d",(int)j+1);
-	pos << PPFNameTag(oname)<<tmvec;
+	if (pos!=NULL) {
+	  *pos << PPFNameTag(oname)<<tmvec;
+	}
+	if (fos!=NULL) {
+	  data_names.push_back(oname);
+	  tmvec.Info()["NOM_OBJET"]=oname;
+	  *fos<<tmvec;
+	}
+      }
+
+      Vector signal = macs.getDataSignal(j);
+      sprintf(oname,"ac_%d_%d",(int)ii+1,(int)j+1);
+      if (pos!=NULL) {
+	*pos << PPFNameTag(oname)<<signal;
       }
+      if (fos!=NULL) {
+	data_names.push_back(oname);
+	signal.Info()["NOM_OBJET"]=oname;
+	*fos<<signal;
+      }
+
+      Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
+      sprintf(oname,"simac_%d_%d",(int)ii+1,(int)j+1);  
+      if (pos!=NULL) {
+	*pos << PPFNameTag(oname)<<expsignal;
+      }
+      if (fos!=NULL) {
+        expsignal.Info()["NOM_OBJET"]=oname;
+	data_names.push_back(oname);
+	*fos<<expsignal;
+      }      
     }
   } 
+  if (pos!=NULL) delete(pos);
+  if (fos!=NULL) {
+    //*fos<<data_names;
+    //TVector <string> tata(data_names);
+    //*fos<<tata ;
+    
+    cout<<" number of extenstions "<<data_names.size()<<endl;
+    for (int k=0 ; k<data_names.size() ; k++){
+      cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
+    }
+
+    delete(fos);
+  }
   return 0;
 }
 
@@ -919,7 +976,19 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac, bool fg_B0, bool fg
 int ACxSetFitter::saveExpectedCx(string outcheckfilename)
 {
   cout << "ACxSetFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
-  POutPersist pox(outcheckfilename);
+  size_t ll = outcheckfilename.length();
+  size_t pp = outcheckfilename.rfind('.');
+  FitsInOutFile *fox=NULL;
+  POutPersist *pox=NULL;
+  vector <string> data_names;
+  
+  if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits"))  {
+    cout << "----- ACxSetFitter::saveExpectedCx : Saving  to FITS file "<<outcheckfilename<<endl;
+    fox=new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
+  } else {
+    pox = new POutPersist (outcheckfilename);
+  }
+
   size_t NB_CXCORS=acxd_.getNbCrossCor();
   size_t NTRK = acxd_.NbTrk();
 
@@ -930,21 +999,69 @@ int ACxSetFitter::saveExpectedCx(string outcheckfilename)
     MyCxSignal cxsig( acxd_.v_time_data, acxd_.vv_cxdata, acxd_.vv_cxerr, acxd_.v_freqs, acxd_.v_noCx, 
 		      tks_.v_interp_elev, tks_.v_interp_sazim, cxbeam, ii);
     for(size_t j=0; j<NTRK; j++) {
-      TVector< complex<double> >  signal = cxsig.getDataSignal(j);
+
+      if (ii==0)  {
+	Vector tmvec = cxsig.getTimeVec(j);
+	sprintf(oname,"tim_%d",(int)j+1);
+	if (pox!=NULL)
+	  *pox << PPFNameTag(oname)<<tmvec;
+	if (fox!=NULL){
+	  data_names.push_back(oname);
+	  tmvec.Info()["NOM_OBJET"]=oname;
+	  *fox << tmvec;
+	}
+      }
+
+      TVector < complex<double> >  signal = cxsig.getDataSignal(j);
       sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
-      pox << PPFNameTag(oname)<<signal;
+
+      if (pox!=NULL)
+	*pox << PPFNameTag(oname)<<signal;
+      if (fox!=NULL){
+	sprintf(oname,"cx_%d_%d_R",(int)ii+1,(int)j+1);
+	data_names.push_back(oname);
+	Vector tmpv=real(signal);
+	tmpv.Info()["NOM_OBJET"]=oname;
+ 	*fox << tmpv ;
+	sprintf(oname,"cx_%d_%d_I",(int)ii+1,(int)j+1);
+	data_names.push_back(oname);
+	tmpv=imag(signal);
+	tmpv.Info()["NOM_OBJET"]=oname;
+ 	*fox << tmpv ;
+      }
       //DBG      cout << " *DBG* getPhase4Freq() phi0="<<acxd_.v_phi_0[ii]<<" a_phi="<<acxd_.v_a_phi[ii]<<" freq="<<acxd_.v_freqs[j]<<endl;
       double phase=cxsig.getPhase4Freq(acxd_.v_phi_0[ii],acxd_.v_a_phi[ii],acxd_.v_freqs[j]);
       TVector< complex<double> >  expsignal = cxsig.getExpectedSignal(j, phase, v_Acx[ii][j]);
       sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
-      pox << PPFNameTag(oname)<<expsignal;
-      if (ii==0)  {
-	Vector tmvec = cxsig.getTimeVec(j);
-	sprintf(oname,"tim_%d",(int)j+1);
-	pox << PPFNameTag(oname)<<tmvec;
+      if (pox!=NULL)
+	*pox << PPFNameTag(oname)<<expsignal;
+      if (fox!=NULL){
+	sprintf(oname,"simcx_%d_%d_R",(int)ii+1,(int)j+1);
+	data_names.push_back(oname);
+	Vector tmpv = real(expsignal);
+	tmpv.Info()["NOM_OBJET"]=oname;
+ 	*fox << tmpv ;
+
+	sprintf(oname,"simcx_%d_%d_I",(int)ii+1,(int)j+1);
+	data_names.push_back(oname);
+	tmpv = imag(expsignal);
+	tmpv.Info()["NOM_OBJET"]=oname;
+
+ 	*fox << tmpv ;
       }
     }
   }
+  if (pox!=NULL) delete(pox);
+  if (fox!=NULL) {
+
+    cout<<" number of extensions "<<data_names.size()<<endl;
+    for (int k=0 ; k<data_names.size() ; k++){
+      cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
+    }
+    
+    delete(fox);
+
+  }
   return 0;
 }  
 
@@ -1272,7 +1389,18 @@ int CxBaselineFitter::doCheck(int fg_fixebaseline)
 int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
 {
   cout << "CxBaselineFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
-  POutPersist pox(outcheckfilename);
+  FitsInOutFile *fox=NULL;
+  POutPersist *pox=NULL;
+  vector <string> data_names;
+  size_t ll = outcheckfilename.length();
+  size_t pp = outcheckfilename.rfind('.');
+  if ((pp<ll)&&(outcheckfilename.substr(pp)==".fits"))  {
+    cout << "----- CxBaselineFitter::saveExpectedCx : Saving  to FITS file "<<outcheckfilename<<endl;
+    fox = new FitsInOutFile (outcheckfilename,FitsInOutFile::Fits_Create);
+  } else {
+    pox = new POutPersist (outcheckfilename);
+  }
+  
   size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
   char oname[48];
 
@@ -1283,16 +1411,61 @@ int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
     for(size_t j=0; j<v_acxd[i].NbTrk(); j++) {
       Vector tmvec = cxsigb.getTimeVec(i,j);
       sprintf(oname,"tim_%d_%d",(int)j+1,(int)i+1);
-      pox << PPFNameTag(oname)<<tmvec;
+      if (pox!=NULL)
+	*pox << PPFNameTag(oname)<<tmvec;
+      if (fox!=NULL){
+	data_names.push_back(oname);
+	tmvec.Info()["NOM_OBJET"]=oname;
+	*fox << tmvec;
+      }
       for(size_t kcx=0; kcx<NB_CXCORS; kcx++) {
 	TVector< complex<double> >  signal = cxsigb.getDataSignal(i,j,kcx);
 	sprintf(oname,"cx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
-	pox << PPFNameTag(oname)<<signal;
+	if (pox!=NULL){
+	  *pox << PPFNameTag(oname)<<signal;
+	}
+	if (fox!=NULL){
+	  sprintf(oname,"cx_%d_%d_%d_R",(int)kcx+1,(int)j+1,(int)i+1);
+	  data_names.push_back(oname);
+	  Vector vtmp=real(signal);
+	  vtmp.Info()["NOM_OBJET"]=oname;
+	  *fox << vtmp ;
+	  sprintf(oname,"cx_%d_%d_%d_I",(int)kcx+1,(int)j+1,(int)i+1);
+	  data_names.push_back(oname);
+	  vtmp=imag(signal);	  
+	  vtmp.Info()["NOM_OBJET"]=oname;
+	  *fox << imag(signal);
+	}
+	
+
 	TVector< complex<double> >  expsignal = cxsigb.getExpectedSignal(i,j,kcx,bestfitparam);
 	sprintf(oname,"simcx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
-	pox << PPFNameTag(oname)<<expsignal;
+	if (pox!=NULL){
+	  *pox << PPFNameTag(oname)<<expsignal;
+	}
+	if (fox!=NULL){
+	  sprintf(oname,"simcx_%d_%d_%d_R",(int)kcx+1,(int)j+1,(int)i+1);
+	  data_names.push_back(oname);
+	  Vector vtmp=real(expsignal);
+	  vtmp.Info()["NOM_OBJET"]=oname;
+	  *fox << vtmp;
+	  sprintf(oname,"simcx_%d_%d_%d_I",(int)kcx+1,(int)j+1,(int)i+1);
+	  data_names.push_back(oname);
+	  vtmp=imag(expsignal);
+	  vtmp.Info()["NOM_OBJET"]=oname;
+	  *fox << vtmp ;
+	}
+
       }
     }
   }
+  if (pox!=NULL) delete(pox);
+  if (fox!=NULL) {
+    cout<<" number of extenstions "<<data_names.size()<<endl;
+    for (int k=0 ; k<data_names.size() ; k++){
+      cout<<" no "<<k<<" ->"<<data_names[k]<<endl;
+    }
+    delete(fox);
+  }
   return 0;
 }