Commit 13d6684d authored by perdereau's avatar perdereau

visiavg : bugS corection ; viewtfmao -> python3

parent f9936cf1
...@@ -7,6 +7,8 @@ from matplotlib import pyplot as plt ...@@ -7,6 +7,8 @@ from matplotlib import pyplot as plt
import itertools as itert import itertools as itert
from skimage import data, img_as_float from skimage import data, img_as_float
from skimage import exposure from skimage import exposure
from matplotlib.ticker import StrMethodFormatter
from matplotlib.ticker import FormatStrFormatter
import mynormalize as mynorm import mynormalize as mynorm
from matplotlib.widgets import Slider, Button, RadioButtons from matplotlib.widgets import Slider, Button, RadioButtons
import time import time
...@@ -15,7 +17,10 @@ from colorsys import hsv_to_rgb, hls_to_rgb ...@@ -15,7 +17,10 @@ from colorsys import hsv_to_rgb, hls_to_rgb
# values of the freqs for the slices # values of the freqs for the slices
rfreq = [1347.,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" normf=[1325.,1465.] # "normal frequencies"
names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H', names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H',
...@@ -83,7 +88,7 @@ matplotlib.rc('font', **font) ...@@ -83,7 +88,7 @@ matplotlib.rc('font', **font)
# #
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False, 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,submed=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) : ramod=False,plras=[],nmras=[],mask=False,logsc=False,snratio=False,usecm="",extaxes=False,choose_cmap=False,nau=True) :
if verb : if verb :
...@@ -166,6 +171,12 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F ...@@ -166,6 +171,12 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
print (np.shape(img)) print (np.shape(img))
for idx in np.arange(np.shape(img)[0]): for idx in np.arange(np.shape(img)[0]):
img[idx]=img[idx]-meds[idx] 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 mod :
if num>7 : if num>7 :
...@@ -378,14 +389,21 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F ...@@ -378,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 if raplot : # does plots of signal vs RA at some frequencies
fig,ax1=plt.subplots(figsize=winwid) 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.xlabel("RA (hour)")
plt.ylabel("Signal ") plt.ylabel("Signal ")
plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max()) plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
leg = plt.legend() leg = plt.legend()
for legobj in leg.legendHandles: for legobj in leg.legendHandles:
legobj.set_linewidth(2.0) 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.)) plt.xlim((0.,24.))
if mtitle =="": if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08) plt.title(filename + ' ' + names[num],y=1.08)
...@@ -402,20 +420,26 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F ...@@ -402,20 +420,26 @@ 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 if timeplot : # does plots of sig vs time at some frequencies
fig,ax1=plt.subplots(figsize=winwid) fig,ax1=plt.subplots(figsize=winwid)
[ plt.plot(times/3600., img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ] [ plt.plot(times/3600., img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
plt.xlabel('Time (TU, hour)') plt.xticks(fontsize=28)
plt.ylabel("Signal ") 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()) plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
#legend with thicker linewidth #legend with thicker linewidth
leg = plt.legend() leg = plt.legend(prop={'size': 28},loc=1 )
for legobj in leg.legendHandles: for legobj in leg.legendHandles:
legobj.set_linewidth(2.0) legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
plt.ylim((mvmin,mvmax))
plt.tight_layout()
if mtitle =="": if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08) plt.title(filename + ' ' + names[num],y=1.08)
else : else :
plt.title(mtitle+ " - " + names[num] ,y=1.08) plt.title(mtitle+ " - " + names[num] ,y=1.08,size=28)
plt.tight_layout() plt.tight_layout()
if save: if save:
if mask : if mask :
......
...@@ -195,7 +195,7 @@ int main(int narg, const char* arg[]) ...@@ -195,7 +195,7 @@ int main(int narg, const char* arg[])
double freqmin = P4FreqBand::freqstart_; double freqmin = P4FreqBand::freqstart_;
double freqmax = P4FreqBand::freqend_; double freqmax = P4FreqBand::freqend_;
sa_size_t ifrmin = 0 ; sa_size_t ifrmin = 0 ;
sa_size_t ifrmax = P4FreqBand::nbnufft_ ; sa_size_t ifrmax = P4FreqBand::nbnufft_ - 1 ; // ne pas oublier le -1 !
while (fgok) { while (fgok) {
//reads next visibility matrix window //reads next visibility matrix window
...@@ -212,7 +212,7 @@ int main(int narg, const char* arg[]) ...@@ -212,7 +212,7 @@ int main(int narg, const char* arg[])
vismtx = wreader.getAverageVisMtx(cfdate); vismtx = wreader.getAverageVisMtx(cfdate);
///cout << vismtx << endl; ///cout << vismtx << endl;
sa_size_t NTFCols = vismtx.NCols(); // bin en freq a ecrire 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 if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
tfmSX=wreader.getTotalNbWindows()/deltaIavg; tfmSX=wreader.getTotalNbWindows()/deltaIavg;
...@@ -288,8 +288,7 @@ int main(int narg, const char* arg[]) ...@@ -288,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) ); for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
} }
} }
// HV if needed
if (params.doHV_){ if (params.doHV_){
//allocating 16 Cross-Corr H-V time-frequency maps //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; cout<<"visiavg/Info: allocating 16 H-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
...@@ -306,7 +305,8 @@ int main(int narg, const char* arg[]) ...@@ -306,7 +305,8 @@ int main(int narg, const char* arg[])
datestart = TimeStamp(cfdate.DaysPart(),0.); datestart = TimeStamp(cfdate.DaysPart(),0.);
} // end if cnt==0 } // end if cnt==0
// HV if needed
if (I==0) { // start filling a new time bin if (I==0) { // start filling a new time bin
dateobs=cfdate; dateobs=cfdate;
if (prtlev>0) if (prtlev>0)
...@@ -317,7 +317,8 @@ int main(int narg, const char* arg[]) ...@@ -317,7 +317,8 @@ int main(int narg, const char* arg[])
acsum_sq = 0.; acsum_sq = 0.;
cxsum_sq_rp = 0.; cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.; cxsum_sq_ip = 0.;
} }
// VV if needed // VV if needed
if (params.doVV_){ if (params.doVV_){
cxsum_vv = complex<r_4>(0.,0.); cxsum_vv = complex<r_4>(0.,0.);
...@@ -326,6 +327,7 @@ int main(int narg, const char* arg[]) ...@@ -326,6 +327,7 @@ int main(int narg, const char* arg[])
cxsum_vv_sq_ip = 0.; cxsum_vv_sq_ip = 0.;
} }
} }
// HV if needed // HV if needed
if (params.doHV_){ if (params.doHV_){
...@@ -335,13 +337,11 @@ int main(int narg, const char* arg[]) ...@@ -335,13 +337,11 @@ int main(int narg, const char* arg[])
cxsum_hv_sq_ip = 0.; cxsum_hv_sq_ip = 0.;
} }
} }
}// end if(I==0) }// end if(I==0)
// sum (integration) along the time axis // sum (integration) along the time axis
for(size_t k=0; k<KVAC.size(); k++){ for(size_t k=0; k<KVAC.size(); k++){
//cout << vismtx.Row(KVAC[k]) << endl;
//cout << (vismtx.Row(KVAC[k])).SubMatrix(Range::all(),Range(ifrmin,ifrmax))<< endl;
//cout <<acsum.Row(k) << endl;
acsum.Row(k) += (vismtx.Row(KVAC[k])).SubMatrix( Range::all(),Range(ifrmin,ifrmax)); acsum.Row(k) += (vismtx.Row(KVAC[k])).SubMatrix( Range::all(),Range(ifrmin,ifrmax));
} }
...@@ -352,7 +352,7 @@ int main(int narg, const char* arg[]) ...@@ -352,7 +352,7 @@ int main(int narg, const char* arg[])
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ; acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
} }
} }
//cout << " end ac "<<endl;
for(size_t k=0; k<KVCXHH.size(); k++){ for(size_t k=0; k<KVCXHH.size(); k++){
cxsum.Row(k) += vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)); // (Range(ifrmin,ifrmax)); // 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 if (params.doSigma_) { // for computing Sigma, if required
...@@ -704,20 +704,23 @@ int main(int narg, const char* arg[]) ...@@ -704,20 +704,23 @@ int main(int narg, const char* arg[])
// --- FIN sauvegarde cartes temps-frequence // --- FIN sauvegarde cartes temps-frequence
// resu.Update(); // resu.Update();
double frbase = myp4fre.getP4Frequency(0);
size_t szFreqs = myp4fre.getP4NbFreqChannels()/TFMfbin; size_t szFreqs = myp4fre.getP4NbFreqChannels()/TFMfbin;
if (params.fbands_.size()>0 ) cout << szFreqs << " taille des freqs " << endl;
if (params.fbands_.size()>0 ){
szFreqs = (params.fbands_[0].jfmax_ - params.fbands_[0].jfmin_)/TFMfbin; szFreqs = (params.fbands_[0].jfmax_ - params.fbands_[0].jfmin_)/TFMfbin;
frbase = myp4fre.getP4Frequency(params.fbands_[0].jfmin_) + myp4fre.getP4FreqResolution()/2. ;
}
TVector <double> avg_freqs( szFreqs ); TVector <double> avg_freqs( szFreqs );
double frbase = myp4fre.getP4Frequency(params.fbands_[0].jfmin_) + myp4fre.getP4FreqResolution()/2. ;
cout << " base freq "<<frbase <<endl;
for (int kf=0 ; kf< szFreqs ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin ) for (int kf=0 ; kf< szFreqs ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
avg_freqs(kf) = frbase ; avg_freqs(kf) = frbase ;
potfm << PPFNameTag("FreqVec") << avg_freqs; potfm << PPFNameTag("FreqVec") << avg_freqs;
cout << " freq ecrites"<<endl;
if (fos != NULL) { if (fos != NULL) {
ext_names.push_back("FreqsLims"); ext_names.push_back("FreqsLims");
(*fos)<< lim_freq ; (*fos)<< lim_freq ;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment