Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4

Merge apres modification de trkfit.h .cc , Reza

Python/README

+scripts in this folder : 
+viewtfmap.py (+mynormalize) : reading and plotting from t-f map 
+(or t-ra) fits files written by visiavg
+
+cmp_simddat.py : to plot fitted and input tracks (visi vs time at some freq) used in the fitting of pointing & phase with satellites (and sources) 
+

Python/cmp_simddat.py

+from astropy.io import fits
+import numpy as np
+from matplotlib import pyplot as plt
+
+# use make_fig_scan(filename,numero_antenne,nombre_de_traces)
+# to get plots of all autocor visi comparision between 
+# data and fitted results from trkacxfit
+# this will also plot the difference (in gray, scale on the right)
+# unused tracks are alsor plotted
+#numero_antenne = 1,2,.. pour 1H, 2H,...
+
+# make_fig_scan_xc(filename,cross,nombre_de_traces)
+# plots of all cross-cor visi comparision between 
+# data and fitted results from trkacxfitc (real and imag part) 
+# this will also plot the difference (in gray, scale on the right)
+# for the real part 
+# unused tracks are alsor plotted
+# cross = 1,...,6 pour 1Hx2H, 1Hx3H, ...,3Hx4H
+
+
+names_ac=['1H','2H','3H','4H']
+names_xc=['1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H']
+#cx_1_check_gal.fits
+keynam='HIERARCH NOM_OBJET'
+
+
+
+def myplot_xcor(hdul,i):
+    plt.clf()
+    plt.plot(hdul[i].data,label='Real')
+    plt.plot(hdul[i+1].data,label='Imag')
+    plt.plot(hdul[i+2].data,label='simu(real)')
+    plt.plot(hdul[i+3].data,label='simu(imag)')
+    plt.legend()
+    
+
+def myplot_ac(hdul,i):
+    plt.clf()
+    print hdul[i].header[keynam]
+    print hdul[i+1].header[keynam]
+    plt.plot(hdul[i].data,label='Signal')
+
+    plt.plot(hdul[i+1].data,label='Simu')
+
+    plt.legend()
+    
+def get_num(nant,ntrk,tot):
+    ntims=tot
+    ntrks=2*tot
+    oft=(nant-1)*(ntrks)+ntims
+    print nant 
+    print 'oft ='+str(oft)
+    idt=2*(ntrk-1) # track no starts at 1
+    print "idt="+str(idt)
+    if nant==1 : #antenna no starts at 1
+        print idt
+        
+        idt = idt+ntrk #tims
+        print idt
+    else :
+        idt = idt+oft
+        print idt
+    return idt
+
+def get_num_cx(nant,ntrk,tot):
+    ntims=tot
+    ntrks=4*tot
+    oft=(nant-1)*(ntrks)+ntims
+    print nant 
+    print 'oft ='+str(oft)
+    idt=4*(ntrk-1) # track no starts at 1
+    print idt
+    if nant==1 : #antenna no starts at 1
+        print idt
+        
+        idt = idt+ntrk #tims
+        print idt
+    else :
+        idt = idt+oft
+        print idt
+    return idt
+    
+    
+def myplot_acn(hdul,nant,ntrk,tot):
+    i=get_num(nant,ntrk,tot)
+    plt.clf()
+    print hdul[i].header[keynam]
+    print hdul[i+1].header[keynam]
+    plt.plot(hdul[i].data,label='Signal')
+
+    plt.plot(hdul[i+1].data,label='Simu')
+
+    plt.legend()
+    
+def myplot_xcorn(hdul,nant,ntrk,tot):
+    i=get_num_cx(nant,ntrk,tot)
+    print hdul[i].header[keynam]
+    print hdul[i+1].header[keynam]
+    print hdul[i+2].header[keynam]
+    print hdul[i+3].header[keynam]
+    plt.clf()
+    plt.plot(hdul[i].data,label='Real')
+    plt.plot(hdul[i+1].data,label='Imag')
+    plt.plot(hdul[i+2].data,label='simu(real)',linestyle=":")
+    plt.plot(hdul[i+3].data,label='simu(imag)',linestyle=":")
+    plt.legend()
+
+ 
+def make_fig_scan(ffnam,ant,ntot):
+    hda=fits.open(ffnam)
+    nrows=np.int(np.sqrt(ntot))
+    ncols=ntot/nrows
+    ncols=ncols+ntot-ncols*nrows
+    print str(nrows)+ "  " +str(ncols)
+    #plt.figure(figsize=(20,17))
+    f, axarr = plt.subplots(nrows,ncols,figsize=(ncols*4,nrows*4))
+    axarf=axarr.flatten()
+    f.suptitle('Auto '+names_ac[ant-1])
+    for ip in np.arange(ntot):
+        i=get_num(ant,ip+1,ntot)
+        
+        print hda[i].header[keynam]
+        print hda[i+1].header[keynam]
+        axarf[ip ].plot(hda[i].data,label='Signal')
+        axarf[ip ].plot(hda[i+1].data,label='Simu')        
+        axarf[ip ].legend()
+        vscale=np.max(hda[i].data)-1.
+        vmin=1-0.2*vscale
+        vmax=(1.+vscale)*1.1
+        
+        axarf[ip ].set_ylim((vmin,vmax))
+        axarf[ip].set_title('Track '+str(ip))
+        axarf[ip].set(xlabel='time (mn)', ylabel='signal (AU)')
+        ax2 = axarf[ip].twinx()
+        ax2.plot(hda[i].data-hda[i+1].data,label='Data-Simu',color='gray')
+        ax2.set_ylim((-1,9))
+        ax2.legend(loc=2)
+        
+    plt.tight_layout()
+    plt.subplots_adjust(top=0.9)
+
+def make_fig_scan_xc(ffnam,cross,ntot):
+    hda=fits.open(ffnam)
+    nrows=np.int(np.sqrt(ntot))
+    ncols=ntot/nrows
+    ncols=ncols+ntot-ncols*nrows
+    print str(nrows)+ "  " +str(ncols)
+
+    f, axarr = plt.subplots(nrows,ncols,figsize=(ncols*4,nrows*4))
+    axarf=axarr.flatten()
+    
+    f.suptitle('CrossCor '+names_xc[cross-1])
+    for ip in np.arange(ntot):
+
+        i=get_num_cx(cross,ip+1,ntot)
+        print hda[i].header[keynam]
+        print hda[i+1].header[keynam]
+        #axarf[ip ].plot(hda[i].data,label='Signal')
+        #axarf[ip ].plot(hda[i+1].data,label='Simu')
+        #axarf[ip ].legend()
+        axarf[ip ].plot(hda[i].data,label='Real')
+        axarf[ip ].plot(hda[i+1].data,label='Imag')
+        axarf[ip ].plot(hda[i+2].data,label='simu(real)',linestyle=":")
+        axarf[ip ].plot(hda[i+3].data,label='simu(imag)',linestyle=":")
+        axarf[ip ].legend()
+
+        vmaxr=np.max(hda[i].data)*1.2
+        vminr=np.min(hda[i].data)*1.2
+        vmaxi=np.max(hda[i+1].data)*1.2
+        vmini=np.min(hda[i+1].data)*1.2
+        vmax=np.max((vmaxr,vmaxi))
+        vmin=np.min((vminr,vmini))
+        
+        axarf[ip ].set_ylim((vmin,vmax))
+        axarf[ip].set_title('Track '+str(ip))
+        axarf[ip].set(xlabel='time (mn)', ylabel='signal (AU)')
+        ax2 = axarf[ip].twinx()
+        ax2.plot((hda[i].data-hda[i+2].data)/(vmax-vmin),label='Data-Simu (real)',color='gray')
+        ax2.set_ylim((-1,1))
+        ax2.legend(loc=3)
+        
+    plt.tight_layout()
+    plt.subplots_adjust(top=0.9)
+
+def get_good_xctrks(datac):
+    flist = open(datac ,"r")
+    ntrk =0
+    ngtrk=0
+    goods=[]
+    for line in flist :
+        iok=0
+        if "@trk" in line :
+            print 'found track '
+            ntrk=ntrk+1
+            if "NOCX" in line :
+                goods.append(iok)
+                continue
+            if "NOACX" in line :
+                goods.append(iok)
+                continue
+            ngtrk=ngtrk+1
+            goods.append(1  )
+
+    return ntrk,ngtrk,goods
+def get_good_actrks(datac):
+    flist = open(datac ,"r")
+    ntrk =0
+    ngtrk=0
+    goods=[]
+    for line in flist :
+        iok=0
+        if "@trk" in line :
+            print 'found track '
+            ntrk=ntrk+1
+            if "NOAC" in line :
+                goods.append(iok)
+                continue
+            if "NOACX" in line :
+                goods.append(iok)
+                continue
+            ngtrk=ngtrk+1
+            goods.append(1  )
+
+    return ntrk,ngtrk,goods
+
+def make_fig_scan_xc_sel(ffnam,cross,datac):
+    hda=fits.open(ffnam)
+
+    ntot,nok,goods=get_good_xctrks(datac)
+    print ntot,nok,goods
+    nrows=np.int(np.sqrt(nok))
+    ncols=nok/ nrows
+    ncols=ncols+nok-ncols*nrows
+    print str(nrows)+ "  " +str(ncols)
+
+    f, axarr = plt.subplots(nrows,ncols,figsize=(ncols*4,nrows*4))
+    axarf=axarr.flatten()
+    f.suptitle('CrossCor '+names_xc[cross-1])
+    ipl=0
+    for ip in np.arange(ntot):
+        if goods[ip]==0 :
+            continue
+        i=get_num_cx(cross,ip+1,ntot)
+        print hda[i].header[keynam]
+        print hda[i+1].header[keynam]
+        #axarf[ip ].plot(hda[i].data,label='Signal')
+        #axarf[ip ].plot(hda[i+1].data,label='Simu')
+        #axarf[ip ].legend()
+        axarf[ipl].plot(hda[i].data,label='Real')
+        axarf[ipl].plot(hda[i+1].data,label='Imag')
+        axarf[ipl].plot(hda[i+2].data,label='simu(real)',linestyle=":")
+        axarf[ipl].plot(hda[i+3].data,label='simu(imag)',linestyle=":")
+        axarf[ipl].legend()
+
+        vmaxr=np.max(hda[i].data)*1.2
+        vminr=np.min(hda[i].data)*1.2
+        vmaxi=np.max(hda[i+1].data)*1.2
+        vmini=np.min(hda[i+1].data)*1.2
+        vmax=np.max((vmaxr,vmaxi))
+        vmin=np.min((vminr,vmini))
+        
+        axarf[ipl].set_ylim((vmin,vmax))
+        axarf[ipl].set_title('Track '+str(ip))
+        axarf[ipl].set(xlabel='time (mn)', ylabel='signal (AU)')
+        ax2 = axarf[ipl].twinx()
+        ax2.plot((hda[i].data-hda[i+2].data)/(vmax-vmin),label='Data-Simu (real)',color='gray')
+        ax2.set_ylim((-1,1))
+        ax2.legend(loc=3)
+        ipl=ipl+1
+    plt.tight_layout()
+    plt.subplots_adjust(top=0.9)
+    
+ 
+def make_fig_scan_sel(ffnam,ant, datac   ):
+    ntot,nok,goods=get_good_actrks(datac)
+    print ntot,nok,goods
+    hda=fits.open(ffnam)
+    nrows=np.int(np.sqrt(nok))
+    ncols=nok/ nrows
+    ncols=ncols+nok-ncols*nrows
+    print str(nrows)+ "  " +str(ncols)
+    #plt.figure(figsize=(20,17))
+    f, axarr = plt.subplots(nrows,ncols,figsize=(ncols*4,nrows*4))
+    axarf=axarr.flatten()
+    
+    f.suptitle('Auto '+names_ac[ant-1])
+    ipl=0
+    for ip in np.arange(ntot):
+        if goods[ip]==0 :
+            continue
+        i=get_num(ant,ip+1,ntot)
+        
+        print hda[i].header[keynam]
+        print hda[i+1].header[keynam]
+        axarf[ipl ].plot(hda[i].data,label='Signal')
+        axarf[ipl ].plot(hda[i+1].data,label='Simu')        
+        axarf[ipl ].legend()
+        vscale=np.max(hda[i].data)-1.
+        vmin=1-0.2*vscale
+        vmax=(1.+vscale)*1.1
+        
+        axarf[ipl ].set_ylim((vmin,vmax))
+        axarf[ipl].set_title('Track '+str(ip))
+        axarf[ipl].set(xlabel='time (mn)', ylabel='signal (AU)')
+        ax2 = axarf[ipl].twinx()
+        ax2.plot(hda[i].data-hda[i+1].data,label='Data-Simu',color='gray')
+        ax2.set_ylim((-1,9))
+        ax2.legend(loc=2)
+        ipl=ipl+1
+
+    plt.tight_layout()
+    plt.subplots_adjust(top=0.9)
+

Python/viewtfmap.py

@@ -29,8 +29,8 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
 # the for the cross-correlations - the TF amplitudes are complex : in each case one HDU with real part, one with imaginary part averages
 # and two for their respective variances
 # NB the variance computation is optional 
-# at the end of the fits files, acddition HDUs contain 1D data (in reverse order) : times (center of bins), frequencies and RA
-# (omitted in the case od RA-freq maps 
+# at the end of the fits files, acddition HDUs contain 1D data (in reverse order) : times (center of bins, in  seconds),
+# frequencies (MHz) and RA (hours, omitted in the case of RA-freq maps 
 # 
 
 # --------------------------------

p4autils.cc

@@ -619,9 +619,15 @@ static double  paon_lat = 47.3820  ;   // 47deg 22' 55.1''
 static double paon_lon =  2.19964;     // 2deg 11' 58.7''
 
 #define PAON4NBANT  4
+//---- position des antennes d'apres mesures de geometre 
 static double paon4_posx[PAON4NBANT]={0., -5.994, 4.38, 4.383};
 static double paon4_posy[PAON4NBANT]={0., 0.001, -5.997, 5.996};
 static double paon4_posz[PAON4NBANT]={0., 0., 0., 0.};
+//---- Correction de la position des antennes, determinees par ajustementd e la variation des phases (satellites Galileo)
+//--- Fit des variation des phases avec l'angle zenitale de -35 deg  a +12 deg - Reza + Olivier, Juillet 2019 
+static double paon4_cor_posx[PAON4NBANT]={0., 0., 0., 0.};
+static double paon4_cor_posy[PAON4NBANT]={0., -0.014, +0.019, -0.004};
+static double paon4_cor_posz[PAON4NBANT]={0., -0.057, -0.034, -0.002};
 
 double P4Coords::getLon(){return(paon_lon);}
 double P4Coords::getLat(){return(paon_lat);}
@@ -638,7 +644,16 @@ Vector3d P4Coords::getBaseline(sa_size_t num1, sa_size_t num2)
   if ((num1<1)||(num1>4))  throw ParmError("P4Coords::getBaseline() out of range antenna position num1 <1 or>4 !");
   if ((num2<1)||(num2>4))  throw ParmError("P4Coords::getBaseline() out of range antenna position num2 <1 or>4 !");
   num1--;  num2--;
-  return Vector3d(paon4_posx[num2]-paon4_posx[num1],paon4_posy[num2]-paon4_posy[num1],paon4_posz[num2]-paon4_posz[num1]);
+  // Implementation de la correction des positions d'antenne 
+  //  return Vector3d(paon4_posx[num2]-paon4_posx[num1],paon4_posy[num2]-paon4_posy[num1],paon4_posz[num2]-paon4_posz[num1]);
+  double posx_1=paon4_posx[num1]+paon4_cor_posx[num1];
+  double posy_1=paon4_posy[num1]+paon4_cor_posy[num1];
+  double posz_1=paon4_posz[num1]+paon4_cor_posz[num1];
+  
+  double posx_2=paon4_posx[num2]+paon4_cor_posx[num2];
+  double posy_2=paon4_posy[num2]+paon4_cor_posy[num2];
+  double posz_2=paon4_posz[num2]+paon4_cor_posz[num2];
+  return Vector3d(posx_2-posx_1, posy_2-posy_1, posz_2-posz_1); 
 }