diff --git a/bin/plotrapthor b/bin/plotrapthor
index e203a94186332cb9f2ac5154057e03f3304206d4..aa95d404d9474fa0986d70b69b503fb8c76422fe 100755
--- a/bin/plotrapthor
+++ b/bin/plotrapthor
@@ -7,23 +7,33 @@ from losoto.operations import plot
 import argparse
 from argparse import RawTextHelpFormatter
 from rapthor.lib import miscellaneous as misc
-import sys
-from losoto.lib_operations import *
+from losoto.lib_operations import normalize_phase, multiprocManager
 from losoto.operations.directionscreen import _calc_piercepoint
-from losoto._logging import logger as logging
+import logging
+import matplotlib as mpl
+mpl.use("Agg")
+import matplotlib.pyplot as plt  # after setting "Agg" to speed up
+from matplotlib.colors import ListedColormap
+from losoto.operations.stationscreen import _makeWCS, _circ_chi2
+from numpy import concatenate, newaxis
+from numpy.linalg import norm
+import numpy as np
+import os
+
+log = logging.getLogger('plotrapthor')
 
 
 def main(h5file, soltype, root=None, refstat=None, soltab=None, dir=None, ant=None,
          freq=None, pol='XX'):
     """
-    Plot solutions vs. time
+    Plot solutions
 
     Parameters
     ----------
     h5file : str
         Name of solution h5parm file
     soltype : str
-        Type of solution to plot: scalarphase, phase, amplitude, phasescreen, or ampscreen
+        Type of solution to plot: phase, amplitude, phasescreen, or ampscreen
     root : str, optional
         Root name for output plots. If None, the soltype is used
     refstat : str, optional
@@ -32,6 +42,14 @@ def main(h5file, soltype, root=None, refstat=None, soltab=None, dir=None, ant=No
         Name of soltab to use. If None, the default for the given soltype is used
     dir : str, optional
         Name of direction to use. If None, all directions are used
+    ant : list of str, optional
+        Name of antennas to use for screen plots. If None, RS210HBA is used
+    freq : float, optional
+        Frequency in Hz to use for screen plots. If None, the first frequency of
+        the input screen solutions is used
+    pol : str, optional
+        Name of the polarization to use for screen plots. If None, XX is used
+
     """
     if 'screen' in soltype:
         h = h5parm(h5file, readonly=False)
@@ -52,60 +70,62 @@ def main(h5file, soltype, root=None, refstat=None, soltab=None, dir=None, ant=No
                 st = ss.getSoltab(soltab)
                 st_resid = ss.getSoltab(soltab+'resid')
         else:
-            print('ERROR: solution type "{}" not understood. Must be one of scalarphase, '
-                  'phase, amplitude, phasescreen, or ampscreen'.format(soltype))
-            sys.exit(1)
+            raise ValueError('Solution type "{}" not understood. Must be one of phase, '
+                             'amplitude, phasescreen, or ampscreen'.format(soltype))
 
         if root is None:
             root = soltype + '_'
-        print('Plotting {} solutions...'.format(soltype))
+        log.info('Plotting {} solutions...'.format(soltype))
         if freq is None:
             freq = st.freq[0]
-            print('Frequency not specified, using {} MHz...'.format(freq/1e6))
+            log.info('Frequency not specified, using {} MHz...'.format(freq/1e6))
         if ant is None:
             ant = ['RS210HBA']
-            print('Antenna not specified, using {}...'.format(ant[0]))
+            log.info('Antenna not specified, using {}...'.format(ant[0]))
         else:
             ant = misc.string2list(ant)
         if soltype == 'phasescreen':
             st.setSelection(freq=freq, ant=ant)
             st_resid.setSelection(freq=freq, ant=ant)
             misc.remove_soltabs(ss, 'phase_screen000resid_filt')
-            tr_st = ss.makeSoltab('phase_screen', 'phase_screen000resid_filt',
-                                  axesNames=['time', 'freq', 'ant', 'dir'],
-                                  axesVals=[st_resid.time, st_resid.freq, st_resid.ant,
-                                            st_resid.dir],
-                                  vals=st_resid.val, weights=st_resid.weight)
+            _ = ss.makeSoltab('phase_screen', 'phase_screen000resid_filt',
+                              axesNames=['time', 'freq', 'ant', 'dir'],
+                              axesVals=[st_resid.time, st_resid.freq, st_resid.ant,
+                                        st_resid.dir],
+                              vals=st_resid.val, weights=st_resid.weight)
             runplots(st, prefix='phase', resSoltab='phase_screen000resid_filt',
-                           minZ=-3.1, maxZ=3.1)
+                     minZ=-3.1, maxZ=3.1)
             misc.remove_soltabs(ss, 'phase_screen000resid_filt')
         else:
             if pol is None:
                 pol = 'XX'
-                print('Polarization not specified, using {}...'.format(pol))
+                log.info('Polarization not specified, using {}...'.format(pol))
             st.setSelection(freq=freq, ant=ant, pol=pol)
             st_resid.setSelection(freq=freq, ant=ant, pol=pol)
             misc.remove_soltabs(ss, 'amplitude_screen000resid_filt')
-            tr_st = ss.makeSoltab('amplitude_screen', 'amplitude_screen000resid_filt',
-                                  axesNames=['time', 'freq', 'ant', 'dir', 'pol'],
-                                  axesVals=[st_resid.time, st_resid.freq, st_resid.ant,
-                                            st_resid.dir, st_resid.pol],
-                                  vals=st_resid.val, weights=st_resid.weight)
+            _ = ss.makeSoltab('amplitude_screen', 'amplitude_screen000resid_filt',
+                              axesNames=['time', 'freq', 'ant', 'dir', 'pol'],
+                              axesVals=[st_resid.time, st_resid.freq, st_resid.ant,
+                                        st_resid.dir, st_resid.pol],
+                              vals=st_resid.val, weights=st_resid.weight)
             runplots(st, prefix='amplitude', resSoltab='amplitude_screen000resid_filt',
-                           minZ=0.0, maxZ=2.0)
+                     minZ=0.0, maxZ=2.0)
             misc.remove_soltabs(ss, 'amplitude_screen000resid_filt')
         h.close()
     else:
         h = h5parm(h5file)
         ss = h.getSolset('sol000')
-        if soltype == 'scalarphase':
+        if soltype == 'phase':
             if soltab is None:
                 st = ss.getSoltab('phase000')
             else:
                 st = ss.getSoltab(soltab)
             ref = st.ant[0]
             ncol = 0
-            color = ''
+            if 'pol' in st.axesNames:
+                color = 'pol'
+            else:
+                color = ''
             minmax = [-3.2, 3.2]
         elif soltype == 'amplitude':
             if soltab is None:
