added new script to l2com/ for commissioning

l2com/get_rfi_fraction.py

+#!/usr/bin/env python
+#
+# Reads the input PSRFITS file to calculate the RFI fraction
+# based on the samples' weights
+# 
+# Vlad Kondratiev (c) - 11.04.2018
+#
+import os, sys
+import numpy as np
+import optparse as opt
+import psrchive as pc
+
+# Main
+if __name__=="__main__":
+        #
+        # Parsing the command line options
+        #
+        usage = "Usage: %prog [options] <ar-file>"
+        cmdline = opt.OptionParser(usage)
+        cmdline.add_option('-v', '--verbose', dest='is_verbose', action="store_true", help="Verbose output, print extra info", default=False)
+
+        # reading cmd options
+        (opts,args) = cmdline.parse_args()
+
+        # check if input file is given
+        if len(args) == 0:
+                cmdline.print_usage()
+                sys.exit(0)
+
+        # input file
+        infile = args[0]
+
+        # reading input ar-file
+        try:
+                raw = pc.Archive_load(infile)
+        except:
+                print ("Can't open input file: '%s'" % (infile))
+                sys.exit(1)
+        if not raw.get_dedispersed():
+                raw.dedisperse()
+        # raw.remove_baseline() - we subtract it outselves
+        raw.pscrunch()
+
+        orig_nsubint = raw.get_nsubint()
+        orig_nchan = raw.get_nchan()
+        if opts.is_verbose: print ("nsubint = %d   nchan = %d" % (orig_nsubint, orig_nchan))
+
+        # checking the weights and re-normalizing them if needed
+        weights=raw.get_weights()
+        subint_durs = [raw.get_Integration(ii).get_duration() for ii in range(orig_nsubint)]
+        max_weight = np.max(weights)
+        if opts.is_verbose: print ("max weight = %f" % (max_weight))
+        max_subint_dur = np.max(subint_durs)
+        badpoints = 0
+        for ii in range(orig_nsubint):
+                isub=raw.get_Integration(ii)
+                dur=isub.get_duration()
+                for ch in range(orig_nchan):
+                        w=isub.get_weight(ch)
+                        wold = w
+                        if w<max_weight:
+                                w *= (max_subint_dur/dur)
+                                if w>max_weight: w=max_weight
+                        w /= max_weight
+                        isub.set_weight(ch, float(w))
+                        if w < max_weight: badpoints += 1
+                        #if opts.is_verbose: print ("sub=%d  chan=%d  wold=%f   w=%f" % (ii, ch, wold, w))
+
+        # re-reading weights again and calculate the RFI fraction
+        weights=raw.get_weights()
+        rfi_fraction=1.-np.sum(weights)/(orig_nsubint*orig_nchan)
+        if opts.is_verbose: print ("RFI fraction = %f   Bad subints/chans = %d (%d)" % (rfi_fraction, badpoints, orig_nsubint * orig_nchan))
+        else: print (rfi_fraction)

l2com/h5modify.py

+#!/usr/bin/env python3
+#
+import os, sys
+import numpy as np
+import matplotlib.pyplot as plt
+import argparse
+
+import h5py
+import psutil
+
+import warnings
+warnings.simplefilter('ignore')
+
+if len(sys.argv) == 1:
+        print ("No input files!")
+        sys.exit(0)
+
+fname = sys.argv[1] # h5 file
+
+h5 = h5py.File(fname, "r+")
+sap_ids = [key for key in h5.keys() if "SUB_ARRAY_POINTING" in key]
+for sap_id in sap_ids:
+    beam_ids = [key for key in h5[sap_id].keys() if "BEAM" in key]
+    for beam_id in beam_ids:
+        group_name = f"{sap_id}/{beam_id}"
+        h5[group_name].attrs["OBSERVATION_NOF_STOKES"] = 1
+        print ("OBSERVATION_NOF_STOKES = ", h5[group_name].attrs["OBSERVATION_NOF_STOKES"])
+        stokes_id = [key for key in h5[group_name].keys() if "STOKES_" in key][0]
+        print ("STOKES GROUP: ", stokes_id)
+        if stokes_id not in 'STOKES_0':
+            h5[group_name]['STOKES_0'] = h5[group_name][stokes_id]
+            print ("%s copied to STOKES_0" % (stokes_id))
+        break
+
+#print (h5[group_name].attrs["OBSERVATION_NOF_STOKES"])

l2com/lofar2jones.py

+#!/usr/bin/env python
+
+# Copyright (C) 2015 by Maciej Serylak
+# Licensed under the Academic Free License version 3.0
+# This program comes with ABSOLUTELY NO WARRANTY.
+# You are free to modify and redistribute this code as long
+# as you do not remove the above attribution and reasonably
+# inform receipients that you have modified the original work.
+
+import os
+import optparse as opt
+import sys
+import numpy as np
+
+from datetime import datetime, timedelta
+from pyrap.measures import measures
+from pyrap.quanta import quantity
+from antennaJones import getJonesByAntFld
+
+
+# Main
+if __name__=="__main__":
+  #
+  # Parsing the command line options
+  #
+  usage = "Usage: %prog --parset <parset_file>"
+  cmdline = opt.OptionParser(usage)
+  cmdline.formatter.max_help_position = 50 # increase space reserved for option flags (default 24), trick to make the help more readable
+  cmdline.formatter.width = 200 # increase help width from 120 to 200
+  cmdline.add_option('--parset', type = 'str', dest = 'pfile', default = '', metavar = '*.parset', help = 'Parset file with observation details.')
+  cmdline.add_option('--subint', type = 'int', dest = 'subint', default = '5', help = 'Length of subintegration.')
+  cmdline.add_option('--output', type = 'str', dest = 'output', help = 'Output file.')
+
+  # reading cmd options
+  (opts, args) = cmdline.parse_args()
+
+  if not opts.pfile:
+    cmdline.print_usage()
+    sys.exit(0)
+
+  # Define constans and variables.
+  subbandNumber = 512   # Number of subbands (fixed at 512)
+  clockFrequency = 200.0   # Sample clock freq. MHz (typically 200)
+  nyquist_zone = 0
+  parsetFile = opts.pfile
+  subintLength = opts.subint
+  outputFile = opts.output
+  if not outputFile:
+    outputFile = parsetFile.split('.')[0] + ".jones"
+
+  with open(parsetFile, 'r') as searchFile:
+      for line in searchFile:
+          if 'Observation.bandFilter' in line:
+              bandFilter = line.split('=')[1]
+              if 'LBA_10_90' in bandFilter or 'LBA_30_90' in bandFilter:
+                 nyquist_zone = 0
+              if 'HBA_110_190' in bandFilter:
+                 nyquist_zone = 1
+              if 'HBA_210_250' in bandFilter:
+                 nyquist_zone = 2
+          if 'Observation.startTime' in line:
+              startTime = line.split('=')[1].strip(' []\'\n')
+              bTime = datetime.strptime(startTime, "%Y-%m-%d %H:%M:%S")
+              stepTime = timedelta(0, float(subintLength))
+          if 'Observation.stopTime' in line:
+              stopTime = line.split('=')[1].strip(' []\'\n')
+              eTime = datetime.strptime(stopTime, "%Y-%m-%d %H:%M:%S")
+              taskDuration=(eTime - bTime).total_seconds()
+          if 'Observation.VirtualInstrument.stationList' in line:
+          #if 'Cobalt.BeamFormer.stationList' in line:
+              stationList = line.split('=')[1].strip(' []\'\n').split(',')
+              stationList.sort()
+              stationList = [item.replace('HBA0', '') for item in stationList]
+              stationList = [item.replace('HBA1', '') for item in stationList]
+              #newStationList = []
+              #for item in stationList:
+              #    if item not in newStationList:
+              #        newStationList.append(item)
+          if 'Observation.Beam[0].TiedArrayBeam[0].angle1' in line:
+              targetRA = line.split('=')[1].strip(' []\'\n')
+          if 'Observation.Beam[0].TiedArrayBeam[0].angle2' in line:
+              targetDec = line.split('=')[1].strip(' []\'\n')
+          if 'Observation.Beam[0].subbandList' in line:
+              subbandList = line.split('=')[1].strip(' []\'\n').replace(' ', '')
+              subbandList = sorted([int(ss) for ss in subbandList.split(',')])
+              subbandStart = subbandList[0]
+              subbandEnd = subbandList[-1]
+
+  direction=measures().direction('J2000', str(float(targetRA,))+'rad', str(float(targetDec,))+'rad')
+
+  Times=[]
+  for sec in range(int(float(taskDuration)/subintLength)):
+      Times.append( (quantity((bTime + sec * stepTime).isoformat())).get_value() )
+
+  obsTimes=quantity(Times,'d')
+  obsTimesArr=obsTimes.get_value()
+  jones=[] # array of jones matrices in freq
+
+  newsubStart = 0
+  for subband in range(subbandStart, subbandEnd + 1, 1):
+      subbandFrequency = (nyquist_zone*100 + (subband * clockFrequency / 2.0 / subbandNumber)) * 1000000
+      # we need to exclude all freqs below 10 MHz, they are not supported:
+      # Unsupported frequency value: 4.6875e+06
+      # frequency needs to be within 10e6 and 300 e6 Hz
+#      print (subband, subbandFrequency)
+      if subbandFrequency < 10.0e6 or subbandFrequency >=300e6:
+          continue
+      if newsubStart == 0: newsubStart = subband
+      print ("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%3d %12.1f Hz" % (subband, subbandFrequency), end='', flush=True)
+      Jn=getJonesByAntFld("Hamaker", obsTimes, stationList[0], direction, subbandFrequency)
+      jones.append(Jn)
+
+  print ()
+  out = open(outputFile, 'w')
+  for ii in range(len(Times)):
+      for jj in range(len(jones)):
+          subband = newsubStart + jj
+          subbandFrequency = (nyquist_zone*100 + (subband * clockFrequency / 2.0 / subbandNumber)) * 1000000
+          out.write("%f %f %f %f %f %f %f %f %f %f\n" % (obsTimesArr[ii], subbandFrequency, np.real(jones[jj][ii,0,0]), np.imag(jones[jj][ii,0,0]), \
+		np.real(jones[jj][ii,0,1]), np.imag(jones[jj][ii,0,1]), np.real(jones[jj][ii,1,0]), np.imag(jones[jj][ii,1,0]), np.real(jones[jj][ii,1,1]), np.imag(jones[jj][ii,1,1])))
+  out.close()

l2com/lofar_bf_h5info.py

 #!/usr/bin/env python3
 #
-# Prints out the keys and attributes of the root group (/) of HDF5 metadata
+# Prints out the keys and attributes of all groups of HDF5 metadata
 # of the LOFAR BF .h5 file
 #
 # Vlad Kondratiev, ASTRON (c) - 28.05.2025
@@ -13,16 +13,74 @@ import h5py
 if __name__ == "__main__":
     # Read command line arguments
     parser = argparse.ArgumentParser(prog='lofar_bf_h5info.py', usage='%(prog)s [options] <h5_file>', \
-        description="Prints out keys/attrs of the Root group of HDF5 metadata")
+        description="Prints out keys/attrs of all groups of HDF5 metadata")
     parser.add_argument('h5_file', nargs='?', help='input .h5 file')
+    parser.add_argument("--show-process-history", help="print out keys/attrs for the PROCESS_HISTORY group",
+                        action="store_true")
     args = parser.parse_args(args=None if sys.argv[1:] else ['--help'])
 
     try:
-        h5 = h5py.File(args.h5_file)
+        h5 = h5py.File(args.h5_file, 'r')
         group = h5["/"]
         keys = sorted(["%s" % item for item in sorted(list(group.attrs))])
         for key in keys:
             print(key + " = " + str(group.attrs[key]))
+
+        sap_ids = [key for key in h5.keys() if "SUB_ARRAY_POINTING" in key]
+        for sap_id in sap_ids:
+            print ("-" * 32)
+            print (sap_id)
+            group_name = f"{sap_id}"
+            skeys = sorted(["%s" % item for item in sorted(list(h5[group_name].attrs))])
+            for skey in skeys:
+                print(skey + " = " + str(h5[group_name].attrs[skey]))
+
+            beam_ids = [key for key in h5[group_name].keys() if "BEAM" in key]
+            for beam_id in beam_ids:
+                print ("-" * 32)
+                print (beam_id)
+                group_name = f"{sap_id}/{beam_id}"
+                bkeys = sorted(["%s" % item for item in sorted(list(h5[group_name].attrs))])
+                for bkey in bkeys:
+                    print(bkey + " = " + str(h5[group_name].attrs[bkey]))
+
+                stokes_ids = [key for key in h5[group_name].keys() if "STOKES_" in key]
+                for stokes_id in stokes_ids:
+                    print ("-" * 32)
+                    print (stokes_id)
+                    stokes_name = f"{sap_id}/{beam_id}/{stokes_id}" 
+                    stkeys = sorted(["%s" % item for item in sorted(list(h5[stokes_name].attrs))])
+                    for stkey in stkeys:
+                        print(stkey + " = " + str(h5[stokes_name].attrs[stkey]))
+
+                coord_ids = [key for key in h5[group_name].keys() if "COORDINATES" in key]
+                for coord_id in coord_ids:
+                    print ("-" * 32)
+                    print (coord_id)
+                    coord_name = f"{sap_id}/{beam_id}/{coord_id}" 
+                    ckeys = sorted(["%s" % item for item in sorted(list(h5[coord_name].attrs))])
+                    for ckey in ckeys:
+                        print(ckey + " = " + str(h5[coord_name].attrs[ckey]))
+
+                    time_freq_ids = [key for key in h5[coord_name].keys()]
+                    for tf_id in time_freq_ids:
+                        print ("-" * 32)
+                        print (tf_id)
+                        time_freq_name = f"{sap_id}/{beam_id}/{coord_id}/{tf_id}"
+                        tfkeys = sorted(["%s" % item for item in sorted(list(h5[time_freq_name].attrs))])
+                        for tfkey in tfkeys:
+                            print(tfkey + " = " + str(h5[time_freq_name].attrs[tfkey]))
+
+                if (args.show_process_history):
+                    hist_ids = [key for key in h5[group_name].keys() if "HISTORY" in key]
+                    for hist_id in hist_ids:
+                        print ("-" * 32)
+                        print (hist_id)
+                        hist_name = f"{sap_id}/{beam_id}/{hist_id}" 
+                        hkeys = sorted(["%s" % item for item in sorted(list(h5[hist_name].attrs))])
+                        for hkey in hkeys:
+                            print(hkey + " = " + str(h5[hist_name].attrs[hkey]))
+
     except FileNotFoundError as err:
         print ("File not found:\n%s" % (err))
         sys.exit(1)

l2com/pol-xxyy-stokes-uncal.py

+#!/usr/bin/env python
+#
+import os, sys
+import numpy as np
+import psrchive as pc
+import scipy.stats as sc
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as colors
+from matplotlib.ticker import *
+import argparse
+
+tick_label_size = 10
+axes_label_size = 10
+#print mpl.rcParams.keys()
+mpl.rcParams['xtick.labelsize'] = tick_label_size
+mpl.rcParams['ytick.labelsize'] = tick_label_size
+mpl.rcParams['axes.labelsize'] = axes_label_size
+
+# normalizing each channel
+def normalize(x, nchan):
+    for ii in range(nchan):
+        osm, osr = sc.probplot(x[ii], sparams=(), dist='norm', fit=0)
+        q_max = np.min(np.where(osm > 1.0))
+        q_min = np.max(np.where(osm < -1.0))
+        rms, mean = np.polyfit(osm[q_min:q_max], osr[q_min:q_max], 1)
+        x[ii] -= mean
+        if rms == 0.0: x[ii] = 0.0
+        else: x[ii] /= rms
+        crit=np.isfinite(x[ii])
+        x[ii][np.logical_not(crit)] = 0.0
+    return x
+
+if __name__ == "__main__":
+    # Read command line arguments
+    parser = argparse.ArgumentParser(description="Plots Stokes components of input .ar file")
+    parser.add_argument("arfile", type=str, help="input .ar file")
+    parser.add_argument("-f", "--fscrunch", help="scrunch factor in frequency (default: 1)", default=1, type=int)
+    parser.add_argument("--psr", type=str, help="pulsar name for the plot's title")
+    parser.add_argument("--noshow", help="Only save dynamic spectrum to file without opening GUI", action="store_true")
+    args = parser.parse_args()
+
+    note="Frequency vs. Phase"
+    binshift=20
+    stokes=["I", "Q", "U", "V"]
+
+    # Read data
+    raw = pc.Archive_load(args.arfile)
+    if not(raw.get_dedispersed()):
+        raw.dedisperse()
+    raw.remove_baseline()
+    raw.centre_max_bin()
+    nchan = raw.get_nchan()
+    nsubint = raw.get_nsubint()
+    if args.fscrunch > 1: raw.fscrunch(args.fscrunch)
+    if nsubint > 1: raw.tscrunch()
+    nbins = raw.get_nbin()
+    r = raw.get_data()
+    #time stokes f phase
+    data = r[0,:,:,:]
+
+    fig = plt.figure(figsize=(8,6), dpi=100)
+    fig.suptitle("%s ,  %s" % (args.psr, note))
+    for pol in range(4):
+        ax = fig.add_subplot(2,2,pol+1)
+        lbin=nbins//2-binshift
+        rbin=nbins//2+binshift
+        sp=normalize(data[pol,:,lbin:rbin], nchan)
+        gmax = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+        ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax, vmin=-gmax)
+        ax.set_xticklabels([])
+        ax.set_yticklabels([])
+        #plt.xlabel("Phase")
+        plt.title(stokes[pol])
+
+    plt.tight_layout()
+    plt.subplots_adjust(top=0.9)
+    plt.gca().axes.get_yaxis().set_visible(False)
+
+    outpng=args.arfile.split(".ar")[0] + ".iquv.png"
+    plt.savefig(outpng, bbox_inches="tight")
+    if not args.noshow:
+        plt.show()
+    plt.close()

l2com/polcomp3n-iquv.py

+#!/usr/bin/env python
+#
+import os, sys
+import numpy as np
+import psrchive as pc
+import scipy.stats as sc
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as colors
+from matplotlib.ticker import *
+
+
+tick_label_size = 10
+axes_label_size = 10
+#print mpl.rcParams.keys()
+mpl.rcParams['xtick.labelsize'] = tick_label_size
+mpl.rcParams['ytick.labelsize'] = tick_label_size
+mpl.rcParams['axes.labelsize'] = axes_label_size
+
+
+psr="B1919+21"
+note="Frequency vs. Phase"
+#note2="(Raw, uncalibrated)"
+#note2="(Calibrated)"
+note2=''
+binshift=50
+
+# normalizing each channel
+def normalize(x, nchan):
+	for ii in range(nchan):
+		osm, osr = sc.probplot(x[ii], sparams=(), dist='norm', fit=0)
+		q_max = np.min(np.where(osm > 1.0))
+		q_min = np.max(np.where(osm < -1.0))
+		rms, mean = np.polyfit(osm[q_min:q_max], osr[q_min:q_max], 1)
+		x[ii] -= mean
+		if rms == 0.0: x[ii] = 0.0
+		else: x[ii] /= rms
+		crit=np.isfinite(x[ii])
+		x[ii][np.logical_not(crit)] = 0.0
+	return x
+
+
+if len(sys.argv) == 1:
+	print ("No input files!")
+	sys.exit(0)
+infile1=sys.argv[1] # S0
+infile2=sys.argv[2] # S1
+infile3=sys.argv[3] # S2
+infile4=sys.argv[4] # S3
+
+stokes=["I", "Q", "U", "V"]
+
+# reading 1st file
+raw = pc.Archive_load(infile1)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan1 = raw.get_nchan()
+nsubint1 = raw.get_nsubint()
+if nsubint1 > 1: raw.tscrunch()
+nbins1 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data1 = r[0,:,:,:]
+
+# reading 2nd file
+raw = pc.Archive_load(infile2)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan2 = raw.get_nchan()
+nsubint2 = raw.get_nsubint()
+if nsubint2 > 1: raw.tscrunch()
+nbins2 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data2 = r[0,:,:,:]
+
+# reading 3rd file
+raw = pc.Archive_load(infile3)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan3 = raw.get_nchan()
+nsubint3 = raw.get_nsubint()
+if nsubint3 > 1: raw.tscrunch()
+nbins3 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data3 = r[0,:,:,:]
+
+# reading 4th file
+raw = pc.Archive_load(infile4)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan4 = raw.get_nchan()
+nsubint4 = raw.get_nsubint()
+if nsubint4 > 1: raw.tscrunch()
+nbins4 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data4 = r[0,:,:,:]
+
+fig = plt.figure(figsize=(8,6), dpi=100)
+fig.suptitle("%s ,  %s %s" % (psr, note, note2))
+
+pol=0
+ax = fig.add_subplot(1,4,1)
+lbin=nbins1//2-binshift
+rbin=nbins1//2+binshift
+sp=normalize(data1[0,:,lbin:rbin], nchan1)
+gmax = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax, vmin=-gmax)
+ax.set_xticklabels([])
+ax.set_yticklabels([])
+#if pol==0: plt.ylabel("Uncalibrated")
+plt.title(stokes[pol])
+
+pol=1
+ax = fig.add_subplot(1,4,2)
+lbin=nbins2//2-binshift
+rbin=nbins2//2+binshift
+sp=normalize(data2[0,:,lbin:rbin], nchan2)
+gmax2 = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax2, vmin=-gmax2)
+ax.set_xticklabels([])
+ax.set_yticklabels([])
+#if pol==0: plt.ylabel("Calibrated")
+plt.title(stokes[pol])
+
+pol=2
+ax = fig.add_subplot(1,4,3)
+lbin=nbins3//2-binshift
+rbin=nbins3//2+binshift
+sp=normalize(data3[0,:,lbin:rbin], nchan3)
+gmax3 = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax3, vmin=-gmax3)
+ax.set_xticklabels([])
+ax.set_yticklabels([])
+#if pol==0: plt.ylabel("Calibrated")
+plt.title(stokes[pol])
+
+pol=3
+ax = fig.add_subplot(1,4,4)
+lbin=nbins4//2-binshift
+rbin=nbins4//2+binshift
+sp=normalize(data4[0,:,lbin:rbin], nchan4)
+gmax4 = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax4, vmin=-gmax4)
+ax.set_xticklabels([])
+ax.set_yticklabels([])
+#if pol==0: plt.ylabel("Calibrated")
+plt.title(stokes[pol])
+
+plt.tight_layout()
+plt.subplots_adjust(top=0.9)
+plt.gca().axes.get_yaxis().set_visible(False)
+plt.show()

