Gitlab is now running v13.9.0 - More info -> here <-

Commit 7155929a authored by Reza  ANSARI's avatar Reza ANSARI

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

  Merge protection sur RA avec modifs Olivier ds p4autils.cc - Reza
parents 9fefb5a6 13d6684d
......@@ -395,12 +395,13 @@ r_8 Beam::Value(r_8 const xp[], r_8 const* Par) {
} else {
#if TYPE_BEAM == 1
lob = exp(-0.220544*x*x); //gaussian equivalent of J1 formula normalise a 1
//BUG JEC 20/9/19 the magic coeff was determined on spherical_bessel fonction and not the Bessel function lob = exp(-0.220544*x*x); //gaussian equivalent of J1 formula normalise a 1
lob = exp(-0.280092*x*x); //gaussian equivalent of J1 formula normalise a 1
#else
if(x==0.) {
lob = 1.;
} else {
lob = 3. * boost::math::cyl_bessel_j(1, x)/x;
lob = 2. * boost::math::bessel_j(1, x)/x;
lob *= lob;
}
#endif
......
......@@ -2,11 +2,13 @@ import math as math
import sys, platform, os
import numpy as np
from astropy.io import fits
import matplotlib as matplotlib
from matplotlib import pyplot as plt
import itertools as itert
from skimage import data, img_as_float
from skimage import exposure
from matplotlib.ticker import StrMethodFormatter
from matplotlib.ticker import FormatStrFormatter
import mynormalize as mynorm
from matplotlib.widgets import Slider, Button, RadioButtons
import time
......@@ -14,14 +16,23 @@ from colorsys import hsv_to_rgb, hls_to_rgb
# Choices: ['auto', 'gtk', 'gtk3', 'inline', 'nbagg', 'notebook', 'osx', 'qt', 'qt4', 'qt5', 'tk', 'wx']
# values of the freqs for the slices
rfreq = [1381.,1390.,1420.5,1430.]
rfreq = [1347.,1390.,1420.5,1430.]
#rfreq = [1420.30,1420.35,1420.4,1420.45,1420.5,1420.55,1420.6,1420.65,1420.7,1420.75]
rfreq = [1420.30,1420.35,1420.4,1420.45,1420.5,1420.55,1420.6,1420.65,1420.7,1420.75]
rfreq=[1420.,1420.4,1420.8,1421.2,1421.6,1421.8,1422.,1422.25,1422.5]
rfreq=[1420.4,1422.8,1421.8,1421.06,1419.7,1419.0]
normf=[1325.,1465.] # "normal frequencies"
names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H',
'1Vx2V','1Vx3V','1Vx4V','2Vx3V','2Vx4V','3Vx4V','1Hx1V','1Hx2V','1Hx3V','1Hx4V','2Hx1V','2Hx2V','2Hx3V','2Hx4V',
'3Hx1V','3Hx2V','3Hx3V','3Hx4V','4Hx1V','4Hx2V','4Hx3V','4Hx4V']
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 20}
matplotlib.rc('font', **font)
# viewtfmap reads and returns data stored in the TF maps fits files build by visiavg (code in AnaPAON4)
# all TF maps are stored as hdu in the fits file, with some additional data at the end.
......@@ -56,7 +67,8 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
# hasvar : to handle the cas where there are no variances
# ramod : to handle the case where the time-freq maps are ra-frequency maps [deprecated]
# mask : if true, also looks for a 'mask' file (same name+'_mask') to mask RFIs and satellites
# yrcm : to use another color map (
# yrcm : to use another color map
# usecm : choice of colormap
# raplot : also shows slices at some frequencies [see array rfreq above] as a fct of RA [deprecated - don't ew
# timeplot : also shows slices at soem frequencies as a fct of time [deprecated - don't use]
......@@ -76,11 +88,11 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
#
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False,
arg=False,cmplx=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True,
ramod=False,plras=[],nmras=[],mask=False,logsc=False,snratio=False):
arg=False,cmplx=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True,submed=False,subhmed=False,
ramod=False,plras=[],nmras=[],mask=False,logsc=False,snratio=False,usecm="",extaxes=False,choose_cmap=False,nau=True) :
if verb :
print len(cor)
print (len(cor))
filename = folder
if noplo :
clean=False
......@@ -113,14 +125,14 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
i=num # devrait etre juste la pour les sorties des voies th et rfi
if (verb) :
print hdulist[i].header
print( hdulist[i].header)
# get the bin indexes for the slices & the "normal range"
bfreq = [ (np.abs(freqs-rf)).argmin() for rf in rfreq ]
bfnor = [(np.abs(freqs-rf)).argmin() for rf in normf ]
if (verb) :
print ' number of hdus ' + str(len(hdulist))
print (' number of hdus ' + str(len(hdulist)))
nfre = np.size(freqs)
......@@ -129,7 +141,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)
print (' timin/max = '+str(timin)+', ' +str(timax))
if ramod :
timin = min(ras)
timax= max(ras)
......@@ -148,13 +160,32 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if(cimag):
img = hdulist[ i+1].data
if (verb) :
print 'Image shape :'
print np.shape(img)
# module
print ('Image shape :')
print (np.shape(img))
# module
meds=np.zeros(0)
if submed :
meds=np.median(img,axis=1)
print (len(meds))
print (np.shape(img))
for idx in np.arange(np.shape(img)[0]):
img[idx]=img[idx]-meds[idx]
if subhmed :
meds=np.median(img,axis=0)
print (len(meds))
print (np.shape(img))
for idx in np.arange(np.shape(img)[1]):
img[:,idx]=img[:,idx]-meds[idx]
if mod :
if num>7 :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
img = np.sqrt(img**2+imgi**2)
if len(cor) ==2 : # for correlation coeff
imca = hdulist[cor[0]].data
......@@ -164,17 +195,25 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if arg :
if num>7 :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
img = np.arctan2(imgi,img) * 180. / np.pi
mvmin=-180.
mvmax=180.
mtitle= names[num]+" arg.(deg)"
else :
print "arg demande mais auto selectionne : ERREUR "
print ("arg demande mais auto selectionne : ERREUR ")
exit
mskdat=img*0.
# complex number representation
if cmplx :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
zmod = np.sqrt(img**2+imgi**2)
arg = np.arctan2(imgi,img)
......@@ -182,9 +221,9 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
aa = (aa+0.5)%1
bb = 1./(abs(zmod)**0.3 +1.)
bb = -bb+1.
print np.shape(aa)
cc = np.ndarray(np.shape(aa))
print np.shape(cc)
cc.fill(1.)
iimg =[]
for ii in np.arange(np.shape(aa)[0]):
......@@ -197,15 +236,16 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if mask : # reads mask file
if cmplx :
print "mask & cmplx not yet compatible "
print ("mask & cmplx not yet compatible ")
exit
msknam = filename +'_masks.fits'
print msknam+ " " + str(num)
hdmskl = fits.open(msknam)
mskdat = np.fabs(hdmskl[num+1].data-1.)
img = img*mskdat
filename=filename+"_masked"
if logsc :
filename=filename+" (log)"
if mod :
......@@ -224,48 +264,55 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if mvmax !=0. :
vmax = mvmax
if verb :
print vmin, np.asarray(img[300:900,]).min()
print vmax,np.asarray(img[300:900,]).max()
print (vmin, np.asarray(img[300:900,]).min())
print (vmax,np.asarray(img[300:900,]).max())
if (not noplo) :
# plots TF image - dual axes with bin no and values
im = ax1.imshow(img,aspect='auto' , interpolation='none',vmin=vmin,vmax=vmax)
im = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), 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')
im.set_cmap(usecm)
cbar = plt.colorbar( im ) #,fraction=0.021,pad=.05)
cbar.set_norm(mynorm.MyNormalize(vmin=vmin,vmax=vmax,stretch='linear'))
cm_cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
cm_index = cm_cycle.index(cbar.get_cmap().name)
if nau :
cbar.ax.set_ylabel('signal (NAU) ', rotation=270,labelpad=20)
else:
cbar.ax.set_ylabel('signal ', rotation=270,labelpad=20)
plt.show()
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
plt.title(mtitle+ " - " + names[num] ,y=1.08)
plt.ylabel("Frequency (Mhz)")
if ramod:
plt.xlabel("RA (h)")
else:
plt.xlabel("Temps (h)")
plt.xlabel("Time (h)")
if ramod :
for rap in plras:
ax1.plot([rap,rap],[frmax,frmin],color='red',linestyle='--')
#ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='center',verticalalignment='bottom')
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(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( timin,timax ,ignore=True)
ax3.set_xlabel("time bin",x=0.85)
if extaxes :
ax2=ax1.twinx()
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( timin,timax ,ignore=True)
#ax3.set_xlabel("time bin",x=0.85)
plt.tight_layout()
if save: # saves png with some name
ext=""
if creal :
......@@ -280,13 +327,8 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
else :
plt.savefig(filename+"_"+ names[num]+"_TFM"+ext+".png")
# ------------ THIS IS TO SWITCH FROM ONE CMAP TO THE OTHER
axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
bnext = Button(axnext,cm_cycle[cm_index+1] )
axprev = fig.add_axes([0.7, 0.04, 0.05, 0.05])
bprev = Button(axprev,cm_cycle[cm_index-1] )
axcur = fig.add_axes([0.65, 0.04, 0.05, 0.05], )
bcur = Button(axcur,cm_cycle[cm_index])
# ------------ THIS IS TO SWITCH FROM ONE CMAP TO THE OTHER
def cm_next(val):
cm_cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
cm_index = cm_cycle.index(cbar.get_cmap().name)
......@@ -327,9 +369,16 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
indnext=0
bnext.label.set_text(cm_cycle[indnext])
bprev.label.set_text(cm_cycle[indprev])
if choose_cmap :
axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
bnext = Button(axnext,cm_cycle[cm_index+1] )
axprev = fig.add_axes([0.7, 0.04, 0.05, 0.05])
bprev = Button(axprev,cm_cycle[cm_index-1] )
axcur = fig.add_axes([0.65, 0.04, 0.05, 0.05], )
bcur = Button(axcur,cm_cycle[cm_index])
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
#----------------- END CMAP SWITCH
#plt.tight_layout()
......@@ -340,14 +389,21 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if raplot : # does plots of signal vs RA at some frequencies
fig,ax1=plt.subplots(figsize=winwid)
[ plt.plot(ras, img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
#[ plt.plot(ras, img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
# [ plt.scatter(np.mod(ras,24.)[np.argsort(np.mod(ras,24.))], img[b,np.argsort(np.mod(ras,24.))],label="F=%6.2f MHz" %(freqs[b]) ,marker='.', s=.5) for b in bfreq ]
[ plt.plot(np.mod(ras,24.)[np.argsort(np.mod(ras,24.))], img[b,np.argsort(np.mod(ras,24.))],label="F=%6.2f MHz" %(freqs[b])) for b in bfreq ]
plt.xlabel("RA (hour)")
plt.ylabel("Signal ")
plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
if( num<8) :
plt.ylim((.95,1.45))
else:
plt.ylim((-.05,.05))
plt.xlim((0.,24.))
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
......@@ -364,29 +420,35 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if timeplot : # does plots of sig vs time at some frequencies
fig,ax1=plt.subplots(figsize=winwid)
[ plt.plot(times/3600., img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
plt.xlabel('Time (TU, hour)')
plt.ylabel("Signal ")
plt.xticks(fontsize=28)
plt.yticks(fontsize=28)
#plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')) # No decimal places
#plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%g'))
plt.gca().ticklabel_format(axis='y', style='sci', scilimits=(5,5))
plt.xlabel('Time (TU, hour)',size=28)
plt.ylabel("Signal ",size=28)
plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
#legend with thicker linewidth
leg = plt.legend()
leg = plt.legend(prop={'size': 28},loc=1 )
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
plt.ylim((mvmin,mvmax))
plt.tight_layout()
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
plt.title(mtitle+ " - " + names[num] ,y=1.08,size=28)
plt.tight_layout()
if save:
if mask :
plt.savefig(filename+"_"+ names[num]+"_TFM_mask_freq.png")
else:
plt.savefig(filename+"_"+ names[num]+"_TFM_freq.png")
if (verb) :
print times.size
if (verb) :
print np.shape(img)
print (times.size())
print (np.shape(img))
#return values
......@@ -395,7 +457,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,mskdat
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,mskdat
......
......@@ -141,14 +141,15 @@ public:
double Value(UnitVector const& uvo)
{
// circular beam response
double alp=acos(bdir_.Psc(uvo))*DoL_;
// remplacer acos(...) par sin(acos(...)) ds l'expression ci-dessus pour avoir
double alp=acos(bdir_.Psc(uvo))*DoL_;
// remplacer acos(...) par sin(acos(...)) ds l'expression ci-dessus pour avoir l'expression exacte qui est sin(theta)
if (fggauss_) {
// Reza / Sep 2019 : le facteur 2.61854 assure la meme largeur FWHM entre gaussian et Bessel
// On pose : a = D/lambda * gamma ( ou D/lambda * sin(gamma)
// Beam Bessel : beamJ(a) = ( 2 J1(pi a) / (pi a) )^2 Beam Gaussian : beamG(a) = Exp[ - K a^2 ]
// BeamJ(a=1.238098) = 0.5 et beamG(a=1.238098) = 0.5 pour K = 2.61854
// K=2.61854 au lieu de K=2.1767 (valeur de JEC , mais avec SphericalBesselJ[1...]
// Jean-Eric donne 0.28 pi^2 = 2.76 avec un fit de BesselJ
//---- return ( exp(-2.1767*alp*alp)*NormFac_ ); // Jiao : Check the factor 2 in the exp( ) -> exp(-2.*alp*alp)
return ( exp(-2.61854*alp*alp)*NormFac_ ); // Jiao : Check the factor 2 in the exp( ) -> exp(-2.*alp*alp)
}
......
......@@ -357,11 +357,13 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
else if (fbo=="-freq") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
double f0,df;
cout << " selecting freqs "<<endl;
sscanf(arg[2],"%lg,%lg",&f0,&df); fbands_.push_back(P4FreqBand(f0,df)); fgdtable_=true;
arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-killfreq") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
cout << " killing freqs "<<endl;
double f0,df;
sscanf(arg[2],"%lg,%lg",&f0,&df); kill_fbands_.push_back(P4FreqBand(f0,df));
arg+=2; narg-=2; lastargs_.clear();
......@@ -598,8 +600,8 @@ ostream& P4AnaParams::Print(ostream& os) const
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;
if (kill_fbands_.size()==0) {
else for(size_t i=0; i<fbands_.size(); i++) os<<"Selected FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
if (kill_fbands_.size()>0) {
os<<" Frequency bands to mask/kill defined :"<<endl;
for(size_t i=0; i<fbands_.size(); i++) os<<"KillFreqBand["<<i<<"]: "<<fbands_[i]<<endl;
}
......
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2017
Projet BAORadio/PAON4 - (C) LAL 2008-2017
Classes et fonctions utilitaires pour les programmes d'analyse PAON4
R. Ansari, C. Magneville Fevrier-Mars 2017
......
......@@ -87,7 +87,8 @@ int main(int narg, const char* arg[])
cout <<"visiavg/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<" DeltaIAvg="<<deltaIavg<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
<< " fband # "<<params.fbands_.size()
<<" outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
if (!FgTFM) {
cout<<" visiavg/parameter error : specify Time-Frequency map parameter with -tfm "<<endl;
......@@ -131,7 +132,7 @@ int main(int narg, const char* arg[])
bool fgok=true;
// vecteur de noms
vector <string> ext_names;
// un vecteur avec les temps
// un vecteur avec les temps & un avec RA
TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg);
TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg);
TMatrix< complex<r_4> > vismtx;
......@@ -190,6 +191,12 @@ int main(int narg, const char* arg[])
sa_size_t tfmSX, tfmSY;
long serialnum;
size_t rdidx, lastrdidx=0;
double freqmin = P4FreqBand::freqstart_;
double freqmax = P4FreqBand::freqend_;
sa_size_t ifrmin = 0 ;
sa_size_t ifrmax = P4FreqBand::nbnufft_ - 1 ; // ne pas oublier le -1 !
while (fgok) {
//reads next visibility matrix window
fgok = wreader.Shift(rdidx, serialnum);
......@@ -203,39 +210,58 @@ int main(int narg, const char* arg[])
lastrdidx=rdidx;
vismtx = wreader.getAverageVisMtx(cfdate);
///cout << vismtx << endl;
sa_size_t NTFCols = vismtx.NCols(); // bin en freq a ecrire
//cout << "nbre de cols "<<NTFCols << " " << ifrmin<< " "<< ifrmax << endl;
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
// default
if ( params.fbands_.size()>0 ){
sa_size_t XtfmSY = (sa_size_t) floor(params.fbands_[0].df_/params.fbands_[0].deltanufft_ /TFMfbin ) ;
cout << " taille finale en freq "<< XtfmSY << " bin "<< TFMfbin <<endl;
tfmSY=XtfmSY;
ifrmin = (sa_size_t) params.fbands_[0].jfmin_ ;
ifrmax = (sa_size_t) params.fbands_[0].jfmax_ ;
freqmin = params.fbands_[0].getP4Frequency(ifrmin);
freqmax = params.fbands_[0].getP4Frequency(ifrmax);
NTFCols = ifrmax - ifrmin +1 ; //## CHECK THIS
cout << " fmin,max "<< freqmin<<","<< freqmax<<" ident "<<ifrmin<<" "<< ifrmax << endl;
}
acsum.SetSize(8, NTFCols);
if (params.doSigma_) // for computing Sigma, if required
acsum_sq.SetSize(8, vismtx.NCols());
acsum_sq.SetSize(8, NTFCols);
cxsum.SetSize(6, vismtx.NCols());
cxsum.SetSize(6, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
cxsum_sq_rp.SetSize(6, NTFCols);
cxsum_sq_ip.SetSize(6, NTFCols);
}
// VV if needed
if (params.doVV_){
cxsum_vv.SetSize(6, vismtx.NCols());
cxsum_vv.SetSize(6, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
cxsum_vv_sq_rp.SetSize(6, NTFCols);
cxsum_vv_sq_ip.SetSize(6, NTFCols);
}
}
// HV if needed
if (params.doHV_){
cxsum_hv.SetSize(16, vismtx.NCols());
cxsum_hv.SetSize(16, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
cxsum_hv_sq_rp.SetSize(16, NTFCols);
cxsum_hv_sq_ip.SetSize(16, NTFCols);
}
}
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
//allocating 8 Auto-Corr time-frequency maps
cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) );
......@@ -262,8 +288,7 @@ int main(int narg, const char* arg[])
for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
}
// HV if needed
if (params.doHV_){
//allocating 16 Cross-Corr H-V time-frequency maps
cout<<"visiavg/Info: allocating 16 H-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
......@@ -280,7 +305,8 @@ int main(int narg, const char* arg[])
datestart = TimeStamp(cfdate.DaysPart(),0.);
} // end if cnt==0
// HV if needed
if (I==0) { // start filling a new time bin
dateobs=cfdate;
if (prtlev>0)
......@@ -291,7 +317,8 @@ int main(int narg, const char* arg[])
acsum_sq = 0.;
cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.;
}
}
// VV if needed
if (params.doVV_){
cxsum_vv = complex<r_4>(0.,0.);
......@@ -300,6 +327,7 @@ int main(int narg, const char* arg[])
cxsum_vv_sq_ip = 0.;
}
}
// HV if needed
if (params.doHV_){
......@@ -309,34 +337,39 @@ int main(int narg, const char* arg[])
cxsum_hv_sq_ip = 0.;
}
}
}// end if(I==0)
// sum (integration) along the time axis
for(size_t k=0; k<KVAC.size(); k++) acsum.Row(k) += vismtx.Row(KVAC[k]); // Les auto-correlations
for(size_t k=0; k<KVAC.size(); k++){
acsum.Row(k) += (vismtx.Row(KVAC[k])).SubMatrix( Range::all(),Range(ifrmin,ifrmax));
}
//cout << " end ac "<<endl;
if (params.doSigma_) { // for computing Sigma, if required
for(size_t k=0; k<KVAC.size(); k++) {
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax))) ; //(Range(ifrmin,ifrmax)));
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
for(size_t k=0; k<KVCXHH.size(); k++){
cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
cxsum.Row(k) += vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)); // (Range(ifrmin,ifrmax)); // les cross-correlations
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));