@@ -114,31 +134,36 @@ def main(h5file, soltype, root=None, refstat=None, soltab=None, dir=None, ant=No
                 st = ss.getSoltab(soltab)
             ref = ''
             ncol = 0
-            color = 'pol'
-            minmax = [0, 0]
-        elif soltype == 'phase':
-            if soltab is None:
-                st = ss.getSoltab('phase000')
+            if 'pol' in st.axesNames:
+                color = 'pol'
             else:
-                st = ss.getSoltab(soltab)
-            ref = st.ant[0]
-            ncol = 0
-            color = 'pol'
-            minmax = [-3.2, 3.2]
+                color = ''
+            minmax = [0, 0]
         else:
-            print('ERROR: solution type "{}" not understood. Must be one of scalarphase, '
-                  'phase, amplitude, phasescreen, or ampscreen'.format(soltype))
-            sys.exit(1)
+            raise ValueError('Solution type "{}" not understood. Must be one of phase, '
+                             'amplitude, phasescreen, or ampscreen'.format(soltype))
 
         if root is None:
-            root = soltype + '_'
+            if 'pol' not in st.axesNames:
+                root = 'scalar' + soltype + '_'
+            else:
+                root = soltype + '_'
         if refstat is not None:
             ref = refstat
-        print('Plotting {} solutions...'.format(soltype))
+        log.info('Plotting {} solutions...'.format(soltype))
         if dir is not None:
             st.setSelection(dir=dir)
-        plot.run(st, ['time', 'freq'], axisInTable='ant', axisInCol=color, NColFig=ncol, refAnt=ref,
-                 prefix=root, minmax=minmax, plotFlag=True, markerSize=4)
+        if len(st.time) > 1 and len(st.freq) > 1:
+            plot.run(st, ['time', 'freq'], axisInTable='ant', axisInCol=color, NColFig=ncol, refAnt=ref,
+                     prefix=root, minmax=minmax, plotFlag=True, markerSize=4)
+        elif len(st.time) > 1:
+            plot.run(st, ['time'], axisInTable='ant', axisInCol=color, NColFig=ncol, refAnt=ref,
+                     prefix=root, minmax=minmax, plotFlag=True, markerSize=4)
+        elif len(st.freq) > 1:
+            plot.run(st, ['freq'], axisInTable='ant', axisInCol=color, NColFig=ncol, refAnt=ref,
+                     prefix=root, minmax=minmax, plotFlag=True, markerSize=4)
+        else:
+            log.warn('Solution table contains only a single time and frequency. No plots made.')
         h.close()
 
 
@@ -146,271 +171,269 @@ def _phase_cm():
     """
     Returns colormap for phase plots
     """
-    from matplotlib.colors import ListedColormap
-
-    cm_data = [[ 0.65830839, 0.46993917, 0.04941288],
-               [ 0.66433742, 0.4662019 , 0.05766473],
-               [ 0.67020869, 0.46248014, 0.0653456 ],
-               [ 0.67604299, 0.45869838, 0.07273174],
-               [ 0.68175228, 0.45491407, 0.07979262],
-               [ 0.6874028 , 0.45108417, 0.08667103],
-               [ 0.6929505 , 0.44723893, 0.09335869],
-               [ 0.69842619, 0.44335768, 0.09992839],
-               [ 0.7038123 , 0.43945328, 0.1063871 ],
-               [ 0.70912069, 0.43551765, 0.11277174],
-               [ 0.71434524, 0.43155576, 0.11909348],
-               [ 0.71949289, 0.42756272, 0.12537606],
-               [ 0.72455619, 0.4235447 , 0.13162325],
-               [ 0.72954895, 0.41949098, 0.13786305],
-               [ 0.73445172, 0.41541774, 0.14408039],
-               [ 0.73929496, 0.41129973, 0.15032217],
-               [ 0.74403834, 0.40717158, 0.15654335],
-               [ 0.74873695, 0.40298519, 0.16282282],
-               [ 0.75332319, 0.39880107, 0.16907566],
-               [ 0.75788083, 0.39454245, 0.17542179],
-               [ 0.7623326 , 0.39028096, 0.18175915],
-               [ 0.76673205, 0.38596549, 0.18816819],
-               [ 0.77105247, 0.38162141, 0.19461532],
-               [ 0.77529528, 0.37724732, 0.20110652],
-               [ 0.77948666, 0.37281509, 0.2076873 ],
-               [ 0.78358534, 0.36836772, 0.21429736],
-               [ 0.78763763, 0.363854  , 0.22101648],
-               [ 0.79161134, 0.35930804, 0.2277974 ],
-               [ 0.79550606, 0.3547299 , 0.23464353],
-               [ 0.79935398, 0.35007959, 0.24161832],
-               [ 0.80311671, 0.34540152, 0.24865892],
-               [ 0.80681033, 0.34067452, 0.25580075],
-               [ 0.8104452 , 0.33588248, 0.26307222],
-               [ 0.8139968 , 0.33105538, 0.27043183],
-               [ 0.81747689, 0.32617526, 0.27791096],
-               [ 0.82089415, 0.32122629, 0.28553846],
-               [ 0.82422713, 0.3162362 , 0.29327617],
-               [ 0.82747661, 0.31120154, 0.30113388],
-               [ 0.83066399, 0.30608459, 0.30917579],
-               [ 0.83376307, 0.30092244, 0.31734921],
-               [ 0.83677286, 0.29571346, 0.32566199],
-               [ 0.83969693, 0.29044723, 0.33413665],
-               [ 0.84253873, 0.28511151, 0.34279962],
-               [ 0.84528297, 0.27972917, 0.35162078],
-               [ 0.84792704, 0.27430045, 0.36060681],
-               [ 0.85046793, 0.26882624, 0.36976395],
-               [ 0.85291056, 0.26328859, 0.37913116],
-               [ 0.855242  , 0.25770888, 0.38868217],
-               [ 0.85745673, 0.25209367, 0.39841601],
-               [ 0.85955023, 0.24644737, 0.40833625],
-               [ 0.86151767, 0.24077563, 0.41844557],
-               [ 0.86335392, 0.23508521, 0.42874606],
-               [ 0.86505685, 0.22937288, 0.43926008],
-               [ 0.86661606, 0.22366308, 0.44996127],
-               [ 0.86802578, 0.21796785, 0.46084758],
-               [ 0.86928003, 0.21230132, 0.47191554],
-               [ 0.87037274, 0.20667988, 0.48316015],
-               [ 0.87129781, 0.2011224 , 0.49457479],
-               [ 0.87204914, 0.19565041, 0.50615118],
-               [ 0.87262076, 0.19028829, 0.51787932],
-               [ 0.87300686, 0.18506334, 0.5297475 ],
-               [ 0.8732019 , 0.18000588, 0.54174232],
-               [ 0.87320066, 0.1751492 , 0.55384874],
-               [ 0.87299833, 0.17052942, 0.56605016],
-               [ 0.87259058, 0.16618514, 0.57832856],
-               [ 0.87197361, 0.16215698, 0.59066466],
-               [ 0.87114414, 0.15848667, 0.60303881],
-               [ 0.87009966, 0.15521687, 0.61542844],
-               [ 0.86883823, 0.15238892, 0.62781175],
-               [ 0.86735858, 0.15004199, 0.64016651],
-               [ 0.8656601 , 0.14821149, 0.65247022],
-               [ 0.86374282, 0.14692762, 0.66470043],
-               [ 0.86160744, 0.14621386, 0.67683495],
-               [ 0.85925523, 0.14608582, 0.68885204],
-               [ 0.85668805, 0.14655046, 0.70073065],
-               [ 0.85390829, 0.14760576, 0.71245054],
-               [ 0.85091881, 0.14924094, 0.7239925 ],
-               [ 0.84772287, 0.15143717, 0.73533849],
-               [ 0.84432409, 0.15416865, 0.74647174],
-               [ 0.84072639, 0.15740403, 0.75737678],
-               [ 0.83693394, 0.16110786, 0.76803952],
-               [ 0.83295108, 0.16524205, 0.77844723],
-               [ 0.82878232, 0.16976729, 0.78858858],
-               [ 0.82443225, 0.17464414, 0.7984536 ],
-               [ 0.81990551, 0.179834  , 0.80803365],
-               [ 0.81520674, 0.18529984, 0.8173214 ],
-               [ 0.81034059, 0.19100664, 0.82631073],
-               [ 0.80531176, 0.1969216 , 0.83499645],
-               [ 0.80012467, 0.20301465, 0.84337486],
-               [ 0.79478367, 0.20925826, 0.8514432 ],
-               [ 0.78929302, 0.21562737, 0.85919957],
-               [ 0.78365681, 0.22209936, 0.86664294],
-               [ 0.77787898, 0.22865386, 0.87377308],
-               [ 0.7719633 , 0.23527265, 0.88059043],
-               [ 0.76591335, 0.24193947, 0.88709606],
-               [ 0.7597325 , 0.24863985, 0.89329158],
-               [ 0.75342394, 0.25536094, 0.89917908],
-               [ 0.74699063, 0.26209137, 0.90476105],
-               [ 0.74043533, 0.2688211 , 0.91004033],
-               [ 0.73376055, 0.27554128, 0.91502   ],
-               [ 0.72696862, 0.28224415, 0.91970339],
-               [ 0.7200616 , 0.2889229 , 0.92409395],
-               [ 0.71304134, 0.29557159, 0.92819525],
-               [ 0.70590945, 0.30218508, 0.9320109 ],
-               [ 0.69866732, 0.30875887, 0.93554451],
-               [ 0.69131609, 0.31528914, 0.93879964],
-               [ 0.68385669, 0.32177259, 0.94177976],
-               [ 0.6762898 , 0.32820641, 0.94448822],
-               [ 0.6686159 , 0.33458824, 0.94692818],
-               [ 0.66083524, 0.3409161 , 0.94910264],
-               [ 0.65294785, 0.34718834, 0.95101432],
-               [ 0.64495358, 0.35340362, 0.95266571],
-               [ 0.63685208, 0.35956083, 0.954059  ],
-               [ 0.62864284, 0.3656591 , 0.95519608],
-               [ 0.62032517, 0.3716977 , 0.95607853],
-               [ 0.61189825, 0.37767607, 0.95670757],
-               [ 0.60336117, 0.38359374, 0.95708408],
-               [ 0.59471291, 0.3894503 , 0.95720861],
-               [ 0.58595242, 0.39524541, 0.95708134],
-               [ 0.5770786 , 0.40097871, 0.95670212],
-               [ 0.56809041, 0.40664983, 0.95607045],
-               [ 0.55898686, 0.41225834, 0.95518556],
-               [ 0.54976709, 0.41780374, 0.95404636],
-               [ 0.5404304 , 0.42328541, 0.95265153],
-               [ 0.53097635, 0.42870263, 0.95099953],
-               [ 0.52140479, 0.43405447, 0.94908866],
-               [ 0.51171597, 0.43933988, 0.94691713],
-               [ 0.50191056, 0.44455757, 0.94448311],
-               [ 0.49198981, 0.44970607, 0.94178481],
-               [ 0.48195555, 0.45478367, 0.93882055],
-               [ 0.47181035, 0.45978843, 0.93558888],
-               [ 0.46155756, 0.46471821, 0.93208866],
-               [ 0.45119801, 0.46957218, 0.92831786],
-               [ 0.44073852, 0.47434688, 0.92427669],
-               [ 0.43018722, 0.47903864, 0.9199662 ],
-               [ 0.41955166, 0.4836444 , 0.91538759],
-               [ 0.40884063, 0.48816094, 0.91054293],
-               [ 0.39806421, 0.49258494, 0.90543523],
-               [ 0.38723377, 0.49691301, 0.90006852],
-               [ 0.37636206, 0.50114173, 0.89444794],
-               [ 0.36546127, 0.5052684 , 0.88857877],
-               [ 0.35454654, 0.5092898 , 0.88246819],
-               [ 0.34363779, 0.51320158, 0.87612664],
-               [ 0.33275309, 0.51700082, 0.86956409],
-               [ 0.32191166, 0.52068487, 0.86279166],
-               [ 0.31113372, 0.52425144, 0.85582152],
-               [ 0.3004404 , 0.52769862, 0.84866679],
-               [ 0.28985326, 0.53102505, 0.84134123],
-               [ 0.27939616, 0.53422931, 0.83386051],
-               [ 0.26909181, 0.53731099, 0.82623984],
-               [ 0.258963  , 0.5402702 , 0.81849475],
-               [ 0.24903239, 0.54310763, 0.8106409 ],
-               [ 0.23932229, 0.54582448, 0.80269392],
-               [ 0.22985664, 0.54842189, 0.79467122],
-               [ 0.2206551 , 0.55090241, 0.78658706],
-               [ 0.21173641, 0.55326901, 0.77845533],
-               [ 0.20311843, 0.55552489, 0.77028973],
-               [ 0.1948172 , 0.55767365, 0.76210318],
-               [ 0.1868466 , 0.55971922, 0.75390763],
-               [ 0.17921799, 0.56166586, 0.74571407],
-               [ 0.1719422 , 0.56351747, 0.73753498],
-               [ 0.16502295, 0.56527915, 0.72937754],
-               [ 0.15846116, 0.566956  , 0.72124819],
-               [ 0.15225499, 0.56855297, 0.71315321],
-               [ 0.14639876, 0.57007506, 0.70509769],
-               [ 0.14088284, 0.57152729, 0.69708554],
-               [ 0.13569366, 0.57291467, 0.68911948],
-               [ 0.13081385, 0.57424211, 0.68120108],
-               [ 0.12622247, 0.57551447, 0.67333078],
-               [ 0.12189539, 0.57673644, 0.66550792],
-               [ 0.11780654, 0.57791235, 0.65773233],
-               [ 0.11392613, 0.5790468 , 0.64999984],
-               [ 0.11022348, 0.58014398, 0.64230637],
-               [ 0.10666732, 0.58120782, 0.63464733],
-               [ 0.10322631, 0.58224198, 0.62701729],
-               [ 0.0998697 , 0.58324982, 0.61941001],
-               [ 0.09656813, 0.58423445, 0.61181853],
-               [ 0.09329429, 0.58519864, 0.60423523],
-               [ 0.09002364, 0.58614483, 0.5966519 ],
-               [ 0.08673514, 0.58707512, 0.58905979],
-               [ 0.08341199, 0.58799127, 0.58144971],
-               [ 0.08004245, 0.58889466, 0.57381211],
-               [ 0.07662083, 0.58978633, 0.56613714],
-               [ 0.07314852, 0.59066692, 0.55841474],
-               [ 0.06963541, 0.5915367 , 0.55063471],
-               [ 0.06610144, 0.59239556, 0.54278681],
-               [ 0.06257861, 0.59324304, 0.53486082],
-               [ 0.05911304, 0.59407833, 0.52684614],
-               [ 0.05576765, 0.5949003 , 0.5187322 ],
-               [ 0.05262511, 0.59570732, 0.51050978],
-               [ 0.04978881, 0.5964975 , 0.50216936],
-               [ 0.04738319, 0.59726862, 0.49370174],
-               [ 0.04555067, 0.59801813, 0.48509809],
-               [ 0.04444396, 0.59874316, 0.47635   ],
-               [ 0.04421323, 0.59944056, 0.46744951],
-               [ 0.04498918, 0.60010687, 0.45838913],
-               [ 0.04686604, 0.60073837, 0.44916187],
-               [ 0.04988979, 0.60133103, 0.43976125],
-               [ 0.05405573, 0.60188055, 0.4301812 ],
-               [ 0.05932209, 0.60238289, 0.42040543],
-               [ 0.06560774, 0.60283258, 0.41043772],
-               [ 0.07281962, 0.60322442, 0.40027363],
-               [ 0.08086177, 0.60355283, 0.38990941],
-               [ 0.08964366, 0.60381194, 0.37934208],
-               [ 0.09908952, 0.60399554, 0.36856412],
-               [ 0.10914617, 0.60409695, 0.35755799],
-               [ 0.11974119, 0.60410858, 0.34634096],
-               [ 0.13082746, 0.6040228 , 0.33491416],
-               [ 0.14238003, 0.60383119, 0.323267  ],
-               [ 0.1543847 , 0.60352425, 0.31138823],
-               [ 0.16679093, 0.60309301, 0.29931029],
-               [ 0.17959757, 0.60252668, 0.2870237 ],
-               [ 0.19279966, 0.60181364, 0.27452964],
-               [ 0.20634465, 0.60094466, 0.2618794 ],
-               [ 0.22027287, 0.5999043 , 0.24904251],
-               [ 0.23449833, 0.59868591, 0.23611022],
-               [ 0.24904416, 0.5972746 , 0.2230778 ],
-               [ 0.26382006, 0.59566656, 0.21004673],
-               [ 0.2788104 , 0.5938521 , 0.19705484],
-               [ 0.29391494, 0.59183348, 0.18421621],
-               [ 0.3090634 , 0.58961302, 0.17161942],
-               [ 0.32415577, 0.58720132, 0.15937753],
-               [ 0.3391059 , 0.58461164, 0.14759012],
-               [ 0.35379624, 0.58186793, 0.13637734],
-               [ 0.36817905, 0.5789861 , 0.12580054],
-               [ 0.38215966, 0.57599512, 0.1159504 ],
-               [ 0.39572824, 0.57290928, 0.10685038],
-               [ 0.40881926, 0.56975727, 0.09855521],
-               [ 0.42148106, 0.56654159, 0.09104002],
-               [ 0.43364953, 0.56329296, 0.08434116],
-               [ 0.44538908, 0.56000859, 0.07841305],
-               [ 0.45672421, 0.5566943 , 0.07322913],
-               [ 0.46765017, 0.55336373, 0.06876762],
-               [ 0.47819138, 0.5500213 , 0.06498436],
-               [ 0.48839686, 0.54666195, 0.06182163],
-               [ 0.49828924, 0.5432874 , 0.05922726],
-               [ 0.50789114, 0.53989827, 0.05714466],
-               [ 0.51722475, 0.53649429, 0.05551476],
-               [ 0.5263115 , 0.53307443, 0.05427793],
-               [ 0.53517186, 0.52963707, 0.05337567],
-               [ 0.54382515, 0.52618009, 0.05275208],
-               [ 0.55228947, 0.52270103, 0.05235479],
-               [ 0.56058163, 0.51919713, 0.0521356 ],
-               [ 0.56871719, 0.51566545, 0.05205062],
-               [ 0.57671045, 0.51210292, 0.0520602 ],
-               [ 0.5845745 , 0.50850636, 0.05212851],
-               [ 0.59232129, 0.50487256, 0.05222299],
-               [ 0.5999617 , 0.50119827, 0.05231367],
-               [ 0.60750568, 0.49748022, 0.05237234],
-               [ 0.61496232, 0.49371512, 0.05237168],
-               [ 0.62233999, 0.48989963, 0.05228423],
-               [ 0.62964652, 0.48603032, 0.05208127],
-               [ 0.63688935, 0.48210362, 0.05173155],
-               [ 0.64407572, 0.4781157 , 0.0511996 ],
-               [ 0.65121289, 0.47406244, 0.05044367],
-               [ 0.65830839, 0.46993917, 0.04941288]]
+    cm_data = [[0.65830839, 0.46993917, 0.04941288],
+               [0.66433742, 0.46620190, 0.05766473],
+               [0.67020869, 0.46248014, 0.06534560],
+               [0.67604299, 0.45869838, 0.07273174],
+               [0.68175228, 0.45491407, 0.07979262],
+               [0.68740280, 0.45108417, 0.08667103],
+               [0.69295050, 0.44723893, 0.09335869],
+               [0.69842619, 0.44335768, 0.09992839],
+               [0.70381230, 0.43945328, 0.10638710],
+               [0.70912069, 0.43551765, 0.11277174],
+               [0.71434524, 0.43155576, 0.11909348],
+               [0.71949289, 0.42756272, 0.12537606],
+               [0.72455619, 0.42354470, 0.13162325],
+               [0.72954895, 0.41949098, 0.13786305],
+               [0.73445172, 0.41541774, 0.14408039],
+               [0.73929496, 0.41129973, 0.15032217],
+               [0.74403834, 0.40717158, 0.15654335],
+               [0.74873695, 0.40298519, 0.16282282],
+               [0.75332319, 0.39880107, 0.16907566],
+               [0.75788083, 0.39454245, 0.17542179],
+               [0.76233260, 0.39028096, 0.18175915],
+               [0.76673205, 0.38596549, 0.18816819],
+               [0.77105247, 0.38162141, 0.19461532],
+               [0.77529528, 0.37724732, 0.20110652],
+               [0.77948666, 0.37281509, 0.20768730],
+               [0.78358534, 0.36836772, 0.21429736],
+               [0.78763763, 0.36385400, 0.22101648],
+               [0.79161134, 0.35930804, 0.22779740],
+               [0.79550606, 0.35472990, 0.23464353],
+               [0.79935398, 0.35007959, 0.24161832],
+               [0.80311671, 0.34540152, 0.24865892],
+               [0.80681033, 0.34067452, 0.25580075],
+               [0.81044520, 0.33588248, 0.26307222],
+               [0.81399680, 0.33105538, 0.27043183],
+               [0.81747689, 0.32617526, 0.27791096],
+               [0.82089415, 0.32122629, 0.28553846],
+               [0.82422713, 0.31623620, 0.29327617],
+               [0.82747661, 0.31120154, 0.30113388],
+               [0.83066399, 0.30608459, 0.30917579],
+               [0.83376307, 0.30092244, 0.31734921],
+               [0.83677286, 0.29571346, 0.32566199],
+               [0.83969693, 0.29044723, 0.33413665],
+               [0.84253873, 0.28511151, 0.34279962],
+               [0.84528297, 0.27972917, 0.35162078],
+               [0.84792704, 0.27430045, 0.36060681],
+               [0.85046793, 0.26882624, 0.36976395],
+               [0.85291056, 0.26328859, 0.37913116],
+               [0.85524200, 0.25770888, 0.38868217],
+               [0.85745673, 0.25209367, 0.39841601],
+               [0.85955023, 0.24644737, 0.40833625],
+               [0.86151767, 0.24077563, 0.41844557],
+               [0.86335392, 0.23508521, 0.42874606],
+               [0.86505685, 0.22937288, 0.43926008],
+               [0.86661606, 0.22366308, 0.44996127],
+               [0.86802578, 0.21796785, 0.46084758],
+               [0.86928003, 0.21230132, 0.47191554],
+               [0.87037274, 0.20667988, 0.48316015],
+               [0.87129781, 0.20112240, 0.49457479],
+               [0.87204914, 0.19565041, 0.50615118],
+               [0.87262076, 0.19028829, 0.51787932],
+               [0.87300686, 0.18506334, 0.52974750],
+               [0.87320190, 0.18000588, 0.54174232],
+               [0.87320066, 0.17514920, 0.55384874],
+               [0.87299833, 0.17052942, 0.56605016],
+               [0.87259058, 0.16618514, 0.57832856],
+               [0.87197361, 0.16215698, 0.59066466],
+               [0.87114414, 0.15848667, 0.60303881],
+               [0.87009966, 0.15521687, 0.61542844],
+               [0.86883823, 0.15238892, 0.62781175],
+               [0.86735858, 0.15004199, 0.64016651],
+               [0.86566010, 0.14821149, 0.65247022],
+               [0.86374282, 0.14692762, 0.66470043],
+               [0.86160744, 0.14621386, 0.67683495],
+               [0.85925523, 0.14608582, 0.68885204],
+               [0.85668805, 0.14655046, 0.70073065],
+               [0.85390829, 0.14760576, 0.71245054],
+               [0.85091881, 0.14924094, 0.72399250],
+               [0.84772287, 0.15143717, 0.73533849],
+               [0.84432409, 0.15416865, 0.74647174],
+               [0.84072639, 0.15740403, 0.75737678],
+               [0.83693394, 0.16110786, 0.76803952],
+               [0.83295108, 0.16524205, 0.77844723],
+               [0.82878232, 0.16976729, 0.78858858],
+               [0.82443225, 0.17464414, 0.79845360],
+               [0.81990551, 0.17983400, 0.80803365],
+               [0.81520674, 0.18529984, 0.81732140],
+               [0.81034059, 0.19100664, 0.82631073],
+               [0.80531176, 0.19692160, 0.83499645],
+               [0.80012467, 0.20301465, 0.84337486],
+               [0.79478367, 0.20925826, 0.85144320],
+               [0.78929302, 0.21562737, 0.85919957],
+               [0.78365681, 0.22209936, 0.86664294],
+               [0.77787898, 0.22865386, 0.87377308],
+               [0.77196330, 0.23527265, 0.88059043],
+               [0.76591335, 0.24193947, 0.88709606],
+               [0.75973250, 0.24863985, 0.89329158],
+               [0.75342394, 0.25536094, 0.89917908],
+               [0.74699063, 0.26209137, 0.90476105],
+               [0.74043533, 0.26882110, 0.91004033],
+               [0.73376055, 0.27554128, 0.91502000],
+               [0.72696862, 0.28224415, 0.91970339],
+               [0.72006160, 0.28892290, 0.92409395],
+               [0.71304134, 0.29557159, 0.92819525],
+               [0.70590945, 0.30218508, 0.93201090],
+               [0.69866732, 0.30875887, 0.93554451],
+               [0.69131609, 0.31528914, 0.93879964],
+               [0.68385669, 0.32177259, 0.94177976],
+               [0.67628980, 0.32820641, 0.94448822],
+               [0.66861590, 0.33458824, 0.94692818],
+               [0.66083524, 0.34091610, 0.94910264],
+               [0.65294785, 0.34718834, 0.95101432],
+               [0.64495358, 0.35340362, 0.95266571],
+               [0.63685208, 0.35956083, 0.95405900],
+               [0.62864284, 0.36565910, 0.95519608],
+               [0.62032517, 0.37169770, 0.95607853],
+               [0.61189825, 0.37767607, 0.95670757],
+               [0.60336117, 0.38359374, 0.95708408],
+               [0.59471291, 0.38945030, 0.95720861],
+               [0.58595242, 0.39524541, 0.95708134],
+               [0.57707860, 0.40097871, 0.95670212],
+               [0.56809041, 0.40664983, 0.95607045],
+               [0.55898686, 0.41225834, 0.95518556],
+               [0.54976709, 0.41780374, 0.95404636],
+               [0.54043040, 0.42328541, 0.95265153],
+               [0.53097635, 0.42870263, 0.95099953],
+               [0.52140479, 0.43405447, 0.94908866],
+               [0.51171597, 0.43933988, 0.94691713],
+               [0.50191056, 0.44455757, 0.94448311],
+               [0.49198981, 0.44970607, 0.94178481],
+               [0.48195555, 0.45478367, 0.93882055],
+               [0.47181035, 0.45978843, 0.93558888],
+               [0.46155756, 0.46471821, 0.93208866],
+               [0.45119801, 0.46957218, 0.92831786],
+               [0.44073852, 0.47434688, 0.92427669],
+               [0.43018722, 0.47903864, 0.91996620],
+               [0.41955166, 0.48364440, 0.91538759],
+               [0.40884063, 0.48816094, 0.91054293],
+               [0.39806421, 0.49258494, 0.90543523],
+               [0.38723377, 0.49691301, 0.90006852],
+               [0.37636206, 0.50114173, 0.89444794],
+               [0.36546127, 0.50526840, 0.88857877],
+               [0.35454654, 0.50928980, 0.88246819],
+               [0.34363779, 0.51320158, 0.87612664],
+               [0.33275309, 0.51700082, 0.86956409],
+               [0.32191166, 0.52068487, 0.86279166],
+               [0.31113372, 0.52425144, 0.85582152],
+               [0.30044040, 0.52769862, 0.84866679],
+               [0.28985326, 0.53102505, 0.84134123],
+               [0.27939616, 0.53422931, 0.83386051],
+               [0.26909181, 0.53731099, 0.82623984],
+               [0.25896300, 0.54027020, 0.81849475],
+               [0.24903239, 0.54310763, 0.81064090],
+               [0.23932229, 0.54582448, 0.80269392],
+               [0.22985664, 0.54842189, 0.79467122],
+               [0.22065510, 0.55090241, 0.78658706],
+               [0.21173641, 0.55326901, 0.77845533],
+               [0.20311843, 0.55552489, 0.77028973],
+               [0.19481720, 0.55767365, 0.76210318],
+               [0.18684660, 0.55971922, 0.75390763],
+               [0.17921799, 0.56166586, 0.74571407],
+               [0.17194220, 0.56351747, 0.73753498],
+               [0.16502295, 0.56527915, 0.72937754],
+               [0.15846116, 0.56695600, 0.72124819],
+               [0.15225499, 0.56855297, 0.71315321],
+               [0.14639876, 0.57007506, 0.70509769],
+               [0.14088284, 0.57152729, 0.69708554],
+               [0.13569366, 0.57291467, 0.68911948],
+               [0.13081385, 0.57424211, 0.68120108],
+               [0.12622247, 0.57551447, 0.67333078],
+               [0.12189539, 0.57673644, 0.66550792],
+               [0.11780654, 0.57791235, 0.65773233],
+               [0.11392613, 0.57904680, 0.64999984],
+               [0.11022348, 0.58014398, 0.64230637],
+               [0.10666732, 0.58120782, 0.63464733],
+               [0.10322631, 0.58224198, 0.62701729],
+               [0.09986970, 0.58324982, 0.61941001],
+               [0.09656813, 0.58423445, 0.61181853],
+               [0.09329429, 0.58519864, 0.60423523],
+               [0.09002364, 0.58614483, 0.59665190],
+               [0.08673514, 0.58707512, 0.58905979],
+               [0.08341199, 0.58799127, 0.58144971],
+               [0.08004245, 0.58889466, 0.57381211],
+               [0.07662083, 0.58978633, 0.56613714],
+               [0.07314852, 0.59066692, 0.55841474],
+               [0.06963541, 0.59153670, 0.55063471],
+               [0.06610144, 0.59239556, 0.54278681],
+               [0.06257861, 0.59324304, 0.53486082],
+               [0.05911304, 0.59407833, 0.52684614],
+               [0.05576765, 0.59490030, 0.51873220],
+               [0.05262511, 0.59570732, 0.51050978],
+               [0.04978881, 0.59649750, 0.50216936],
+               [0.04738319, 0.59726862, 0.49370174],
+               [0.04555067, 0.59801813, 0.48509809],
+               [0.04444396, 0.59874316, 0.47635000],
+               [0.04421323, 0.59944056, 0.46744951],
+               [0.04498918, 0.60010687, 0.45838913],
+               [0.04686604, 0.60073837, 0.44916187],
+               [0.04988979, 0.60133103, 0.43976125],
+               [0.05405573, 0.60188055, 0.43018120],
+               [0.05932209, 0.60238289, 0.42040543],
+               [0.06560774, 0.60283258, 0.41043772],
+               [0.07281962, 0.60322442, 0.40027363],
+               [0.08086177, 0.60355283, 0.38990941],
+               [0.08964366, 0.60381194, 0.37934208],
+               [0.09908952, 0.60399554, 0.36856412],
+               [0.10914617, 0.60409695, 0.35755799],
+               [0.11974119, 0.60410858, 0.34634096],
+               [0.13082746, 0.60402280, 0.33491416],
+               [0.14238003, 0.60383119, 0.32326700],
+               [0.15438470, 0.60352425, 0.31138823],
+               [0.16679093, 0.60309301, 0.29931029],
+               [0.17959757, 0.60252668, 0.28702370],
+               [0.19279966, 0.60181364, 0.27452964],
+               [0.20634465, 0.60094466, 0.26187940],
+               [0.22027287, 0.59990430, 0.24904251],
+               [0.23449833, 0.59868591, 0.23611022],
+               [0.24904416, 0.59727460, 0.22307780],
+               [0.26382006, 0.59566656, 0.21004673],
+               [0.27881040, 0.59385210, 0.19705484],
+               [0.29391494, 0.59183348, 0.18421621],
+               [0.30906340, 0.58961302, 0.17161942],
+               [0.32415577, 0.58720132, 0.15937753],
+               [0.33910590, 0.58461164, 0.14759012],
+               [0.35379624, 0.58186793, 0.13637734],
+               [0.36817905, 0.57898610, 0.12580054],
+               [0.38215966, 0.57599512, 0.11595040],
+               [0.39572824, 0.57290928, 0.10685038],
+               [0.40881926, 0.56975727, 0.09855521],
+               [0.42148106, 0.56654159, 0.09104002],
+               [0.43364953, 0.56329296, 0.08434116],
+               [0.44538908, 0.56000859, 0.07841305],
+               [0.45672421, 0.55669430, 0.07322913],
+               [0.46765017, 0.55336373, 0.06876762],
+               [0.47819138, 0.55002130, 0.06498436],
+               [0.48839686, 0.54666195, 0.06182163],
+               [0.49828924, 0.54328740, 0.05922726],
+               [0.50789114, 0.53989827, 0.05714466],
+               [0.51722475, 0.53649429, 0.05551476],
+               [0.52631150, 0.53307443, 0.05427793],
+               [0.53517186, 0.52963707, 0.05337567],
+               [0.54382515, 0.52618009, 0.05275208],
+               [0.55228947, 0.52270103, 0.05235479],
+               [0.56058163, 0.51919713, 0.05213560],
+               [0.56871719, 0.51566545, 0.05205062],
+               [0.57671045, 0.51210292, 0.05206020],
+               [0.58457450, 0.50850636, 0.05212851],
+               [0.59232129, 0.50487256, 0.05222299],
+               [0.59996170, 0.50119827, 0.05231367],
+               [0.60750568, 0.49748022, 0.05237234],
+               [0.61496232, 0.49371512, 0.05237168],
+               [0.62233999, 0.48989963, 0.05228423],
+               [0.62964652, 0.48603032, 0.05208127],
+               [0.63688935, 0.48210362, 0.05173155],
+               [0.64407572, 0.47811570, 0.05119960],
+               [0.65121289, 0.47406244, 0.05044367],
+               [0.65830839, 0.46993917, 0.04941288]]
     cm = ListedColormap(cm_data, name=__file__)
 
     return cm
 
 
 def _calculate_screen(inscreen, residuals, pp, N_piercepoints, k, east, north, up,
-    T, Nx, Ny, sindx, height, beta_val, r_0, is_phase, outQueue):
+                      T, Nx, Ny, sindx, height, beta_val, r_0, is_phase, outQueue):
     """
     Calculates screen images
 
@@ -451,10 +474,6 @@ def _calculate_screen(inscreen, residuals, pp, N_piercepoints, k, east, north, u
         input screen is a phase screen
 
     """
-    from numpy import kron, concatenate, newaxis
-    from numpy.linalg import pinv, norm
-    import numpy as np
-
     screen = np.zeros((Nx, Ny))
 
     if height == 0.0:
@@ -491,7 +510,7 @@ def _calculate_screen(inscreen, residuals, pp, N_piercepoints, k, east, north, u
             else:
                 p, airmass = _calc_piercepoint(np.dot(np.array([xi, yi]), np.array([east, north])), up, height)
             d2 = np.sum(np.square(pp - p), axis=1)
-            c = -(d2 / ( r_0**2 ))**(beta_val / 2.0) / 2.0
+            c = -(d2 / (r_0**2))**(beta_val / 2.0) / 2.0
             if is_phase:
                 screen[i, j] = np.dot(c, f)
             else:
@@ -502,15 +521,15 @@ def _calculate_screen(inscreen, residuals, pp, N_piercepoints, k, east, north, u
         x = pp1[:, 0]
         y = pp1[:, 1]
     else:
-        x = (pp1[:, 0] / 1000.0 - min_xy[0] / 1000.0) / extent[0] # km
-        y = (pp1[:, 1] / 1000.0 - min_xy[1] / 1000.0) / extent[1] # km
+        x = (pp1[:, 0] / 1000.0 - min_xy[0] / 1000.0) / extent[0]  # km
+        y = (pp1[:, 1] / 1000.0 - min_xy[1] / 1000.0) / extent[1]  # km
 
     outQueue.put([k, fitted_vals, screen, x, y])
 
 
 def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
-    upper, vmin, vmax, source_names, show_source_names, station_names, sindx,
-    root_dir, prestr, is_image_plane,  midRA, midDec, order, is_phase, outQueue):
+                upper, vmin, vmax, source_names, show_source_names, station_names, sindx,
+                root_dir, prestr, is_image_plane,  midRA, midDec, order, is_phase, outQueue):
     """
     Plots screen images
 
@@ -548,38 +567,18 @@ def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
         True if screens are phase screens
 
     """
-    if not 'matplotlib' in sys.modules:
-        import matplotlib as mpl
-        mpl.rc('figure.subplot',left=0.05, bottom=0.05, right=0.95, top=0.95,wspace=0.22, hspace=0.22 )
-        mpl.use("Agg")
-    import matplotlib as mpl
-    import matplotlib.pyplot as plt # after setting "Agg" to speed up
-    from losoto.operations.stationscreen import _makeWCS, _circ_chi2
-    import numpy as np
-    try:
-        try:
-            from astropy.visualization.wcsaxes import WCSAxes
-            hasWCSaxes = True
-        except:
-            from wcsaxes import WCSAxes
-            hasWCSaxes = True
-    except:
-        hasWCSaxes = False
-    from matplotlib.colors import LinearSegmentedColormap
-
-
-    fig = plt.figure(figsize=(6,6))
+    mpl.rc('figure.subplot', left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.22, hspace=0.22)
+    fig = plt.figure(figsize=(6, 6))
 
     # Set colormap
     if is_phase:
         cmap = _phase_cm()
     else:
         cmap = plt.cm.jet
-    sm = plt.cm.ScalarMappable(cmap=cmap,
-        norm=mpl.colors.Normalize(vmin=vmin, vmax=vmax))
+    sm = plt.cm.ScalarMappable(cmap=cmap, norm=mpl.colors.Normalize(vmin=vmin, vmax=vmax))
     sm._A = []
 
-    if is_image_plane and hasWCSaxes:
+    if is_image_plane:
         wcs = _makeWCS(midRA, midDec)
         ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=wcs)
     else:
@@ -622,12 +621,9 @@ def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
         lower /= 1000.0
         upper /= 1000.0
 
-    im = ax.imshow(screen.transpose([1, 0])[:, :],
-        cmap = cmap,
-        origin = 'lower',
-        interpolation = 'nearest',
-        extent = (lower[0], upper[0], lower[1], upper[1]),
-        vmin=vmin, vmax=vmax)
+    im = ax.imshow(screen.transpose([1, 0])[:, :], cmap=cmap, origin='lower',
+                   interpolation='nearest', extent=(lower[0], upper[0], lower[1], upper[1]),
+                   vmin=vmin, vmax=vmax)
 
     cbar = plt.colorbar(im)
     cbar.set_label('Value', rotation=270)
@@ -638,16 +634,14 @@ def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
     if show_source_names:
         labels = source_names
         for label, xl, yl in zip(labels, x, y):
-            plt.annotate(
-                label,
-                xy = (xl, yl), xytext = (-2, 2),
-                textcoords = 'offset points', ha = 'right', va = 'bottom')
+            plt.annotate(label, xy=(xl, yl), xytext=(-2, 2), textcoords='offset points',
+                         ha='right', va='bottom')
 
     nsrcs = np.where(weights > 0.0)[0].size
     if is_phase:
-        redchi2 =  _circ_chi2(residuals, weights) / (nsrcs-order)
+        redchi2 = _circ_chi2(residuals, weights) / (nsrcs-order)
     else:
-        redchi2 =  np.sum(np.square(residuals) * weights) / (nsrcs-order)
+        redchi2 = np.sum(np.square(residuals) * weights) / (nsrcs-order)
     if sindx >= 0:
         plt.title('Station {0}, Time {1} (red. chi2 = {2:0.3f})'.format(station_names[sindx], k, redchi2))
     else:
@@ -656,19 +650,13 @@ def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
         ax.set_xlim(lower[0], upper[0])
         ax.set_ylim(lower[1], upper[1])
         ax.set_aspect('equal')
-        if hasWCSaxes:
-            RAAxis = ax.coords['ra']
-            RAAxis.set_axislabel('RA', minpad=0.75)
-            RAAxis.set_major_formatter('hh:mm:ss')
-            DecAxis = ax.coords['dec']
-            DecAxis.set_axislabel('Dec', minpad=0.75)
-            DecAxis.set_major_formatter('dd:mm:ss')
-            ax.coords.grid(color='black', alpha=0.5, linestyle='solid')
-            # plt.xlabel("RA")
-            # plt.ylabel("Dec")
-        else:
-            plt.xlabel("RA (arb. units)")
-            plt.ylabel("Dec (arb. units)")
+        RAAxis = ax.coords['ra']
+        RAAxis.set_axislabel('RA', minpad=0.75)
+        RAAxis.set_major_formatter('hh:mm:ss')
+        DecAxis = ax.coords['dec']
+        DecAxis.set_axislabel('Dec', minpad=0.75)
+        DecAxis.set_major_formatter('dd:mm:ss')
+        ax.coords.grid(color='black', alpha=0.5, linestyle='solid')
     else:
         # Reverse the axis so that RA coord increases to left
         plt.xlim(upper[0], lower[0])
@@ -683,9 +671,10 @@ def _plot_frame(screen, fitted_phase1, residuals, weights, x, y, k, lower,
 
 
 def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
-    station_positions, source_names, times, height, station_order, beta_val,
-    r_0, prefix='frame_', remove_gradient=True, show_source_names=False,
-    min_val=None, max_val=None, is_phase=False, midRA=0.0, midDec=0.0, ncpu=0):
+                       station_positions, source_names, times, height, station_order,
+                       beta_val, r_0, prefix='frame_', remove_gradient=True,
+                       show_source_names=False, min_val=None, max_val=None, is_phase=False,
+                       midRA=0.0, midDec=0.0, ncpu=0):
     """
     Makes plots of screens
 
@@ -731,20 +720,6 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
         Number of CPUs to use
 
     """
-    from numpy import kron, concatenate, newaxis
-    from numpy.linalg import pinv, norm
-    import numpy as np
-    import os
-
-    # avoids error if re-setting "agg" a second run of plot
-    if not 'matplotlib' in sys.modules:
-        import matplotlib as mpl
-        mpl.rc('figure.subplot',left=0.05, bottom=0.05, right=0.95, top=0.95,wspace=0.22, hspace=0.22 )
-        mpl.use("Agg")
-    import matplotlib as mpl
-    import matplotlib.pyplot as plt # after setting "Agg" to speed up
-    from mpl_toolkits.axes_grid1.inset_locator import inset_axes
-
     # input check
     root_dir = os.path.dirname(prefix)
     if root_dir == '':
@@ -756,13 +731,13 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
         pass
 
     if height == 0.0:
-        N_stations = 1 # screens are single-station screens
+        N_stations = 1  # screens are single-station screens
     else:
         N_stations = len(station_names)
     N_sources = len(source_names)
     N_times = len(times)
     N_piercepoints = N_sources * N_stations
-    xp, yp, zp = station_positions[0, :] # use first station
+    xp, yp, zp = station_positions[0, :]  # use first station
     east = np.array([-yp, xp, 0])
     east = east / norm(east)
 
@@ -777,7 +752,7 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
     # Use pierce point locations of first and last time slots to estimate
     # required size of plot in meters
     if height == 0.0:
-        is_image_plane = True # pierce points are image plane coords
+        is_image_plane = True  # pierce points are image plane coords
         pp1_0 = pp[:, 0:2]
         pp1_1 = pp[:, 0:2]
     else:
@@ -799,7 +774,7 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
     im_extent_m = upper - lower
     fitted_phase1 = np.zeros((N_piercepoints, N_times))
 
-    Nx = 60 # set approximate number of pixels in screen
+    Nx = 60  # set approximate number of pixels in screen
     pix_per_m = Nx / im_extent_m[0]
     m_per_pix = 1.0 / pix_per_m
     xr = np.arange(lower[0], upper[0], m_per_pix)
@@ -809,19 +784,19 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
     lower = np.array([xr[0], yr[0]])
     upper = np.array([xr[-1], yr[-1]])
 
-    x = np.zeros((N_times, N_piercepoints)) # plot x pos of piercepoints
-    y = np.zeros((N_times, N_piercepoints)) # plot y pos of piercepoints
+    x = np.zeros((N_times, N_piercepoints))  # plot x pos of piercepoints
+    y = np.zeros((N_times, N_piercepoints))  # plot y pos of piercepoints
     screen = np.zeros((Nx, Ny, N_times))
 
     if height == 0.0:
         for sindx in range(station_positions.shape[0]):
-            logging.info('Calculating screen images...')
+            log.info('Calculating screen images...')
             residuals = inresiduals[:, :, sindx, newaxis].transpose([0, 2, 1]).reshape(N_piercepoints, N_times)
             mpm = multiprocManager(ncpu, _calculate_screen)
             for k in range(N_times):
                 mpm.put([inscreen[:, k, sindx], residuals[:, k], pp,
-                    N_piercepoints, k, east, north, up, T, Nx, Ny, sindx, height,
-                    beta_val, r_0, is_phase])
+                        N_piercepoints, k, east, north, up, T, Nx, Ny, sindx, height,
+                        beta_val, r_0, is_phase])
             mpm.wait()
             for (k, ft, scr, xa, ya) in mpm.get():
                 screen[:, :, k] = scr
@@ -830,7 +805,7 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
                     x[k, :] = xa
                     y[k, :] = ya
                 else:
-                    x[k, :] = xa - np.amin(xa) # remove offsets for each time slot
+                    x[k, :] = xa - np.amin(xa)  # remove offsets for each time slot
                     y[k, :] = ya - np.amin(ya)
 
             # Normalize piercepoint locations to extent calculated above
@@ -847,23 +822,24 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
             else:
                 vmax = max_val
 
-            logging.info('Plotting screens...')
+            log.info('Plotting screens...')
             mpm = multiprocManager(ncpu, _plot_frame)
             for k in range(N_times):
                 mpm.put([screen[:, :, k], fitted_phase1[:, k], residuals[:, k],
-                weights[:, k, sindx], x[k, :], y[k, :], k, lower, upper, vmin, vmax,
-                source_names, show_source_names, station_names, sindx, root_dir,
-                prestr, is_image_plane, midRA, midDec, station_order[0, k, sindx], is_phase])
+                        weights[:, k, sindx], x[k, :], y[k, :], k, lower, upper, vmin, vmax,
+                        source_names, show_source_names, station_names, sindx, root_dir,
+                        prestr, is_image_plane, midRA, midDec, station_order[0, k, sindx],
+                        is_phase])
             mpm.wait()
     else:
-        logging.info('Calculating screen images...')
+        log.info('Calculating screen images...')
         residuals = inresiduals.transpose([0, 2, 1]).reshape(N_piercepoints, N_times)
         weights = weights.transpose([0, 2, 1]).reshape(N_piercepoints, N_times)
         mpm = multiprocManager(ncpu, _calculate_screen)
         for k in range(N_times):
             mpm.put([inscreen[:, k, :], residuals[:, k], pp[k, :, :],
-                N_piercepoints, k, east, north, up, T, Nx, Ny, -1, height,
-                beta_val, r_0, is_phase])
+                    N_piercepoints, k, east, north, up, T, Nx, Ny, -1, height,
+                    beta_val, r_0, is_phase])
         mpm.wait()
         for (k, ft, scr, xa, ya) in mpm.get():
             screen[:, :, k] = scr
@@ -872,7 +848,7 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
                 x[k, :] = xa
                 y[k, :] = ya
             else:
-                x[k, :] = xa - np.amin(xa) # remove offsets for each time slot
+                x[k, :] = xa - np.amin(xa)  # remove offsets for each time slot
                 y[k, :] = ya - np.amin(ya)
 
         # Normalize piercepoint locations to extent calculated above
@@ -889,7 +865,7 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
         else:
             vmax = max_val
 
-        logging.info('Plotting screens...')
+        log.info('Plotting screens...')
         if show_source_names:
             pp_names = []
             for stat in station_names:
@@ -901,9 +877,9 @@ def _make_screen_plots(pp, inscreen, inresiduals, weights, station_names,
         for k in range(N_times):
             order = station_order[0]
             mpm.put([screen[:, :, k], fitted_phase1[:, k], residuals[:, k],
-            weights[:, k], x[k, :], y[k, :], k, lower, upper, vmin, vmax,
-            pp_names, show_source_names, station_names, -1, root_dir,
-            prestr, is_image_plane, midRA, midDec, order, is_phase])
+                    weights[:, k], x[k, :], y[k, :], k, lower, upper, vmin, vmax,
+                    pp_names, show_source_names, station_names, -1, root_dir,
+                    prestr, is_image_plane, midRA, midDec, order, is_phase])
         mpm.wait()
 
 
@@ -924,18 +900,17 @@ def _fitPLaneLTSQ(XYZ):
         plane parameters, where Z = aX + bY + c
 
     """
-    import numpy as np
     [rows, cols] = XYZ.shape
     G = np.ones((rows, 3))
-    G[:, 0] = XYZ[:, 0]  #X
-    G[:, 1] = XYZ[:, 1]  #Y
+    G[:, 0] = XYZ[:, 0]  # X
+    G[:, 1] = XYZ[:, 1]  # Y
     Z = XYZ[:, 2]
     (a, b, c), resid, rank, s = np.linalg.lstsq(G, Z)
     return (a, b, c)
 
 
 def runplots(soltab, resSoltab='', minZ=-3.2, maxZ=3.2, prefix='', remove_gradient=False,
-    show_source_names=False, ncpu=0):
+             show_source_names=False, ncpu=0):
     """
     Plot screens (one plot is made per time and per station)
 
@@ -959,15 +934,12 @@ def runplots(soltab, resSoltab='', minZ=-3.2, maxZ=3.2, prefix='', remove_gradie
         Number of CPUs to use. If 0, all are used.
 
     """
-    import os
-    import numpy as np
-
     screen_type = soltab.getType()
     if screen_type == 'phasescreen':
         is_phase = True
     else:
         is_phase = False
-    logging.info('Using input {0} soltab {1}'.format(screen_type, soltab.name))
+    log.info('Using input {0} soltab {1}'.format(screen_type, soltab.name))
 
     # Get values from soltabs
     solset = soltab.getSolset()
@@ -975,13 +947,12 @@ def runplots(soltab, resSoltab='', minZ=-3.2, maxZ=3.2, prefix='', remove_gradie
         try:
             # Look for residual soltab assuming standard naming conventions
             ressoltab = solset.getSoltab(soltab.name+'resid')
-        except:
-            logging.error('Could not find the soltab with associated screen residuals. '
-                'Please specify it with the "resSoltab" argument.')
-            return 1
+        except Exception:
+            raise ValueError('Could not find the soltab with associated screen residuals. '
+                             'Please specify it with the "resSoltab" argument.')
     else:
         ressoltab = solset.getSoltab(resSoltab)
-    logging.info('Using input screen residual soltab: {}'.format(ressoltab.name))
+    log.info('Using input screen residual soltab: {}'.format(ressoltab.name))
     screen = np.array(soltab.val)
     weights = np.array(soltab.weight)
     residuals = np.array(ressoltab.val)
@@ -1045,10 +1016,11 @@ def runplots(soltab, resSoltab='', minZ=-3.2, maxZ=3.2, prefix='', remove_gradie
         max_val = maxZ
 
     _make_screen_plots(pp, screen, residuals, weights, np.array(station_names),
-        np.array(station_positions), np.array(source_names), times,
-        height, orders, beta_val, r_0, prefix=prefix,
-        remove_gradient=remove_gradient, show_source_names=show_source_names, min_val=min_val,
-        max_val=max_val, is_phase=is_phase, midRA=midRA, midDec=midDec, ncpu=ncpu)
+                       np.array(station_positions), np.array(source_names), times,
+                       height, orders, beta_val, r_0, prefix=prefix,
+                       remove_gradient=remove_gradient, show_source_names=show_source_names,
+                       min_val=min_val, max_val=max_val, is_phase=is_phase, midRA=midRA,
+                       midDec=midDec, ncpu=ncpu)
 
     return 0
 
@@ -1057,7 +1029,7 @@ if __name__ == "__main__":
     descriptiontext = "Plot solutions.\n"
     parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
     parser.add_argument('h5file', help="Name of solution h5parm file")
-    parser.add_argument('soltype', help="Type of solution to plot: tec, tecerror, scalarphase, phase, or amplitude")
+    parser.add_argument('soltype', help="Type of solution to plot: phase, amplitude, phasescreen, or ampscreen")
     parser.add_argument('--root', help="Root name for output plots (default: 'soltype_')", default=None)
     parser.add_argument('--refstat', help="Name of referance station (default: first)", default=None)
     parser.add_argument('--soltab', help="Name of solution table (default: default for soltype)", default=None)
@@ -1067,5 +1039,8 @@ if __name__ == "__main__":
     parser.add_argument('--pol', help="Name of direction (default: all)", type=str, default=None)
 
     args = parser.parse_args()
-    main(args.h5file, args.soltype, args.root, args.refstat, args.soltab, args.dir, args.ant,
-         args.freq, args.pol)
+    try:
+        main(args.h5file, args.soltype, args.root, args.refstat, args.soltab, args.dir, args.ant,
+             args.freq, args.pol)
+    except Exception as e:
+        log.error('{0}.'.format(e))
diff --git a/docs/source/products.rst b/docs/source/products.rst
index e6e0c23417f92d3278f9ea57ec563262dbeee1f0..2fe40a0901a00275624adcc6d974d30d367b0383 100644
--- a/docs/source/products.rst
+++ b/docs/source/products.rst
@@ -9,7 +9,7 @@ Rapthor produces the following output inside the working directory:
     Directory containing the log files.
 
 ``images/``
-    Directory containing the images.
+    Directory containing the FITS images.
 
 ``pipelines/``
     Directory containing intermediate files of each operation's pipeline.
@@ -17,11 +17,11 @@ Rapthor produces the following output inside the working directory:
 ``regions/``
     Directory containing ds9 region files.
 
-``scratch/``
-    Directory containing temporary files.
-
 ``skymodels/``
     Directory containing sky model files.
 
 ``solutions/``
-    Directory containing the calibration solutions.
+    Directory containing the calibration solution h5parm files.
+
+``plots/``
+    Directory containing the PNG plots of the calibration solutions and images.
diff --git a/rapthor/lib/facet.py b/rapthor/lib/facet.py
index 7de7d5bbfd4100950b755a9dc6fb648139790772..1f5201e0cc3408a1a58feca5b627a15639aa9a28 100644
--- a/rapthor/lib/facet.py
+++ b/rapthor/lib/facet.py
@@ -272,11 +272,11 @@ def make_ds9_region_file(center_coords, facet_polygons, outfile, names=None):
         Decs = vertices.T[1]
         for ra, dec in zip(RAs, Decs):
             radec_list.append('{0}, {1}'.format(ra, dec))
+        lines.append('polygon({0})\n'.format(', '.join(radec_list)))
         if name is None:
-            lines.append('polygon({0})\n'.format(', '.join(radec_list)))
+            lines.append('point({0}, {1})\n'.format(center_coord[0], center_coord[1]))
         else:
-            lines.append('polygon({0}) # text={1}\n'.format(', '.join(radec_list), name))
-        lines.append('point({0}, {1})\n'.format(center_coord[0], center_coord[1]))
+            lines.append('point({0}, {1}) # text={2}\n'.format(center_coord[0], center_coord[1], name))
 
     with open(outfile, 'w') as f:
         f.writelines(lines)
diff --git a/rapthor/lib/parset.py b/rapthor/lib/parset.py
index c043460a1ec7cf1bafe351b98a4185e1ae757bbd..6052fbc37de2496cf1c1b6e62b457dbdc1500bd6 100644
--- a/rapthor/lib/parset.py
+++ b/rapthor/lib/parset.py
@@ -56,7 +56,7 @@ def parset_read(parset_file, use_log_file=True, skip_cluster=False):
     try:
         os.chdir(parset_dict['dir_working'])
         for subdir in ['logs', 'pipelines', 'regions', 'skymodels', 'images',
-                       'solutions']:
+                       'solutions', 'plots']:
             subdir_path = os.path.join(parset_dict['dir_working'], subdir)
             if not os.path.isdir(subdir_path):
                 os.mkdir(subdir_path)
diff --git a/rapthor/operations/calibrate.py b/rapthor/operations/calibrate.py
index 1d3718383cf54105e2f43b39eed8110ecf84e733..caaf013fe0d2ed52ee179d110baa7c581c1f8a1f 100644
--- a/rapthor/operations/calibrate.py
+++ b/rapthor/operations/calibrate.py
@@ -4,6 +4,7 @@ Module that holds the Calibrate class
 import os
 import logging
 import shutil
+import glob
 from rapthor.lib.operation import Operation
 from rapthor.lib import miscellaneous as misc
 from rapthor.lib.cwl import CWLFile, CWLDir
@@ -110,7 +111,6 @@ class Calibrate(Operation):
             fast_antennaconstraint = '[]'
         slow_antennaconstraint = '[[{}]]'.format(','.join(self.field.stations))
 
-
         # Get various DDECal solver parameters
         llssolver = self.field.llssolver
         maxiter = self.field.maxiter
@@ -120,15 +120,10 @@ class Calibrate(Operation):
         stepsize = self.field.stepsize
         tolerance = self.field.tolerance
         uvlambdamin = self.field.solve_min_uv_lambda
-        parallelbaselines=self.field.parallelbaselines
-        if solveralgorithm=='lbfgs':
-           solverlbfgs_dof=float(self.field.solverlbfgs_dof)
-           solverlbfgs_iter=int(self.field.solverlbfgs_iter)
-           solverlbfgs_minibatches=int(self.field.solverlbfgs_minibatches)
-        else:
-           solverlbfgs_dof=200.0
-           solverlbfgs_iter=4
-           solverlbfgs_minibatches=1
+        parallelbaselines = self.field.parallelbaselines
+        solverlbfgs_dof = self.field.solverlbfgs_dof
+        solverlbfgs_iter = self.field.solverlbfgs_iter
+        solverlbfgs_minibatches = self.field.solverlbfgs_minibatches
 
         # Get the size of the imaging area (for use in making the a-term images)
         sector_bounds_deg = '{}'.format(self.field.sector_bounds_deg)
@@ -198,10 +193,10 @@ class Calibrate(Operation):
                             'combined_slow_h5parm1': combined_slow_h5parm1,
                             'combined_slow_h5parm2': combined_slow_h5parm2,
                             'combined_h5parms1': combined_h5parms1,
-                            'combined_h5parms2': combined_h5parms2}
-        self.input_parms.update({'solverlbfgs_dof':solverlbfgs_dof,
-                 'solverlbfgs_iter':solverlbfgs_iter,
-                 'solverlbfgs_minibatches':solverlbfgs_minibatches})
+                            'combined_h5parms2': combined_h5parms2,
+                            'solverlbfgs_dof': solverlbfgs_dof,
+                            'solverlbfgs_iter': solverlbfgs_iter,
+                            'solverlbfgs_minibatches': solverlbfgs_minibatches}
 
         if self.field.debug:
             output_slow_h5parm_debug = ['slow_gain_{}_debug.h5parm'.format(i)
@@ -292,7 +287,7 @@ class Calibrate(Operation):
             self.field.aterm_image_filenames = [os.path.join(self.pipeline_working_dir, af.strip())
                                                 for af in self.field.aterm_image_filenames]
 
-        # Save the solutions
+        # Copy the solutions (h5parm files)
         dst_dir = os.path.join(self.parset['dir_working'], 'solutions', 'calibrate_{}'.format(self.index))
         misc.create_directory(dst_dir)
         self.field.h5parm_filename = os.path.join(dst_dir, 'field-solutions.h5')
@@ -304,3 +299,13 @@ class Calibrate(Operation):
         else:
             shutil.copy(os.path.join(self.pipeline_working_dir, self.combined_fast_h5parm),
                         os.path.join(dst_dir, self.field.h5parm_filename))
+
+        # Copy the plots (PNG files)
+        dst_dir = os.path.join(self.parset['dir_working'], 'plots', 'calibrate_{}'.format(self.index))
+        misc.create_directory(dst_dir)
+        plot_filenames = glob.glob(os.path.join(self.pipeline_working_dir, '*.png'))
+        for plot_filename in plot_filenames:
+            dst_filename = os.path.join(dst_dir, os.path.basename(plot_filename))
+            if os.path.exists(dst_filename):
+                os.remove(dst_filename)
+            shutil.copy(plot_filename, dst_filename)
diff --git a/rapthor/operations/image.py b/rapthor/operations/image.py
index 2f3e8d38f787f3ceedd98b72a2c532f464c33ac6..b38f6cf3026d7c84eee076158e79cdc676ec26e8 100644
--- a/rapthor/operations/image.py
+++ b/rapthor/operations/image.py
@@ -3,6 +3,8 @@ Module that holds the Image class
 """
 import os
 import logging
+import shutil
+from rapthor.lib import miscellaneous as misc
 from rapthor.lib.operation import Operation
 from rapthor.lib.cwl import CWLFile, CWLDir
 
@@ -207,7 +209,7 @@ class Image(Operation):
         """
         Finalize this operation
         """
-        # Save output FITS image and model for each sector
+        # Save output FITS image, sky model, and ds9 region file for each sector
         # NOTE: currently, -save-source-list only works with pol=I -- when it works with other
         # pols, save them all
         for sector in self.field.imaging_sectors:
@@ -230,6 +232,17 @@ class Image(Operation):
             sector.image_skymodel_file_true_sky = image_root + '.true_sky.txt'
             sector.image_skymodel_file_apparent_sky = image_root + '.apparent_sky.txt'
 
+            # The ds9 region file, if made
+            if self.field.dde_method == 'facets':
+                dst_dir = os.path.join(self.parset['dir_working'], 'regions', 'image_{}'.format(self.index))
+                misc.create_directory(dst_dir)
+                region_filename = '{}_facets_ds9.reg'.format(sector.name)
+                src_filename = os.path.join(self.pipeline_working_dir, region_filename)
+                dst_filename = os.path.join(dst_dir, region_filename)
+                if os.path.exists(dst_filename):
+                    os.remove(dst_filename)
+                shutil.copy(src_filename, dst_filename)
+
         # Symlink to datasets and remove old ones
 #         dst_dir = os.path.join(self.parset['dir_working'], 'datasets', self.direction.name)
 #         misc.create_directory(dst_dir)
diff --git a/rapthor/pipeline/parsets/calibrate_pipeline.cwl b/rapthor/pipeline/parsets/calibrate_pipeline.cwl
index 182a7b0e33f1394aa1903ac8b125b637391be389..3764db47fd25676ac6ded1d18582d1be8305a3ed 100644
--- a/rapthor/pipeline/parsets/calibrate_pipeline.cwl
+++ b/rapthor/pipeline/parsets/calibrate_pipeline.cwl
@@ -10,7 +10,7 @@ doc: |
   further unconstrained slow gain calibration to correct for station-to-station
   differences. Steps (2) and (3) are skipped if the calibration is phase-only.
   This calibration scheme works for both HBA and LBA data. The final products of
-  this pipeline are solution tables (h5parm files) and a-term screens (FITS
+  this pipeline are solution tables (h5parm files), plots, and a-term screens (FITS
   files).
 
 requirements:
@@ -369,12 +369,26 @@ outputs:
     outputSource:
       - adjust_h5parm_sources/adjustedh5parm
     type: File
+  - id: fast_phase_plots
+    outputSource:
+      - plot_fast_phase_solutions/plots
+    type: File[]
 {% if use_screens %}
   - id: diagonal_aterms
     outputSource:
       - merge_aterm_files/output
     type: File[]
 {% endif %}
+{% if do_slowgain_solve %}
+  - id: slow_phase_plots
+    outputSource:
+      - plot_slow_phase_solutions/plots
+    type: File[]
+  - id: slow_amp_plots
+    outputSource:
+      - plot_slow_amp_solutions/plots
+    type: File[]
+{% endif %}
 
 
 steps:
@@ -465,6 +479,19 @@ steps:
     out:
       - id: outh5parm
 
+  - id: plot_fast_phase_solutions
+    label: Plot fast phase solutions
+    doc: |
+      This step makes plots of the fast phase solutions.
+    run: {{ rapthor_pipeline_dir }}/steps/plot_solutions.cwl
+    in:
+      - id: h5parm
+        source: combine_fast_phases/outh5parm
+      - id: soltype
+        valueFrom: 'phase'
+    out:
+      - id: plots
+
 {% if do_slowgain_solve %}
 # start do_slowgain_solve
 
@@ -728,6 +755,32 @@ steps:
     out:
       - id: outh5parm
 
+  - id: plot_slow_phase_solutions
+    label: Plot slow phase solutions
+    doc: |
+      This step makes plots of the slow phase solutions.
+    run: {{ rapthor_pipeline_dir }}/steps/plot_solutions.cwl
+    in:
+      - id: h5parm
+        source: normalize_slow_amplitudes/outh5parm
+      - id: soltype
+        valueFrom: 'phase'
+    out:
+      - id: plots
+
+  - id: plot_slow_amp_solutions
+    label: Plot slow amp solutions
+    doc: |
+      This step makes plots of the slow amplitude solutions.
+    run: {{ rapthor_pipeline_dir }}/steps/plot_solutions.cwl
+    in:
+      - id: h5parm
+        source: normalize_slow_amplitudes/outh5parm
+      - id: soltype
+        valueFrom: 'amplitude'
+    out:
+      - id: plots
+
   - id: combine_fast_and_slow_h5parms2
     label: Combine fast-phase and slow-gain solutions
     doc: |
diff --git a/rapthor/pipeline/parsets/image_pipeline.cwl b/rapthor/pipeline/parsets/image_pipeline.cwl
index 9d36954b32b8ddcfb584efe16a10297851ef1ef7..2dddb4d9eabb3eb65d93edde64ec491626c02894 100644
--- a/rapthor/pipeline/parsets/image_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_pipeline.cwl
@@ -354,6 +354,15 @@ outputs:
       items:
         type: array
         items: File
+{% if use_facets %}
+  - id: region_file
+    outputSource:
+      - image_sector/region_file
+    type:
+      type: array
+      items: File
+{% endif %}
+
 
 steps:
   - id: image_sector
@@ -512,3 +521,4 @@ steps:
     out:
       - id: filtered_skymodels
       - id: sector_images
+      - id: region_file
diff --git a/rapthor/pipeline/parsets/image_sector_pipeline.cwl b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
index 95adff176eae2ae615b014c142e7cb5fecd16ad5..cf189ac8774271d4dac7a0c9ccb648f865387d4b 100644
--- a/rapthor/pipeline/parsets/image_sector_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
@@ -319,6 +319,12 @@ outputs:
       - image/skymodel_nonpb
       - image/skymodel_pb
     type: File[]
+{% if use_facets %}
+  - id: region_file
+    outputSource:
+      - make_region_file/region_file
+    type: File
+{% endif %}
 
 steps:
   - id: prepare_imaging_data
diff --git a/rapthor/pipeline/steps/make_region_file.cwl b/rapthor/pipeline/steps/make_region_file.cwl
index 486ce856febd1ca43e7ace91c8d283d515e3a074..0cad4d6d4cfa3a705af429421d6486e157805bb5 100644
--- a/rapthor/pipeline/steps/make_region_file.cwl
+++ b/rapthor/pipeline/steps/make_region_file.cwl
@@ -61,7 +61,6 @@ outputs:
       parameter "outfile".
     type:
       - File
-      - string
     outputBinding:
       glob: $(inputs.outfile)
 
diff --git a/rapthor/pipeline/steps/plot_solutions.cwl b/rapthor/pipeline/steps/plot_solutions.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..c5184d496f45a1d6365febd146afd2cb10ed9683
--- /dev/null
+++ b/rapthor/pipeline/steps/plot_solutions.cwl
@@ -0,0 +1,35 @@
+cwlVersion: v1.2
+class: CommandLineTool
+baseCommand: [plotrapthor]
+label: Plots calibration solutions
+doc: |
+  This tool plots the calibration solutions.
+
+inputs:
+  - id: h5parm
+    label: Input solution table
+    doc: |
+      The filename of the input h5parm file.
+    type: File
+    inputBinding:
+      position: 0
+  - id: soltype
+    label: Solution type
+    doc: |
+      The type of solution to plot (one of phase, amplitude, phasescreen, or ampscreen).
+    type: string
+    inputBinding:
+      position: 1
+
+outputs:
+  - id: plots
+    label: Output plots
+    doc: |
+      The filenames of the output plots.
+    type: File[]
+    outputBinding:
+      glob: '*.png'
+
+hints:
+  - class: DockerRequirement
+    dockerPull: 'astronrd/rapthor'
diff --git a/rapthor/scripts/calculate_bl_lengths.py b/rapthor/scripts/calculate_bl_lengths.py
deleted file mode 100755
index 728d1408f3de6387ce1e29744bfb70e940ebe398..0000000000000000000000000000000000000000
--- a/rapthor/scripts/calculate_bl_lengths.py
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/usr/bin/env python3
-"""
-Calculates the baseline lengths in km for a MS
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import numpy as np
-import sys
-import os
-import itertools
-import casacore.tables as pt
-import pickle
-
-
-def main(ms_input, output_file):
-    """
-    Calculates the baseline lengths in km for a MS
-
-    Parameters
-    ----------
-    ms_input : list or str
-        List of MS filenames, or string with list, or path to a mapfile. The MS
-        files should all have the same frequency
-    output_file : str
-        Filename of output pickle file
-
-    """
-    ms_list = input2strlist(ms_input)
-
-    print('Calculating baseline lengths...')
-    baseline_dict = get_baseline_lengths(ms_list[0])
-
-    with open(output_file, 'wb') as f:
-        pickle.dump(baseline_dict, f)
-
-
-def input2strlist(invar):
-    str_list = None
-    if type(invar) is str:
-        if invar.startswith('[') and invar.endswith(']'):
-            str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
-    elif type(invar) is list:
-        str_list = [str(f).strip(' \'\"') for f in invar]
-    else:
-        raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
-    return str_list
-
-
-def get_baseline_lengths(ms):
-    """
-    Returns dict of baseline lengths in km for all baselines in input dataset
-    """
-    anttab = pt.table(ms+'::ANTENNA', ack=False)
-    antnames = anttab.getcol('NAME')
-    anttab.close()
-
-    t = pt.table(ms, ack=False)
-    ant1 = t.getcol('ANTENNA1')
-    ant2 = t.getcol('ANTENNA2')
-    all_uvw = t.getcol('UVW')
-    t.close()
-
-    baseline_dict = {}
-    for ant in itertools.product(set(ant1), set(ant2)):
-        if ant[0] >= ant[1]:
-            continue
-        sel1 = np.where(ant1 == ant[0])[0]
-        sel2 = np.where(ant2 == ant[1])[0]
-        sel = sorted(list(frozenset(sel1).intersection(sel2)))
-        uvw = all_uvw[sel, :]
-        uvw_dist = np.sqrt(uvw[:, 0]**2 + uvw[:, 1]**2 + uvw[:, 2]**2)
-        baseline_dict['{0}'.format(ant[0])] = antnames[ant[0]]
-        baseline_dict['{0}'.format(ant[1])] = antnames[ant[1]]
-        baseline_dict['{0}-{1}'.format(ant[0], ant[1])] = np.mean(uvw_dist) / 1.e3
-
-    return baseline_dict
diff --git a/rapthor/scripts/compare_image_stats.py b/rapthor/scripts/compare_image_stats.py
index 88db51ebf9f98d8c8da4ecc8af9671f39ccf5ad1..950aa446df03b74018cd1cc9dd03457297232d31 100755
--- a/rapthor/scripts/compare_image_stats.py
+++ b/rapthor/scripts/compare_image_stats.py
@@ -194,15 +194,15 @@ def find_imagenoise(imagename):
     return rms, numpy.abs(numpy.max(image)/rms), minmax
 
 
-def main(im1, im2, count=-1, rapthor=1.0125):
+def main(im1, im2, count=-1, factor=1.0125):
     """
     Compare the dynamic range and min/max of two images and check whether:
 
-        dynamic_range1 / rapthor > dynamic_range2
+        dynamic_range1 / factor > dynamic_range2
 
     or
 
-        abs(min1/max1) * rapthor < abs(min2/max2)
+        abs(min1/max1) * factor < abs(min2/max2)
 
     Typically, im1 is the latest image and im2 the previous one.
 
@@ -215,17 +215,17 @@ def main(im1, im2, count=-1, rapthor=1.0125):
     count : int, optional
         Loop index. If nonzero, the dynamic range check is skipped for count = 0
         only and break is set to False
-    rapthor : float
-        Required improvement rapthor for success (i.e., break = True)
+    factor : float
+        Required improvement factor for success (i.e., break = True)
 
     Returns
     -------
     result : dict
-        Dict with break value (False if dynamic_range1 / rapthor > dynamic_range2;
+        Dict with break value (False if dynamic_range1 / factor > dynamic_range2;
         True otherwise)
 
     """
-    rapthor = float(rapthor)
+    factor = float(factor)
     count = int(count)
 
     rms1, dynamic_range1, minmax1 =  find_imagenoise(im1)
@@ -242,7 +242,7 @@ def main(im1, im2, count=-1, rapthor=1.0125):
     else:
         # Check whether dynamic range is increasing or minmax is decreasing. If
         # so, continue (return False)
-        if (dynamic_range1 / rapthor > dynamic_range2) or (minmax1 * rapthor < minmax2):
+        if (dynamic_range1 / factor > dynamic_range2) or (minmax1 * factor < minmax2):
             return {'break': False}
         else:
             return {'break': True}
@@ -254,7 +254,7 @@ if __name__ == '__main__':
     parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
     parser.add_argument('im1', help='name of image #1')
     parser.add_argument('im2', help='name of image #2')
-    parser.add_argument('rapthor', help='required improvement rapthor for success')
+    parser.add_argument('factor', help='required improvement factor for success')
     args = parser.parse_args()
 
-    main(args.im1, args.im2, rapthor=args.rapthor)
+    main(args.im1, args.im2, factor=args.factor)
diff --git a/rapthor/scripts/flag_residual_gains.py b/rapthor/scripts/flag_residual_gains.py
deleted file mode 100755
index cbc5076fa264e9a509baa2c4866718ab5c26ce8a..0000000000000000000000000000000000000000
--- a/rapthor/scripts/flag_residual_gains.py
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to flag on residual gain solutions
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import losoto.operations as operations
-from losoto.h5parm import h5parm
-
-
-def main(h5parmfile, solsetname='sol000', phasesoltabname='phase000', ampsoltabname='amplitude000'):
-    """
-    Flag residual gains
-
-    Parameters
-    ----------
-    h5parmfile : str
-        Filename of h5parm
-    phasesoltabname : str, optional
-        Name of phase soltab
-    ampsoltabname : str, optional
-        Name of amplitude soltab
-    """
-    # Read in solutions
-    H = h5parm(h5parmfile, readonly=False)
-    solset = H.getSolset(solsetname)
-    phsoltab = solset.getSoltab(phasesoltabname)
-    ampsoltab = solset.getSoltab(ampsoltabname)
-
-    # Flag
-    operations.flagstation.run(phsoltab, 'resid', nSigma=3.0)
-    operations.flagstation.run(ampsoltab, 'resid', nSigma=5.0)
-
-    # Reset phases to 0 and amps to 1, so only the new flags are applied
-    operations.reset.run(phsoltab)
-    operations.reset.run(ampsoltab)
diff --git a/rapthor/scripts/pre_average_freq.py b/rapthor/scripts/pre_average_freq.py
deleted file mode 100755
index 7984898d172e7fbd630879d4cebeb1348f00fa67..0000000000000000000000000000000000000000
--- a/rapthor/scripts/pre_average_freq.py
+++ /dev/null
@@ -1,177 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to pre-average data using a sliding Gaussian kernel in frequency
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import numpy as np
-import sys
-import os
-from scipy.ndimage.filters import gaussian_filter1d as gfilter
-from scipy.special import erf
-import casacore.tables as pt
-import pickle
-import itertools
-
-
-def main(ms_input, input_colname, output_data_colname, output_weights_colname,
-    baseline_file, delta_theta_deg, target_peak_reduction_rapthor=0.99):
-    """
-    Pre-average data using a sliding Gaussian kernel in frequency
-
-    Parameters
-    ----------
-    ms_input : str
-        MS filename
-    input_colname : str
-        Name of the column in the MS from which the data are read
-    output_data_colname : str
-        Name of the column in the MS into which the averaged data are written
-    output_weights_colname : str
-        Name of the column in the MS into which the averaged data weights are
-        written
-    baseline_file : str
-        Filename of pickled baseline lengths
-    delta_theta_deg : float
-        Radius of calibration region in degrees
-    target_peak_reduction_rapthor : float, optional
-        Target reduction in peak flux density. Note: this reduction is in
-        addition to any incurred by earlier averaging
-
-    """
-    if os.path.exists(baseline_file):
-        f = open(baseline_file, 'r')
-        baseline_dict = pickle.load(f)
-        f.close()
-    else:
-        print('Cannot find baseline_file. Exiting...')
-        sys.exit(1)
-    delta_theta_deg = float(delta_theta_deg)
-    target_peak_reduction_rapthor = float(target_peak_reduction_rapthor)
-
-    ms = pt.table(ms_input, readonly=False, ack=False)
-    ant1_list = ms.getcol('ANTENNA1')
-    ant2_list = ms.getcol('ANTENNA2')
-    data_all = ms.getcol(input_colname)
-    weights_all = ms.getcol('WEIGHT_SPECTRUM')
-    flags = ms.getcol('FLAG')
-
-    # Get lowest frequency of MS and channel width
-    sw = pt.table(ms_input+'::SPECTRAL_WINDOW', ack=False)
-    freq_hz = sw.col('CHAN_FREQ')[0][0]
-    chan_width_hz = sw.col('CHAN_WIDTH')[0][0]
-
-    flags[ np.isnan(data_all) ] = True # flag NaNs
-    weights_all = weights_all * ~flags # set weight of flagged data to 0
-
-    # Check that all NaNs are flagged
-    if np.count_nonzero(np.isnan(data_all[~flags])) > 0:
-        logging.error('NaNs in unflagged data in {0}!'.format(ms_input))
-        sys.exit(1)
-
-    # Weight data and set bad data to 0 so nans do not propagate
-    data_all = np.nan_to_num(data_all*weights_all)
-
-    # Iteration on baseline combination
-    for ant in itertools.product(set(ant1_list), set(ant2_list)):
-        if ant[0] >= ant[1]:
-            continue
-        sel1 = np.where(ant1_list == ant[0])[0]
-        sel2 = np.where(ant2_list == ant[1])[0]
-        sel_list = sorted(list(frozenset(sel1).intersection(sel2)))
-
-        data = data_all[sel_list,:,:]
-        weights = weights_all[sel_list,:,:]
-
-        # compute the Gaussian sigma from the max bandwidth over which we
-        # can average and avoid significant bandwidth smearing but limited to
-        # no more than 3 MHz (to avoid smoothing over the beam-induced effects)
-        lambda_km = 299792.458 / freq_hz
-        dist_km = baseline_dict['{0}-{1}'.format(ant[0], ant[1])]
-        resolution_deg = lambda_km / dist_km * 180.0 / np.pi
-        stddev_hz = min(3e6, get_target_bandwidth(freq_hz, delta_theta_deg,
-            resolution_deg, target_peak_reduction_rapthor)/4.0)
-        stddev_nchan = stddev_hz / chan_width_hz * np.sqrt(0.5 / dist_km)
-
-        # smear weighted data and weights
-        dataR = gfilter(np.real(data), stddev_nchan, axis=1)
-        dataI = gfilter(np.imag(data), stddev_nchan, axis=1)
-        weights = gfilter(weights, stddev_nchan, axis=1)
-
-        # re-create data
-        data = (dataR + 1j * dataI)
-        data[(weights != 0)] /= weights[(weights != 0)] # avoid divbyzero
-        data_all[sel_list,:,:] = data
-        weights_all[sel_list,:,:] = weights
-
-    # Add the output columns if needed
-    if output_data_colname not in ms.colnames():
-        desc = ms.getcoldesc(input_colname)
-        desc['name'] = output_data_colname
-        ms.addcols(desc)
-    if output_weights_colname not in ms.colnames():
-        desc = ms.getcoldesc('WEIGHT_SPECTRUM')
-        desc['name'] = output_weights_colname
-        ms.addcols(desc)
-
-    ms.putcol(output_data_colname, data_all)
-    ms.putcol('FLAG', flags) # this saves flags of nans, which is always good
-    ms.putcol(output_weights_colname, weights_all)
-    ms.close()
-
-
-def get_bandwidth_smearing_rapthor(freq, delta_freq, delta_theta, resolution):
-    """
-    Returns peak flux density reduction rapthor due to bandwidth smearing
-
-    Parameters
-    ----------
-    freq : float
-        Frequency at which averaging will be done
-    delta_freq : float
-        Bandwidth over which averaging will be done
-    delta_theta : float
-        Distance from phase center
-    resolution : float
-        Resolution of restoring beam
-
-    Returns
-    -------
-    reduction_rapthor : float
-        Ratio of post-to-pre averaging peak flux density
-
-    """
-    beta = (delta_freq/freq) * (delta_theta/resolution)
-    gamma = 2*(np.log(2)**0.5)
-    reduction_rapthor = ((np.pi**0.5)/(gamma * beta)) * (erf(beta*gamma/2.0))
-
-    return reduction_rapthor
-
-
-def get_target_bandwidth(freq, delta_theta, resolution, reduction_rapthor):
-    """
-    Returns the bandwidth for given peak flux density reduction rapthor
-
-    Parameters
-    ----------
-    freq : float
-        Frequency at which averaging will be done
-    delta_theta : float
-        Distance from phase center
-    resolution : float
-        Resolution of restoring beam
-    reduction_rapthor : float
-        Ratio of post-to-pre averaging peak flux density
-
-    Returns
-    -------
-    delta_freq : float
-        Bandwidth over which averaging will be done
-    """
-    # Increase delta_freq until we drop below target reduction_rapthor
-    delta_freq = 1e-3 * freq
-    while get_bandwidth_smearing_rapthor(freq, delta_freq, delta_theta,
-        resolution) > reduction_rapthor:
-        delta_freq *= 1.1
-
-    return delta_freq
diff --git a/rapthor/scripts/pre_average_time.py b/rapthor/scripts/pre_average_time.py
deleted file mode 100755
index f9c60579b6f13a203402ce8f94c125077f07036c..0000000000000000000000000000000000000000
--- a/rapthor/scripts/pre_average_time.py
+++ /dev/null
@@ -1,428 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to pre-average data using a sliding Gaussian kernel in time
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import numpy as np
-import glob
-import sys
-import os
-import itertools
-import pickle
-from scipy.ndimage.filters import gaussian_filter1d as gfilter
-import casacore.tables as pt
-import lofar.parmdb
-from astropy.stats import median_absolute_deviation
-
-
-def main(ms_input, input_colname, output_data_colname, output_weights_colname,
-    baseline_file, parmdb_input=None, target_rms_rad=None, ionrapthor=None,
-    minutes_per_block=10.0, verbose=True):
-    """
-    Pre-average data using a sliding Gaussian kernel in time
-
-    Parameters
-    ----------
-    ms_input : list or str
-        List of MS filenames, or string with list, or path to a mapfile
-    parmdb_input : list or str
-        List of parmdb filenames, or string with list, or path to a mapfile
-        The resulting list from parmdb_input must match the one from ms_input
-    input_colname : str
-        Name of the column in the MS from which the data are read
-    output_data_colname : str
-        Name of the column in the MS into which the averaged data are written
-    output_weights_colname : str
-        Name of the column in the MS into which the averaged data weights are written
-    target_rms_rad : float (str)
-        The target RMS for the phase noise in the input parmDBs. (Or whatever???)
-    baseline_file : str
-        Filename of pickled baseline lengths
-
-    """
-
-    # convert input to needed types
-    verbose = input2bool(verbose)
-    ms_list = input2strlist(ms_input)
-    if parmdb_input is not None:
-        parmdb_list = input2strlist(parmdb_input)
-        if len(ms_list) != len(parmdb_list):
-            raise ValueError('pre_average_time: Length of MS-list ({0}) and length of parmdb-list ({1}) differ.'.format(len(ms_list),len(parmdb_list)))
-    if type(target_rms_rad) is str:
-        target_rms_rad = float(target_rms_rad)
-    if type(ionrapthor) is str:
-        ionrapthor = float(ionrapthor)
-    if os.path.exists(baseline_file):
-        f = open(baseline_file, 'r')
-        baseline_dict = pickle.load(f)
-        f.close()
-    else:
-        print('Cannot find baseline_file. Exiting...')
-        sys.exit(1)
-
-    # Iterate through time chunks and find the lowest ionrapthor
-    start_times = []
-    end_times = []
-    ionrapthors = []
-    if verbose:
-        print('Determining ionrapthors...')
-    for msind in range(len(ms_list)):
-        tab = pt.table(ms_list[msind], ack=False)
-        start_time = tab[0]['TIME']
-        end_time = tab[-1]['TIME']
-        remaining_time = end_time - start_time # seconds
-        start_times.append(start_time)
-        end_times.append(end_time)
-        tab.close()
-        t_delta = minutes_per_block * 60.0 # seconds
-        t1 = 0.0
-        while remaining_time > 0.0:
-            if remaining_time < 1.5 * t_delta:
-                # If remaining time is too short, just include it all in this chunk
-                t_delta = remaining_time + 10.0
-            remaining_time -= t_delta
-
-            # Find ionrapthor for this period
-            if parmdb_input is not None:
-                ionrapthors.append(find_ionrapthor(parmdb_list[msind], baseline_dict, t1+start_time,
-                                                 t1+start_time+t_delta, target_rms_rad=target_rms_rad))
-                if verbose:
-                    print('    ionrapthor (for timerange {0}-{1} sec) = {2}'.format(t1,
-                          t1+t_delta, ionrapthors[-1]))
-                t1 += t_delta
-            else:
-                ionrapthors = [ionrapthor]
-
-    sorted_ms_tuples = sorted(zip(start_times,end_times, list(range(len(ms_list))), ms_list))
-    sorted_ms_dict = { 'msnames' :[ms for starttime,endtime,index,ms in sorted_ms_tuples],
-                       'starttimes' : [starttime for starttime,endtime,index,ms in sorted_ms_tuples],
-                       'endtimes' : [endtime for starttime,endtime,index,ms in sorted_ms_tuples] }
-
-    # Do pre-averaging using lowest ionrapthor
-    ionrapthor_min = min(ionrapthors)
-    if verbose:
-        print('Using ionrapthor = {}'.format(ionrapthor_min))
-        print('Averaging...')
-    BLavg_multi(sorted_ms_dict, baseline_dict, input_colname, output_data_colname,
-        output_weights_colname, ionrapthor_min)
-
-
-def find_ionrapthor(parmdb_file, baseline_dict, t1, t2, target_rms_rad=0.2):
-    """
-    Finds ionospheric scaling rapthor
-    """
-    pdb_in = lofar.parmdb.parmdb(parmdb_file)
-    parms = pdb_in.getValuesGrid('*')
-
-    # Filter any stations not in both the instrument table and the ms
-    stations_pbd = set([s.split(':')[-1] for s in pdb_in.getNames()])
-    stations_ms = set([s for s in baseline_dict.values() if type(s) is str])
-    stations = sorted(list(stations_pbd.intersection(stations_ms)))
-
-    # Select long baselines only (BL > 10 km), as they will set the ionrapthor scaling
-    ant1 = []
-    ant2 = []
-    dist = []
-    min_length = 10.0
-    for k, v in baseline_dict.items():
-        if type(v) is not str and '-' in k:
-            if v > min_length:
-                s1 = k.split('-')[0]
-                s2 = k.split('-')[1]
-                s1_name = baseline_dict[s1]
-                s2_name = baseline_dict[s2]
-                if s1_name in stations and s2_name in stations:
-                    ant1.append(s1_name)
-                    ant2.append(s2_name)
-                    dist.append(v)
-
-    # Find correlation times
-    rmstimes = []
-    dists = []
-    freq = None
-    for a1, a2, d in zip(ant1, ant2, dist):
-        if freq is None:
-            freq = np.copy(parms['Gain:0:0:Phase:{}'.format(a1)]['freqs'])[0]
-            times = np.copy(parms['Gain:0:0:Phase:{}'.format(a1)]['times'])
-            time_ind = np.where((times >= t1) & (times < t2))[0]
-            timepersolution = np.copy(parms['Gain:0:0:Phase:{}'.format(a1)]['timewidths'])[0]
-        ph1 = np.copy(parms['Gain:0:0:Phase:{}'.format(a1)]['values'])[time_ind]
-        ph2 = np.copy(parms['Gain:0:0:Phase:{}'.format(a2)]['values'])[time_ind]
-
-        # Filter flagged solutions
-        good = np.where((~np.isnan(ph1)) & (~np.isnan(ph2)))[0]
-        if len(good) == 0:
-            continue
-
-        rmstime = None
-        ph = unwrap_fft(ph2[good] - ph1[good])
-
-        step = 1
-        for i in range(1, len(ph)/2, step):
-            p1 = ph[i:]
-            p2 = ph[:-i]
-            mean = np.mean(p1-p2)
-            if mean > target_rms_rad:
-                rmstime = i
-                break
-        if rmstime is None:
-            rmstime = len(ph)/2
-        rmstimes.append(rmstime)
-        dists.append(d)
-
-    # Find the mean ionrapthor assuming that the correlation time goes as
-    # t_corr ~ 1/sqrt(BL). The ionrapthor is defined in BLavg() as:
-    #
-    #     ionrapthor = (t_corr / 30.0 sec) / ( np.sqrt((25.0 / dist_km)) * (freq_hz / 60.e6) )
-    #
-    ionrapthor = np.mean(np.array(rmstimes) / 30.0 / (np.sqrt(25.0 / np.array(dists))
-        * freq / 60.0e6)) * timepersolution
-
-    return ionrapthor
-
-
-def BLavg_multi(sorted_ms_dict, baseline_dict, input_colname, output_data_colname,
-        output_weights_colname, ionrapthor, clobber=True, maxgap_sec=1800,
-        check_files = True):
-    """
-    Averages data using a sliding Gaussian kernel on the weights
-    """
-
-    #### sort msnames into groups with gaps < maxgap_sec
-    nfiles = len(sorted_ms_dict['msnames'])
-    ms_groups = []
-    newgroup = []
-    for msindex in range(nfiles):
-        if msindex+1 == nfiles or sorted_ms_dict['starttimes'][msindex+1] > sorted_ms_dict['endtimes'][msindex] + maxgap_sec:
-            newgroup.append(sorted_ms_dict['msnames'][msindex])
-            ms_groups.append(newgroup)
-            newgroup = []
-        else:
-            newgroup.append(sorted_ms_dict['msnames'][msindex])
-
-    print("BLavg_multi: Working on",len(ms_groups),"groups of measurement sets.")
-    #### loop over all groups
-    msindex = 0
-    for ms_names in ms_groups:
-        ### collect data from all files in this group
-        freqtab = pt.table(ms_names[0] + '::SPECTRAL_WINDOW', ack=False)
-        freq = freqtab.getcell('REF_FREQUENCY',0)
-        freqtab.close()
-        timepersample = None
-        ant1_list        = []
-        ant2_list        = []
-        all_time_list    = []
-        all_data_list    = []
-        all_weights_list = []
-        all_flags_list   = []
-        for msfile in ms_names:
-            if not os.path.exists(msfile):
-                print("Cannot find MS file: {0}.".format(msfile))
-                sys.exit(1)
-            # open input/output MS
-            ms = pt.table(msfile, readonly=True, ack=False)
-            if check_files:
-                freqtab = pt.table(msfile + '::SPECTRAL_WINDOW', ack=False)
-                if freqtab.getcell('REF_FREQUENCY',0) != freq:
-                    print("Different REF_FREQUENCYs: {0} and: {1} in {2}.".format(freq,freqtab.getcell('REF_FREQUENCY',0),msfile))
-                    sys.exit(1)
-                freqtab.close()
-            #wav = 299792458. / freq
-            if timepersample is None:
-                timepersample = ms.getcell('INTERVAL',0)
-            elif check_files:
-                if timepersample != ms.getcell('INTERVAL',0):
-                    print("Different INTERVALs: {0} and: {1} in {2}.".format(timepersample,ms.getcell('INTERVAL',0),msfile))
-                    sys.exit(1)
-            all_time_list.append( ms.getcol('TIME_CENTROID') )
-            ant1_list.append( ms.getcol('ANTENNA1') )
-            ant2_list.append( ms.getcol('ANTENNA2') )
-            all_data_list.append( ms.getcol(input_colname) )
-            all_weights_list.append( ms.getcol('WEIGHT_SPECTRUM') )
-            all_flags_list.append( ms.getcol('FLAG') )
-
-            all_flags_list[-1][ np.isnan(all_data_list[-1]) ] = True # flag NaNs
-            all_weights_list[-1] = all_weights_list[-1] * ~all_flags_list[-1] # set weight of flagged data to 0
-
-            # Check that all NaNs are flagged
-            if np.count_nonzero(np.isnan(all_data_list[-1][~all_flags_list[-1]])) > 0:
-                logging.error('NaNs in unflagged data in {0}!'.format(msfile))
-                sys.exit(1)
-
-        ### iteration on baseline combination
-        for ant in itertools.product(set(ant1_list[0]), set(ant2_list[0])):
-            if ant[0] >= ant[1]:
-                continue
-            sel_list = []
-            weights_list = []
-            data_list = []
-            # select data from all MSs
-            for msindex in range(len(ms_names)):
-                sel1 = np.where(ant1_list[msindex] == ant[0])[0]
-                sel2 = np.where(ant2_list[msindex] == ant[1])[0]
-                sel_list.append( sorted(list(frozenset(sel1).intersection(sel2))) )
-
-                # # get weights and data
-                # weights_list.append( all_weights[sel_list[msindex],:,:] )
-                # data_list.append( all_data[sel_list[msindex],:,:] )
-            # combine data and weights into one array
-            data = all_data_list[0][sel_list[0],:,:]
-            weights = all_weights_list[0][sel_list[0],:,:]
-            fillshape = list(data.shape)
-            startidx = [0]
-            endidx = [data.shape[0]]
-            for msindex in range(1,len(ms_names)):
-                #pad gap between obs
-                filltimes = np.arange(np.max(all_time_list[msindex-1]),np.min(all_time_list[msindex]),timepersample)
-                fillshape[0] = len(filltimes)
-                data = np.concatenate( (data,np.zeros(fillshape)), axis=0 )
-                weights = np.concatenate( (weights,np.zeros(fillshape)), axis=0  )
-                startidx.append(data.shape[0])
-                data = np.concatenate( (data,all_data_list[msindex][sel_list[msindex],:,:]), axis=0  )
-                weights = np.concatenate( (weights,all_weights_list[msindex][sel_list[msindex],:,:]), axis=0  )
-                endidx.append(data.shape[0])
-
-            # compute the FWHM
-            dist = baseline_dict['{0}-{1}'.format(ant[0], ant[1])]
-            stddev = 30.0 * ionrapthor * np.sqrt((25.0 / dist)) * (freq / 60.e6) # in sec
-            stddev = stddev/timepersample # in samples
-
-            # Multiply every element of the data by the weights, convolve both
-            # the scaled data and the weights, and then divide the convolved data
-            # by the convolved weights (translating flagged data into weight=0).
-            # That's basically the equivalent of a running weighted average with
-            # a Gaussian window function.
-
-            # weigth data and set bad data to 0 so nans do not propagate
-            data = np.nan_to_num(data*weights)
-
-            # smear weighted data and weights
-            dataR = gfilter(np.real(data), stddev, axis=0)
-            dataI = gfilter(np.imag(data), stddev, axis=0)
-            weights = gfilter(weights, stddev, axis=0)
-
-            # re-create data
-            data = (dataR + 1j * dataI)
-            data[(weights != 0)] /= weights[(weights != 0)] # avoid divbyzero
-            for msindex in range(len(ms_names)):
-                all_data_list[msindex][sel_list[msindex],:,:] = data[startidx[msindex]:endidx[msindex],:,:]
-                all_weights_list[msindex][sel_list[msindex],:,:] = weights[startidx[msindex]:endidx[msindex],:,:]
-
-        ### write the data back to the files
-        for msindex in range(len(ms_names)):
-            ms = pt.table(ms_names[msindex], readonly=False, ack=False)
-            # Add the output columns if needed
-            if output_data_colname not in ms.colnames():
-                desc = ms.getcoldesc(input_colname)
-                desc['name'] = output_data_colname
-                ms.addcols(desc)
-            if output_weights_colname not in ms.colnames():
-                desc = ms.getcoldesc('WEIGHT_SPECTRUM')
-                desc['name'] = output_weights_colname
-                ms.addcols(desc)
-
-            ms.putcol(output_data_colname, all_data_list[msindex])
-            ms.putcol('FLAG', all_flags_list[msindex]) # this saves flags of nans, which is always good
-            ms.putcol(output_weights_colname, all_weights_list[msindex])
-            ms.close()
-        print("BLavg_multi: Finished one group of measurement sets.")
-
-
-def unwrap_fft(phase, iterations=3):
-    """
-    Unwrap phase using Fourier techniques.
-
-    For details, see:
-    Marvin A. Schofield & Yimei Zhu, Optics Letters, 28, 14 (2003)
-
-    Keyword arguments:
-    phase -- array of phase solutions
-    iterations -- number of iterations to perform
-    """
-    puRadius=lambda x : np.roll( np.roll(
-          np.add.outer( np.arange(-x.shape[0]/2+1,x.shape[0]/2+1)**2.0,
-                        np.arange(-x.shape[1]/2+1,x.shape[1]/2+1)**2.0 ),
-          x.shape[1]/2+1,axis=1), x.shape[0]/2+1,axis=0)+1e-9
-
-    idt,dt=np.fft.ifft2,np.fft.fft2
-    puOp=lambda x : idt( np.where(puRadius(x)==1e-9,1,puRadius(x)**-1.0)*dt(
-          np.cos(x)*idt(puRadius(x)*dt(np.sin(x)))
-         -np.sin(x)*idt(puRadius(x)*dt(np.cos(x))) ) )
-
-    def phaseUnwrapper(ip):
-       mirrored=np.zeros([x*2 for x in ip.shape])
-       mirrored[:ip.shape[0],:ip.shape[1]]=ip
-       mirrored[ip.shape[0]:,:ip.shape[1]]=ip[::-1,:]
-       mirrored[ip.shape[0]:,ip.shape[1]:]=ip[::-1,::-1]
-       mirrored[:ip.shape[0],ip.shape[1]:]=ip[:,::-1]
-
-       return (ip+2*np.pi*
-             np.round((puOp(mirrored).real[:ip.shape[0],:ip.shape[1]]-ip)
-             /2/np.pi))
-
-    i = 0
-    if iterations < 1:
-        interations = 1
-    while i < iterations:
-        i += 1
-        phase = phaseUnwrapper(phase)
-
-    return phase
-
-
-def input2bool(invar):
-    if isinstance(invar, bool):
-        return invar
-    elif isinstance(invar, str):
-        if invar.upper() == 'TRUE' or invar == '1':
-            return True
-        elif invar.upper() == 'FALSE' or invar == '0':
-            return False
-        else:
-            raise ValueError('input2bool: Cannot convert string "'+invar+'" to boolean!')
-    elif isinstance(invar, int) or isinstance(invar, float):
-        return bool(invar)
-    else:
-        raise TypeError('input2bool: Unsupported data type:'+str(type(invar)))
-
-
-def input2strlist(invar):
-    str_list = None
-    if type(invar) is str:
-        if invar.startswith('[') and invar.endswith(']'):
-            str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
-        else:
-            map_in = DataMap.load(invar)
-            map_in.iterator = DataMap.SkipIterator
-            str_list = []
-            for fname in map_in:
-                if fname.startswith('[') and fname.endswith(']'):
-                    for f in fname.strip('[]').split(','):
-                        str_list.append(f.strip(' \'\"'))
-                else:
-                    str_list.append(fname.strip(' \'\"'))
-    elif type(invar) is list:
-        str_list = [str(f).strip(' \'\"') for f in invar]
-    else:
-        raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
-    return str_list
-
-
-if __name__ == '__main__':
-    descriptiontext = "Pre-average data.\n"
-
-    parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
-    parser.add_argument('ms_file_pattern', help='Glob-able filename-pattern of input datasets')
-    parser.add_argument('parmdb_file_pattern', help='Glob-able filename-pattern of input direction-independent selfcal instrument parmdbs')
-    parser.add_argument('input_colname', help='Name of input column to pre-average')
-    parser.add_argument('output_data_colname', help='Name of output column')
-    parser.add_argument('output_weights_colname', help='Name of output column')
-    parser.add_argument('target_rms', help='Target rms in Jy/beam')
-    args = parser.parse_args()
-
-    ms_input = glob.glob(args.ms_file_pattern)
-    parmdb_input = glob.glob(args.parmdb_file_pattern)
-
-    main(ms_input, parmdb_input, args.input_colname, args.output_data_colname,
-        args.output_weights_colname, args.target_rms)
diff --git a/rapthor/scripts/reset_h5parm.py b/rapthor/scripts/reset_h5parm.py
deleted file mode 100755
index e97cd6f40797ff4c12792e862729cc15f08c28ef..0000000000000000000000000000000000000000
--- a/rapthor/scripts/reset_h5parm.py
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to reset solutions for given stations in an h5parm
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-from losoto.h5parm import h5parm
-from losoto.operations import reset, reweight
-import sys
-import os
-
-
-def main(h5in, baseline, solset='sol000', soltab='tec000', dataval=0.0):
-    """
-    Reset solutions
-
-    Parameters
-    ----------
-    h5in : str
-        Filename of h5parm
-    baseline : str
-        Baseline selection used in calibration. E.g.:
-            ""[CR]*&&;!RS106HBA;!RS205HBA;!RS208HBA"
-    solset : str, optional
-        Name of solset
-    soltab : str, optional
-        Name of soltab
-    dataval : float, optional
-        Value to set solutions to
-
-    """
-    h = h5parm(h5in, readonly=False)
-    s = h.getSolset(solset)
-    t = s.getSoltab(soltab)
-
-    # Make list of excluded stations and set selection to them
-    remotelist = [stat for stat in t.ant if '!{}'.format(stat) in baseline]
-    t.setSelection(ant=remotelist)
-
-    # Reset the values
-    reset.run(t, dataVal=float(dataval))
-
-    # Reset the weights
-    reweight.run(t, weightVal=1.0)
-
-    h.close()
-
-
-if __name__ == '__main__':
-    descriptiontext = "Reset solutions.\n"
-
-    parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
-    parser.add_argument('h5', help='name of input h5')
-    parser.add_argument('baseline', help='baseline selection')
-    args = parser.parse_args()
-
-    main(args.h5, args.baseline)
diff --git a/rapthor/scripts/reweight.py b/rapthor/scripts/reweight.py
deleted file mode 100755
index 3eea3fe37fd6471175b886620985ddce459944fc..0000000000000000000000000000000000000000
--- a/rapthor/scripts/reweight.py
+++ /dev/null
@@ -1,323 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to reweight uv data
-
-Based on https://github.com/ebonnassieux/Scripts/blob/master/QualityWeightsLOFAR.py
-"""
-from casacore.tables import table
-import numpy as np
-import sys
-import warnings
-from rapthor.lib import miscellaneous as misc
-
-
-class CovWeights:
-    def __init__(self, MSName, solint_sec, solint_hz, uvcut=[0, 2000], gainfile=None,
-                 phaseonly=False, dirname=None, quiet=True):
-        if MSName[-1] == "/":
-            self.MSName = MSName[0:-1]
-        else:
-            self.MSName = MSName
-        tab = table(self.MSName, ack=False)
-        self.timepersample = tab.getcell('EXPOSURE', 0)
-        self.ntSol = max(1, int(round(solint_sec / self.timepersample)))
-        tab.close()
-        sw = table(self.MSName+'::SPECTRAL_WINDOW', ack=False)
-        self.referencefreq = sw.col('REF_FREQUENCY')[0]
-        self.channelwidth = sw.col('CHAN_WIDTH')[0][0]
-        self.numchannels = sw.col('NUM_CHAN')[0]
-        sw.close()
-        self.nchanSol = max(1, self.get_nearest_frequstep(solint_hz / self.channelwidth))
-        self.uvcut = uvcut
-        self.gainfile = gainfile
-        self.phaseonly = phaseonly
-        self.dirname = dirname
-        self.quiet = quiet
-
-    def FindWeights(self, colname=""):
-        ms = table(self.MSName, ack=False)
-        ants = table(ms.getkeyword("ANTENNA"), ack=False)
-        antnames = ants.getcol("NAME")
-        ants.close()
-        nAnt = len(antnames)
-
-        u, v, _ = ms.getcol("UVW").T
-        A0 = ms.getcol("ANTENNA1")
-        A1 = ms.getcol("ANTENNA2")
-        tarray = ms.getcol("TIME")
-        nbl = np.where(tarray == tarray[0])[0].size
-        warnings.filterwarnings("ignore")
-        warnings.filterwarnings("default")
-        if "RESIDUAL_DATA" not in ms.colnames():
-            print('reweight: RESIDUAL_DATA not found. Exiting...')
-            sys.exit(1)
-        residualdata = ms.getcol("RESIDUAL_DATA")
-        flags = ms.getcol("FLAG")
-        ms.close()
-
-        # apply uvcut
-        c_m_s = 2.99792458e8
-        uvlen = np.sqrt(u**2 + v**2) / c_m_s * self.referencefreq
-        flags[uvlen > self.uvcut[1], :, :] = True
-        flags[uvlen < self.uvcut[0], :, :] = True
-        residualdata[flags] = np.nan
-        residualdata[residualdata == 0] = np.nan
-
-        # initialise
-        nChan = residualdata.shape[1]
-        nPola = residualdata.shape[2]
-        nt = residualdata.shape[0] / nbl
-        residualdata = residualdata.reshape((nt, nbl, nChan, nPola))
-        A0 = A0.reshape((nt, nbl))[0, :]
-        A1 = A1.reshape((nt, nbl))[0, :]
-
-        # make rms and residuals arrays
-        rmsarray = np.zeros((nt, nbl, nChan, 2), dtype=np.complex64)
-        residuals = np.zeros_like(rmsarray, dtype=np.complex64)
-        rmsarray[:, :, :, 0] = residualdata[:, :, :, 1]
-        rmsarray[:, :, :, 1] = residualdata[:, :, :, 2]
-        residuals[:, :, :, 0] = residualdata[:, :, :, 0]
-        residuals[:, :, :, 1] = residualdata[:, :, :, 3]
-
-        # start calculating the weights
-        CoeffArray = np.zeros((nt, nChan, nAnt))
-        ant1 = np.arange(nAnt)
-        num_sols_time = int(np.ceil(float(nt) / self.ntSol))
-        num_sols_freq = int(np.ceil(float(nChan) / self.nchanSol))
-        tcellsize = self.ntSol
-        for t_i in range(0, nt, self.ntSol):
-            if (t_i == nt - self.ntSol) and (nt % self.ntSol > 0):
-                tcellsize = nt % self.ntSol
-            t_e = t_i + tcellsize
-            fcellsize = self.nchanSol
-            for f_i in range(0, nChan, self.nchanSol):
-                if (f_i == nChan - self.nchanSol) and (nChan % self.nchanSol > 0):
-                    fcellsize = nChan % self.nchanSol
-                f_e = f_i + fcellsize
-
-                # build weights for each antenna in the current time-frequency block
-                for ant in ant1:
-                    # set of vis for baselines ant-ant_i
-                    set1 = np.where(A0 == ant)[0]
-                    # set of vis for baselines ant_i-ant
-                    set2 = np.where(A1 == ant)[0]
-                    CoeffArray[t_i:t_e, f_i:f_e, ant] = np.sqrt(
-                        np.nanmean(
-                            np.append(residuals[t_i:t_e, set1, f_i:f_e, :],
-                                      residuals[t_i:t_e, set2, f_i:f_e, :]) *
-                            np.append(residuals[t_i:t_e, set1, f_i:f_e, :],
-                                      residuals[t_i:t_e, set2, f_i:f_e, :]).conj()
-                            ) -
-                        np.nanstd(
-                            np.append(rmsarray[t_i:t_e, set1, f_i:f_e, :],
-                                      rmsarray[t_i:t_e, set2, f_i:f_e, :])
-                            )
-                        )
-            if not self.quiet:
-                PrintProgress(t_i, nt)
-
-        # plot
-#         Nr = 8
-#         Nc = 8
-#         xvals = range(nt)
-#         yvals = range(nChan)
-#         figSize = [10+3*Nc, 8+2*Nr]
-#         figgrid, axa = plt.subplots(Nr, Nc, sharex=True, sharey=True, figsize=figSize)
-#         for i in range(nAnt):
-#             ax = axa.flatten()[i]
-#             bbox = ax.get_window_extent().transformed(figgrid.dpi_scale_trans.inverted())
-#             aspect = ((xvals[-1]-xvals[0])*bbox.height)/((yvals[-1]-yvals[0])*bbox.width)
-#             im = ax.imshow(CoeffArray[:, :, i].transpose(1,0), origin='lower', interpolation="none", cmap=plt.cm.rainbow, norm=None,
-#                             extent=[xvals[0],xvals[-1],yvals[0],yvals[-1]], aspect=str(aspect))
-#         figgrid.colorbar(im, ax=axa.ravel().tolist(), use_gridspec=True, fraction=0.02, pad=0.005, aspect=35)
-#         figgrid.savefig(self.MSName+'_coef1.png', bbox_inches='tight')
-#         plt.close()
-
-        # get rid of NaNs and low values
-        CoeffArray[np.isnan(CoeffArray)] = np.inf
-        for i in range(nAnt):
-            tempars = CoeffArray[:, :, i]
-            thres = 0.25 * np.median(tempars[np.where(np.isfinite(tempars))])
-            CoeffArray[:, :, i][tempars < thres] = thres
-        return CoeffArray
-
-    def SaveWeights(self, CoeffArray, colname=None):
-        ms = table(self.MSName, readonly=False, ack=False)
-        ants = table(ms.getkeyword("ANTENNA"), ack=False)
-        antnames = ants.getcol("NAME")
-        nAnt = len(antnames)
-        tarray = ms.getcol("TIME")
-        darray = ms.getcol("DATA")
-        tvalues = np.array(sorted(list(set(tarray))))
-        nt = tvalues.shape[0]
-        nbl = tarray.shape[0]/nt
-        nchan = darray.shape[1]
-        npol = darray.shape[2]
-        A0 = np.array(ms.getcol("ANTENNA1").reshape((nt, nbl)))
-        A1 = np.array(ms.getcol("ANTENNA2").reshape((nt, nbl)))
-        warnings.filterwarnings("ignore")
-        warnings.filterwarnings("default")
-
-        # initialize weight array
-        w = np.zeros((nt, nbl, nchan, npol))
-        A0ind = A0[0, :]
-        A1ind = A1[0, :]
-
-        # do gains stuff
-        ant1gainarray, ant2gainarray = readGainFile(self.gainfile, ms, nt, nchan, nbl,
-                                                    tarray, nAnt, self.MSName, self.phaseonly,
-                                                    self.dirname)
-        ant1gainarray = ant1gainarray.reshape((nt, nbl, nchan))
-        ant2gainarray = ant2gainarray.reshape((nt, nbl, nchan))
-        for t in range(nt):
-            for i in range(nbl):
-                for j in range(nchan):
-                    w[t, i, j, :] = 1.0 / (CoeffArray[t, j, A0ind[i]] * ant1gainarray[t, i, j] +
-                                           CoeffArray[t, j, A1ind[i]] * ant2gainarray[t, i, j] +
-                                           CoeffArray[t, j, A0ind[i]] * CoeffArray[t, j, A1ind[i]] +
-                                           0.1)
-            if not self.quiet:
-                PrintProgress(t, nt)
-
-        # normalize
-        w = w.reshape(nt*nbl, nchan, npol)
-        w[np.isinf(w)] = np.nan
-        w = w / np.nanmean(w)
-        w[np.isnan(w)] = 0
-
-        # save in weights column
-        if colname is not None:
-            ms.putcol(colname, w)
-        ants.close()
-        ms.close()
-
-
-    def get_nearest_frequstep(self, freqstep):
-        """
-        Gets the nearest frequstep
-
-        Parameters
-        ----------
-        freqstep : int
-            Target frequency step
-
-        Returns
-        -------
-        optimum_step : int
-            Optimum frequency step nearest to target step
-        """
-        # Generate a list of possible values for freqstep
-        if not hasattr(self, 'freq_divisors'):
-            tmp_divisors = []
-            for step in range(self.numchannels, 0, -1):
-                if (self.numchannels % step) == 0:
-                    tmp_divisors.append(step)
-            self.freq_divisors = np.array(tmp_divisors)
-
-        # Find nearest
-        idx = np.argmin(np.abs(self.freq_divisors - freqstep))
-
-        return self.freq_divisors[idx]
-
-
-def readGainFile(gainfile, ms, nt, nchan, nbl, tarray, nAnt, msname, phaseonly, dirname):
-    if phaseonly:
-        ant1gainarray1 = np.ones((nt*nbl, nchan))
-        ant2gainarray1 = np.ones((nt*nbl, nchan))
-    else:
-        import losoto.h5parm
-        solsetName = "sol000"
-        soltabName = "screenamplitude000"
-        try:
-            gfile = losoto.h5parm.openSoltab(gainfile, solsetName=solsetName, soltabName=soltabName)
-        except:
-            print("Could not find amplitude gains in h5parm. Assuming gains of 1 everywhere.")
-            ant1gainarray1 = np.ones((nt*nbl, nchan))
-            ant2gainarray1 = np.ones((nt*nbl, nchan))
-            return ant1gainarray1, ant2gainarray1
-
-        freqs = table(msname+"/SPECTRAL_WINDOW").getcol("CHAN_FREQ")
-        gains = gfile.val  # axes: times, freqs, ants, dirs, pols
-        flagged = np.where(gains == 0.0)
-        gains[flagged] = np.nan
-        gfreqs = gfile.freq
-        times = gfile.time
-        dindx = gfile.dir.tolist().index(dirname)
-        ant1gainarray = np.zeros((nt*nbl, nchan))
-        ant2gainarray = np.zeros((nt*nbl, nchan))
-        A0arr = ms.getcol("ANTENNA1")
-        A1arr = ms.getcol("ANTENNA2")
-        deltime = (times[1] - times[0]) / 2.0
-        delfreq = (gfreqs[1] - gfreqs[0]) / 2.0
-        for i in range(len(times)):
-            timemask = (tarray >= times[i]-deltime) * (tarray < times[i]+deltime)
-            if np.all(~timemask):
-                continue
-            for j in range(nAnt):
-                mask1 = timemask * (A0arr == j)
-                mask2 = timemask * (A1arr == j)
-                for k in range(nchan):
-                    chan_freq = freqs[0, k]
-                    freqmask = np.logical_and(gfreqs >= chan_freq-delfreq,
-                                              gfreqs < chan_freq+delfreq)
-                    if chan_freq < gfreqs[0]:
-                        freqmask[0] = True
-                    if chan_freq > gfreqs[-1]:
-                        freqmask[-1] = True
-                    ant1gainarray[mask1, k] = np.nanmean(gains[i, freqmask, j, dindx, :], axis=(0, 1))
-                    ant2gainarray[mask2, k] = np.nanmean(gains[i, freqmask, j, dindx, :], axis=(0, 1))
-        ant1gainarray1 = ant1gainarray**2
-        ant2gainarray1 = ant2gainarray**2
-
-    return ant1gainarray1, ant2gainarray1
-
-
-def PrintProgress(currentIter, maxIter, msg=""):
-    sys.stdout.flush()
-    if msg == "":
-        msg = "Progress:"
-    sys.stdout.write("\r%s %5.1f %% " % (msg, 100*(currentIter+1.)/maxIter))
-    if currentIter == (maxIter-1):
-        sys.stdout.write("\n")
-
-
-def main(msname, solint_sec, solint_hz, colname="CAL_WEIGHT", gainfile="", uvcut_min=80.0,
-         uvcut_max=1e6, phaseonly=True, dirname=None, quiet=True):
-    """
-    Reweight visibilities
-
-    Parameters
-    ----------
-    filename : str
-        Name of the input measurement set
-    solint_sec : float
-        Solution interval in seconds of calibration
-    solint_hz : float
-        Solution interval in Hz of calibration
-    colname : str, optional
-        Name of the weights column name you want to save the weights to
-    gainfile : str, optional
-        Name of the gain file you want to read to rebuild the calibration quality weights.
-        If no file is given, equivalent to rebuilding weights for phase-only calibration
-    uvcut_min : float, optional
-        Min uvcut in lambda used during calibration
-    uvcut_max : float, optional
-        Max uvcut in lambda used during calibration
-    phaseonly : bool, optional
-        Use if calibration was phase-only; this means that gain information doesn't need
-        to be read
-    dirname : str, optional
-        Name of calibration patch
-    """
-    solint_sec = float(solint_sec)
-    solint_hz = float(solint_hz)
-    uvcut_min = float(uvcut_min)
-    uvcut_max = float(uvcut_max)
-    uvcut = [uvcut_min, uvcut_max]
-    phaseonly = misc.string2bool(phaseonly)
-
-    covweights = CovWeights(MSName=msname, solint_sec=solint_sec, solint_hz=solint_hz,
-                            gainfile=gainfile, uvcut=uvcut, phaseonly=phaseonly,
-                            dirname=dirname, quiet=quiet)
-    coefficients = covweights.FindWeights(colname=colname)
-    covweights.SaveWeights(coefficients, colname=colname)