l2com/polcomp3n.py

+#!/usr/bin/env python
+#
+import os, sys
+import numpy as np
+import psrchive as pc
+import scipy.stats as sc
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as colors
+from matplotlib.ticker import *
+
+
+tick_label_size = 10
+axes_label_size = 10
+#print mpl.rcParams.keys()
+mpl.rcParams['xtick.labelsize'] = tick_label_size
+mpl.rcParams['ytick.labelsize'] = tick_label_size
+mpl.rcParams['axes.labelsize'] = axes_label_size
+
+
+psr="B1919+21"
+note="Frequency vs. Phase"
+#note2="(Raw, uncalibrated)"
+#note2="(Calibrated)"
+note2=''
+binshift=20
+
+# normalizing each channel
+def normalize(x, nchan):
+	for ii in range(nchan):
+		osm, osr = sc.probplot(x[ii], sparams=(), dist='norm', fit=0)
+		q_max = np.min(np.where(osm > 1.0))
+		q_min = np.max(np.where(osm < -1.0))
+		rms, mean = np.polyfit(osm[q_min:q_max], osr[q_min:q_max], 1)
+		x[ii] -= mean
+		if rms == 0.0: x[ii] = 0.0
+		else: x[ii] /= rms
+		crit=np.isfinite(x[ii])
+		x[ii][np.logical_not(crit)] = 0.0
+	return x
+
+
+if len(sys.argv) == 1:
+	print ("No input files!")
+	sys.exit(0)
+infile1=sys.argv[1] # uncal
+infile2=sys.argv[2] # polcal
+
+stokes=["I", "Q", "U", "V"]
+
+# reading 1st file
+raw = pc.Archive_load(infile1)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan1 = raw.get_nchan()
+nsubint1 = raw.get_nsubint()
+#if nchan1 > 1: raw.fscrunch()
+if nsubint1 > 1: raw.tscrunch()
+nbins1 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data1 = r[0,:,:,:]
+
+# reading 2nd file
+raw = pc.Archive_load(infile2)
+if not(raw.get_dedispersed()):
+	raw.dedisperse()
+raw.remove_baseline()
+raw.centre_max_bin()
+nchan2 = raw.get_nchan()
+nsubint2 = raw.get_nsubint()
+#if nchan2 > 1: raw.fscrunch()
+if nsubint2 > 1: raw.tscrunch()
+nbins2 = raw.get_nbin()
+r = raw.get_data()
+#time stokes f phase
+data2 = r[0,:,:,:]
+
+fig = plt.figure(figsize=(8,6), dpi=100)
+fig.suptitle("%s ,  %s %s" % (psr, note, note2))
+for pol in range(4):
+	ax = fig.add_subplot(2,4,pol+1)
+	lbin=nbins1//2-binshift
+	rbin=nbins1//2+binshift
+	sp=normalize(data1[pol,:,lbin:rbin], nchan1)
+	#gmin=np.min(sp)
+	#gmax=np.max(sp)
+	#norm=colors.normalize(gmin, gmax)
+	gmax = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+	ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax, vmin=-gmax)
+	ax.set_xticklabels([])
+	ax.set_yticklabels([])
+	#plt.xlabel("Phase")
+	if pol==0: plt.ylabel("Uncalibrated")
+	plt.title(stokes[pol])
+
+	ax = fig.add_subplot(2,4,pol+5)
+	lbin=nbins2//2-binshift
+	rbin=nbins2//2+binshift
+	sp=normalize(data2[pol,:,lbin:rbin], nchan2)
+	gmax2 = np.max([np.abs(np.min(sp)), np.abs(np.max(sp))])
+	ax.imshow(sp, interpolation='none', aspect='auto', origin='lower', cmap="bwr", vmax=gmax2, vmin=-gmax2)
+	ax.set_xticklabels([])
+	ax.set_yticklabels([])
+	#plt.xlabel("Phase")
+	if pol==0: plt.ylabel("Calibrated")
+
+plt.tight_layout()
+plt.subplots_adjust(top=0.9)
+plt.gca().axes.get_yaxis().set_visible(False)
+plt.show()

l2com/process-iquv.sh

+#!/bin/bash
+#
+obsid=L2050325
+project=l2cb_022
+psr="1919+21"
+psr2=1919
+fscr=256
+parfile="../${psr}.par"
+opts="-t 8 -A -L 10 -b 512 -E $parfile"
+inputdir="cs"
+stem="b${psr2}-hba-low-iquv"
+
+for stokes in `seq -s ' ' 0 1 3`; do
+ stokes_dir=s$stokes
+ echo "STOKES = $stokes_dir"
+ mkdir $stokes_dir
+ echo "cp $inputdir/${obsid}_*_S${stokes}_*_bf.h5 ${stokes_dir}/"
+ cp $inputdir/${obsid}_*_S${stokes}_*_bf.h5 ${stokes_dir}/
+ cd ${stokes_dir}
+ rawfile=`ls -1 ../$inputdir/${obsid}_*_S${stokes}_*_bf.raw`
+ echo "ln -s $rawfile ."
+ ln -s $rawfile .
+ cd ..
+
+ echo "/home/kondratiev/src/h5modify.py ${stokes_dir}/*.h5"
+ /home/kondratiev/src/h5modify.py ${stokes_dir}/*.h5
+ cd ${stokes_dir}
+ h5name=`ls -1 *.h5`
+ h5new=`ls -1 *.h5 | sed -e s/_S${stokes}_/_S0_/`
+ echo "mv $h5name $h5new"
+ mv $h5name $h5new
+ cd ..
+ echo "dspsr $opts -O ${stem}-${stokes_dir} ${stokes_dir}/${obsid}_SAP000_B000_S0_P000_bf.h5"
+ dspsr $opts -O ${stem}-${stokes_dir} ${stokes_dir}/${obsid}_SAP000_B000_S0_P000_bf.h5
+
+ echo "clfd --no-report ${stem}-${stokes_dir}.ar"
+ clfd --no-report ${stem}-${stokes_dir}.ar ; mv ${stem}-${stokes_dir}.ar.clfd ${stem}-${stokes_dir}.clfd.ar
+ echo "pam -f ${fscr} -e f${fscr}.ar ${stem}-${stokes_dir}.clfd.ar"
+ pam -f ${fscr} -e f${fscr}.ar ${stem}-${stokes_dir}.clfd.ar
+done

l2com/process-xxyy.sh

+#!/bin/bash
+#
+obsid=L2050960
+project=l2cb_022
+psr="B1508+55"
+fscr=8
+nsubs=488
+subs_per_file=20
+parfile="../${psr:1}.par"
+antenna="lba" # hba
+ram="-U 32768 -t 1"
+opts="${ram} -A -L 10 -b 512 -E $parfile"
+inputdir="/data/projects/${project}/${obsid}/cs"
+stem="b${psr:1:4}-${antenna}-xxyy"
+
+nparts=$(($nsubs / $subs_per_file))
+rest=$(($nsubs - $nparts * $subs_per_file))
+if [[ $rest -eq 0 ]]; then
+ lastpart=$(($nparts - 1))
+else
+ lastpart=$nparts
+ nparts=$(($lastpart + 1))
+fi
+parts=`seq -s ' ' 0 1 $lastpart`
+echo "PARTs = $parts"
+
+for part in $parts; do
+ strpart=`printf "%03d" $part`
+ echo "PART = $part"
+ if [[ $fscr -ne 1 ]]; then
+  if [[ $antenna == 'lba' ]]; then
+   if [[ $part -eq $lastpart ]]; then
+     nchan=$((($nsubs-($lastpart * $subs_per_file))*$fscr)) 
+    else
+     nchan=$(($fscr * $subs_per_file))
+   fi
+   extra="-F ${nchan}:D"
+  else
+   extra=""
+  fi
+ fi
+ echo "dspsr ${extra} ${opts} -O ${stem}-P${part} $inputdir/${obsid}_SAP000_B000_S0_P${strpart}_bf.h5"
+ dspsr ${extra} ${opts} -O ${stem}-P${part} $inputdir/${obsid}_SAP000_B000_S0_P${strpart}_bf.h5
+done
+
+psradd -R -m time -o ${stem}.ar ${stem}-P*.ar
+
+# for this dataset 'paz' does a better job
+paz -r -e paz.ar ${stem}.ar
+
+# in this instance the clfd leaves more RFIs, and S/N is smaller
+clfd --no-report ${stem}.ar ; mv ${stem}.ar.clfd ${stem}.clfd.ar
+
+# making diagnostic plot
+if [[ $fscr -ne 1 ]]; then
+ ~/src/make-bf-diag.sh ${stem}.clfd.ar j S f${fscr}
+ ~/src/make-bf-diag.sh ${stem}.paz.ar j S f${fscr}
+else
+ ~/src/make-bf-diag.sh ${stem}.clfd.ar j S
+ ~/src/make-bf-diag.sh ${stem}.paz.ar j S
+fi
+#~/src/plot-bf-hdf5-spectrum.py $inputdir/${obsid}_SAP000_B000_S0_P000_bf.h5 -f ${fscr} --noflag --noshow
+
+pam -S -e stokes.uncal.ar ${stem}.paz.ar ${stem}.clfd.ar
+
+if [[ $fscr -ne 1 ]]; then
+ pam -f ${fscr} -e f${fscr}.ar ${stem}.paz.stokes.uncal.ar ${stem}.clfd.stokes.uncal.ar
+fi
+
+if [[ $fscr -ne 1 ]]; then
+ ~/src/pol-xxyy-stokes-uncal.py --psr ${psr} --noshow ${stem}.paz.stokes.uncal.f${fscr}.ar
+ ~/src/pol-xxyy-stokes-uncal.py --psr ${psr} --noshow ${stem}.clfd.stokes.uncal.f${fscr}.ar
+else
+ ~/src/pol-xxyy-stokes-uncal.py --psr ${psr} --noshow ${stem}.paz.stokes.uncal.ar
+ ~/src/pol-xxyy-stokes-uncal.py --psr ${psr} --noshow ${stem}.clfd.stokes.uncal.ar
+fi