Skip to content
Snippets Groups Projects
Commit c8e636c9 authored by Vlad Kondratiev's avatar Vlad Kondratiev
Browse files

lofar2-psr-quickview.py is updated to work also for Stokes I data; the script... 

lofar2-psr-quickview.py is updated to work also for Stokes I data; the script was also removed to lofar2-psr-quickview_l2com.py; get_rfi_fraction.py was also renamed to get_rfi_fraction_l2com.py
parent 0e129342
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
l2com/get_rfi_fraction.py → l2com/get_rfi_fraction_l2com.py
+ 0
− 0
View file @ c8e636c9
File moved
l2com/lofar2-psr-quickview.py → l2com/lofar2-psr-quickview_l2com.py
+ 102
− 57
View file @ c8e636c9
......@@ -254,12 +254,28 @@ if __name__ == "__main__":
elev= str(os.popen(cmd).readlines()[0][:-1])
cmd="azza.pl -ra %s -dec %s -t %s | grep AZ | cut -d ' ' -f 3" % (ra, dec, timangles[key][0])
azim= str(os.popen(cmd).readlines()[0][:-1])
cmd="azza.pl -ra %s -dec %s -t %s | grep HA | cut -d ' ' -f 3-4" % (ra, dec, timangles[key][0])
cmd="azza.pl -ra %s -dec %s -t %s | grep HA | cut -d ' ' -f 3" % (ra, dec, timangles[key][0])
ha= str(os.popen(cmd).readlines()[0][:-1])
timangles[key].extend([elev, azim, ha])
cmd="azza.pl -ra %s -dec %s -t %s | grep HA | cut -d ' ' -f 4" % (ra, dec, timangles[key][0])
haunit= str(os.popen(cmd).readlines()[0][:-1])
# change seconds in haunit to minutes
if haunit == "sec":
ha = ".2f" % (float(ha) / 60.)
haunit = "min"
if haunit == "h": haunit = "hr"
timangles[key].extend([elev, azim, ha, haunit])
# at this point HA units can either be in minutes or hours. So, if they are not the same, then we will convert those in minutes to hours
if timangles["startobs"][4] == timangles["midobs"][4] and timangles["midobs"][4] == timangles["endobs"][4]:
haunit = timangles["midobs"][4] # HA units are the same at the start, mid and end of observation
else:
haunit = "hr"
for key in timangles.keys():
if timangles[key][4] == "min":
timangles[key][3] = ".2f" % (float(timangles[key][3]) / 60.)
timangles[key][4] = "hr"
datainfo['elevation'] = ['Elevation (in deg) at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][1], timangles["midobs"][1], timangles["endobs"][1])]
datainfo['azimuth'] = ['Azimuth (in deg) at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][2], timangles["midobs"][2], timangles["endobs"][2])]
datainfo['transit'] = ['Hour angle at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][3], timangles["midobs"][3], timangles["endobs"][3])]
datainfo['transit'] = ['Hour angle (in %s) at the start / mid-point / end' % (haunit), "%s / %s / %s" % (timangles["startobs"][3], timangles["midobs"][3], timangles["endobs"][3])]
if args.paz_rfi_fraction != "undef" or args.clfd_rfi_fraction != "undef":
datainfo['rfi'] = ['RFI fraction (in %) by paz / clfd', "%s / %s" % (args.paz_rfi_fraction, args.clfd_rfi_fraction)]
......@@ -267,6 +283,7 @@ if __name__ == "__main__":
workraw = pc.Archive_load(infile)
if not(workraw.get_dedispersed()):
workraw.dedisperse()
npol = workraw.get_npol()
workraw.fscrunch(args.fscrunch)
workraw.tscrunch(args.tscrunch)
bandpass=workraw.clone()
......@@ -287,6 +304,7 @@ if __name__ == "__main__":
stokesraw.bscrunch(args.bscrunch)
nbins = stokesraw.get_nbin()
if npol > 1:
stokesraw.convert_state('Stokes')
r = stokesraw.get_data()
#time stokes f phase
......@@ -296,8 +314,11 @@ if __name__ == "__main__":
fig = plt.figure(figsize=(8.27, 11.69), dpi=100, constrained_layout=True)
gs = fig.add_gridspec(5, 1, hspace=0.01, height_ratios=[2,1,1,1, 0.1])
gs0 = gs[0].subgridspec(1,2, wspace=0.02, hspace=0.01, width_ratios=[1, 2])
if npol > 1:
gs1 = gs[1].subgridspec(1,4, wspace=0.02, hspace=0.01, width_ratios=[1.5, 1, 1, 1])
gs2 = gs[2].subgridspec(1,2, wspace=0.13, hspace=0.01)
else:
gs1 = gs[1].subgridspec(1,2, wspace=0.15, hspace=0.01)
gs2 = gs[2].subgridspec(1,2, wspace=0.15, hspace=0.01)
gs3 = gs[3].subgridspec(1,2, wspace=0.02, hspace=0.01)
gs4 = gs[4].subgridspec(1,1, wspace=0.0, hspace=0.0)
......@@ -321,7 +342,7 @@ if __name__ == "__main__":
range_binpeak = np.argmax(range_prof)
range_snrmean = np.mean(range_prof)
range_weq = abs(np.sum(range_prof)/range_snrpeak)
binshift = int(5*range_weq)
binshift = int(10*range_weq)
if binshift > nbins//2: binshift=nbins//2
lbin=nbins//2-binshift
rbin=nbins//2+binshift
......@@ -362,6 +383,7 @@ if __name__ == "__main__":
ax11 = fig.add_subplot(gs0[0, 0]) # profile
ax11.plot(phases, (data1[0,lbin:rbin]-mean)/rms, '-', color='black', lw=1, label='I')
if npol > 1:
linpol=np.array([np.sqrt(data1[1][bin]*data1[1][bin]+data1[2][bin]*data1[2][bin]) for bin in range(nbins)])
ax11.plot(phases, (linpol[lbin:rbin]-mean)/rms, '-', color='red', lw=1, label='L')
ax11.plot(phases, (data1[3,lbin:rbin]-mean)/rms, '-', color='blue', lw=1, label='V')
......@@ -377,6 +399,7 @@ if __name__ == "__main__":
if nsubint > 1: stokes_iquv.tscrunch()
stokes_iquv.bscrunch(args.bscrunch)
nbins = stokes_iquv.get_nbin()
if npol > 1:
stokes_iquv.convert_state('Stokes')
r = stokes_iquv.get_data()
#time stokes f phase
......@@ -391,6 +414,7 @@ if __name__ == "__main__":
plt.xlabel("Pulse phase", labelpad=0)
plt.ylabel("Frequency (MHz)")
if npol > 1:
ax22 = fig.add_subplot(gs1[0, 1]) # stokes Q
sp=normalize(data2[1,:,:], nchan)
v=sp[:,lbin:rbin]
......@@ -432,6 +456,23 @@ if __name__ == "__main__":
ax24.yaxis.set_minor_formatter(FormatStrFormatter(""))
ax24.annotate(" V ", xycoords="axes fraction", xy=(0.95, 0.93), bbox=dict(facecolor='0.95', boxstyle="round"), ha='right', va='top', multialignment="right")
plt.xlabel("Pulse phase", labelpad=0)
else:
ax22 = fig.add_subplot(gs1[0, 1]) # samples histogram
from scipy.optimize import leastsq
fitfunc = lambda p, x: p[0]*np.exp(-0.5*((x-mean)/rms)**2)+p[1]
errfunc = lambda p, x, y: (y - fitfunc(p, x))
init = [1.0, 0.5]
histbins=500
n, bins, patches = ax22.hist(data1[0], bins=histbins, color='skyblue', alpha=0.6)
out = leastsq(errfunc, init, args=(bins[:histbins], n[:histbins]))
ax22.plot(bins, fitfunc(out[0], bins), color='green', linewidth=0.5)
ax22.set_xlim(xmax=30*rms)
ax22.set_ylim(ymin=-1)
ax22.annotate("$\mu$ = %.3g\n$\sigma$ = %.3g" % (mean, rms), xycoords="axes fraction", xy=(0.97, 0.9), bbox=dict(facecolor='0.95', boxstyle="round"), ha='right', va='top', multialignment="right", fontsize=8)
ax22.grid(which='both', ls=':')
plt.xlabel("Sample value", labelpad=0)
plt.ylabel("Number of samples", labelpad=0)
ax31 = fig.add_subplot(gs2[0, 0]) # time-phase
timephase=workraw.clone()
......@@ -455,11 +496,15 @@ if __name__ == "__main__":
data4 = np.mean(np.sum(r[:,:,:,:], axis=3), axis=0)
dB = lambda x: 10 * np.log10(x)
freqs=np.array([lowfreq + chan * chan_bw for chan in range(nchan)])
if npol > 1:
ax32.plot(freqs, dB(data4[0,:]), color="k", label="PP", lw=0.5)
ax32.plot(freqs, dB(data4[1,:]), color="b", label="QQ", lw=0.5)
#ax32.plot(freqs, dB(data4[2,:]), color="g", label="Re[PQ]")
#ax32.plot(freqs, dB(data4[3,:]), color="r", label="Im[PQ]")
ax32.legend(ncol=2)
else:
ax32.plot(freqs, dB(data4[0,:]), color="k", label="Stokes I", lw=1)
ax32.legend()
ax32.grid(which='both', ls=':')
plt.xlabel("Frequency (MHz)", labelpad=0)
plt.ylabel("Power (dB)")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment