diff --git a/rfistrategies/lofar-default.lua b/rfistrategies/lofar-default.lua
new file mode 100644
index 0000000000000000000000000000000000000000..810f306b075447274ded8c663c5dc6a9adff09c6
--- /dev/null
+++ b/rfistrategies/lofar-default.lua
@@ -0,0 +1,79 @@
+--[[
+ This is the LOFAR strategy. It is almost equal to the
+ generic "minimal" AOFlagger strategy, version 2020-06-14.
+ Author: André Offringa
+]]--
+
+function execute(input)
+
+  --
+  -- Generic settings
+  --
+
+  local base_threshold = 1.0  -- lower means more sensitive detection
+  -- How to flag complex values, options are: phase, amplitude, real, imaginary, complex
+  local representation = "amplitude"
+  local iteration_count = 3  -- how many iterations to perform?
+  local threshold_factor_step = 2.0 -- How much to increase the sensitivity each iteration?
+  local frequency_resize_factor = 3.0 -- Amount of "extra" smoothing in frequency direction
+  local transient_threshold_factor = 1.0 -- decreasing this value makes detection of transient RFI more aggressive
+
+  --
+  -- End of generic settings
+  --
+
+  local inpPolarizations = input:get_polarizations()
+
+  input:clear_mask()
+
+  -- For collecting statistics. Note that this is done after clear_mask(),
+  -- so that the statistics ignore any flags in the input data. 
+  local copy_of_input = input:copy()
+  
+  for ipol,polarization in ipairs(inpPolarizations) do
+ 
+    local data = input:convert_to_polarization(polarization)
+
+    data = data:convert_to_complex(representation)
+    local original_data = data:copy()
+
+    for i=1,iteration_count-1 do
+      local threshold_factor = threshold_factor_step^(iteration_count-i)
+
+      local sumthr_level = threshold_factor * base_threshold
+      aoflagger.sumthreshold(data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+ 
+      -- Do timestep & channel flagging
+      local chdata = data:copy()
+      aoflagger.threshold_timestep_rms(data, 3.5)
+      aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true)
+      data:join_mask(chdata)
+
+      -- High pass filtering steps
+      data:set_visibilities(original_data)
+      local resized_data = aoflagger.downsample(data, 3, frequency_resize_factor, true)
+      aoflagger.low_pass_filter(resized_data, 21, 31, 2.5, 5.0)
+      aoflagger.upsample(resized_data, data, 3, frequency_resize_factor)
+
+      local tmp = original_data - data
+      tmp:set_mask(data)
+      data = tmp
+
+      aoflagger.set_progress((ipol-1)*iteration_count+i, #inpPolarizations*iteration_count )
+    end -- end of iterations
+
+    aoflagger.sumthreshold(data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+
+    if input:is_complex() then
+      data = data:convert_to_complex("complex")
+    end
+    input:set_polarization_data(polarization, data)
+
+    aoflagger.set_progress(ipol, #inpPolarizations )
+  end -- end of polarization iterations
+
+  aoflagger.scale_invariant_rank_operator(input, 0.2, 0.2)
+  aoflagger.threshold_timestep_rms(input, 4.0)
+  input:flag_nans()
+end
+
diff --git a/rfistrategies/lofar-lba-wideband.lua b/rfistrategies/lofar-lba-wideband.lua
new file mode 100644
index 0000000000000000000000000000000000000000..d514af8760a6c956c12acd825d769c4db239525a
--- /dev/null
+++ b/rfistrategies/lofar-lba-wideband.lua
@@ -0,0 +1,143 @@
+--[[
+ This is a LBA AOFlagger strategy for combined subbands, version 2021-03-30
+ Author: André Offringa
+]]--
+
+aoflagger.require_min_version("3.1")
+
+function execute(input)
+
+  --
+  -- Generic settings
+  --
+
+  -- What polarizations to flag? Default: input:get_polarizations() (=all that are in the input data)
+  -- Other options are e.g.:
+  -- { 'XY', 'YX' } to flag only XY and YX, or
+  -- { 'I', 'Q' } to flag only on Stokes I and Q
+  local flag_polarizations = input:get_polarizations()
+
+  local base_threshold = 1.2  -- lower means more sensitive detection
+  -- How to flag complex values, options are: phase, amplitude, real, imaginary, complex
+  -- May have multiple values to perform detection multiple times
+  local flag_representations = { "amplitude" }
+  local iteration_count = 5  -- how many iterations to perform?
+  local threshold_factor_step = 2.0 -- How much to increase the sensitivity each iteration?
+  -- If the following variable is true, the strategy will consider existing flags
+  -- as bad data. It will exclude flagged data from detection, and make sure that any existing
+  -- flags on input will be flagged on output. If set to false, existing flags are ignored.
+  local exclude_original_flags = true
+  local transient_threshold_factor = 1.0 -- decreasing this value makes detection of transient RFI more aggressive
+ 
+  --
+  -- End of generic settings
+  --
+
+  local inpPolarizations = input:get_polarizations()
+
+  if(not exclude_original_flags) then
+    input:clear_mask()
+  end
+  -- For collecting statistics. Note that this is done after clear_mask(),
+  -- so that the statistics ignore any flags in the input data. 
+  local copy_of_input = input:copy()
+
+  aoflagger.normalize_bandpass(input)
+  
+  for ipol,polarization in ipairs(flag_polarizations) do
+ 
+    local pol_data = input:convert_to_polarization(polarization)
+    local original_data
+
+    for _,representation in ipairs(flag_representations) do
+
+      data = pol_data:convert_to_complex(representation)
+      original_data = data:copy()
+
+      for i=1,iteration_count-1 do
+        local threshold_factor = math.pow(threshold_factor_step, iteration_count-i)
+
+        local sumthr_level = threshold_factor * base_threshold
+        if(exclude_original_flags) then
+          aoflagger.sumthreshold_masked(data, original_data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+        else
+          aoflagger.sumthreshold(data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+        end
+
+        -- Do timestep & channel flagging
+        local chdata = data:copy()
+        aoflagger.threshold_timestep_rms(data, 3.5)
+        aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true)
+        data:join_mask(chdata)
+
+        -- High pass filtering steps
+        data:set_visibilities(original_data)
+        if(exclude_original_flags) then
+          data:join_mask(original_data)
+        end
+
+        aoflagger.low_pass_filter(data, 21, 31, 1.0, 1.0)
+
+        -- In case this script is run from inside rfigui, calling
+        -- the following visualize function will add the current result
+        -- to the list of displayable visualizations.
+        -- If the script is not running inside rfigui, the call is ignored.
+        aoflagger.visualize(data, "Fit #"..i, i-1)
+
+        local tmp = original_data - data
+        tmp:set_mask(data)
+        data = tmp
+
+        aoflagger.visualize(data, "Residual #"..i, i+iteration_count)
+        aoflagger.set_progress((ipol-1)*iteration_count+i, #flag_polarizations*iteration_count )
+      end -- end of iterations
+
+      if(exclude_original_flags) then
+        aoflagger.sumthreshold_masked(data, original_data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+      else
+        aoflagger.sumthreshold(data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+      end
+    end -- end of complex representation iteration
+
+    if(exclude_original_flags) then
+      data:join_mask(original_data)
+    end
+
+    -- Helper function used below
+    function contains(arr, val)
+      for _,v in ipairs(arr) do
+        if v == val then return true end
+      end
+      return false
+    end
+
+    if contains(inpPolarizations, polarization) then
+      if input:is_complex() then
+        data = data:convert_to_complex("complex")
+      end
+      input:set_polarization_data(polarization, data)
+    else
+      input:join_mask(data)
+    end
+
+    aoflagger.visualize(data, "Residual #"..iteration_count, 2*iteration_count)
+    aoflagger.set_progress(ipol, #flag_polarizations )
+  end -- end of polarization iterations
+
+  if(exclude_original_flags) then
+    aoflagger.scale_invariant_rank_operator_masked(input, copy_of_input, 0.2, 0.2)
+  else
+    aoflagger.scale_invariant_rank_operator(input, 0.2, 0.2)
+  end
+
+  aoflagger.threshold_timestep_rms(input, 4.0)
+
+  if input:is_complex() and input:has_metadata() then
+    -- This command will calculate a few statistics like flag% and stddev over
+    -- time, frequency and baseline and write those to the MS. These can be
+    -- visualized with aoqplot.
+    aoflagger.collect_statistics(input, copy_of_input)
+  end
+  input:flag_nans()
+end
+
diff --git a/scripts/TargetListToCoords.py b/scripts/TargetListToCoords.py
index 5a712f5b5f2220da0af1ec5794b2823215287e09..46654e2b54c0d781a2135cf6831a3f5ba9111b3b 100644
--- a/scripts/TargetListToCoords.py
+++ b/scripts/TargetListToCoords.py
@@ -27,8 +27,8 @@ def plugin_main(**kwargs):
     
     # read in the catalogue to get source_id, RA, and DEC
     t = Table.read(target_file, format='csv')
-    RA_val = t['RA_LOTSS'].data[0]
-    DEC_val = t['DEC_LOTSS'].data[0]
+    RA_val = t['RA'].data[0]
+    DEC_val = t['DEC'].data[0]
     Source_id = t['Source_id'].data[0]
     if str(Source_id)[0:1] == 'I':
         pass
diff --git a/scripts/check_Ateam_separation.py b/scripts/check_Ateam_separation.py
index d4ae692f5bacb7c7c958498845dbe5ceae219cbc..d1549a8cd44bfbf8a38d007a1bb2a6eecadbc8e2 100755
--- a/scripts/check_Ateam_separation.py
+++ b/scripts/check_Ateam_separation.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
 #
diff --git a/scripts/downloadCats.py b/scripts/downloadCats.py
index 3a69db79978f2dcb64a6d5a6569a4cecc25ee2c1..1f7121a185a405d9f904937d0d1144aba2ef53b2 100755
--- a/scripts/downloadCats.py
+++ b/scripts/downloadCats.py
@@ -1,21 +1,13 @@
 #!/usr/bin/env python3
-#from scipy import stats  
 import os, sys, logging, io
 import numpy as np
-#import csv
-#import matplotlib
-#matplotlib.use('Agg')
-#import matplotlib.pyplot as plt
-#from lmfit import SkewedGaussianModel
-
 import pyvo as vo
 import pyrap.tables as pt
 from astropy.table import Table, Column, vstack, unique, hstack
 import argparse
-#from lofarpipe.support.data_map import DataMap
-#from lofarpipe.support.data_map import DataProduct
-import requests
 from time import sleep
+
+import requests
 from astropy.coordinates import SkyCoord
 from astropy import units as u
 
@@ -27,8 +19,8 @@ def sum_digits(n):
         if ii != '-':
             s += int(ii)
             c += 1
-    norm = float(s)/float(c)
-    return( norm )
+    norm = float(s) / float(c)
+    return norm
 
 def count_p(n):
     s = 0
@@ -38,8 +30,8 @@ def count_p(n):
             c += 1
             if ii == 'P':
                 s += 1
-    norm = float(s)/float(c)
-    return(norm)
+    norm = float(s) / float(c)
+    return norm
 
 def count_s(n):
     s = 0
@@ -49,8 +41,8 @@ def count_s(n):
             c += 1
             if ii == 'S':
                 s += 1
-    norm = float(s)/float(c)
-    return(norm)
+    norm = float(s) / float(c)
+    return norm
 
 def count_x(n):
     s = 0
@@ -60,8 +52,8 @@ def count_x(n):
             c += 1
             if ii == 'X':
                 s += 1
-    norm = float(s)/float(c)
-    return(norm)
+    norm = float(s) / float(c)
+    return norm
 
 def grab_coo_MS(MS):
     """
@@ -77,7 +69,7 @@ def grab_coo_MS(MS):
     """
 
     # reading the coordinates ("position") from the MS
-    # NB: they are given in rad,rad (J2000) 
+    # NB: they are given in rad,rad (J2000)
     [[[ra,dec]]] = pt.table(MS+'/FIELD', readonly=True, ack=False).getcol('PHASE_DIR')
     # [[[ra,dec]]] = pt.table(MS, readonly=True, ack=False).getcol('PHASE_DIR')
 
@@ -120,7 +112,7 @@ def my_lotss_catalogue( ms_input, Radius=1.5, bright_limit_Jy=5., outfile='' ):
         String from the list (map) of the target MSs
     Radius : float (default = 1.5)
         Radius for the LOTSS cone search in degrees
-    
+
     """
     ## first check if the file already exists
     if os.path.isfile( outfile ):
@@ -138,7 +130,8 @@ def my_lotss_catalogue( ms_input, Radius=1.5, bright_limit_Jy=5., outfile='' ):
         ## this is the tier 1 database to query
         #url = 'http://vo.astron.nl/lofartier1/q/cone/scs.xml'
         # HETDEX database.
-        url = 'https://vo.astron.nl/hetdex/lotss-dr1/cone/scs.xml'
+        #url = 'https://vo.astron.nl/hetdex/lotss-dr1/cone/scs.xml'
+        url = 'https://vo.astron.nl/lotss_dr2/q/src_cone/scs.xml'
 
         ## query the database
         query = vo.dal.scs.SCSQuery( url )
@@ -161,10 +154,10 @@ def my_lotss_catalogue( ms_input, Radius=1.5, bright_limit_Jy=5., outfile='' ):
         if 'Resolved' not in colnames:
             resolved = np.where(is_resolved(tb_sorted['Total_flux'], tb_sorted['Peak_flux'], tb_sorted['Isl_rms']), 'R', 'U')
             tb_sorted['Resolved'] = resolved
-        ## rename source_id column if necessary   
+        ## rename source_id column if necessary
         if 'Source_Name' in tb_sorted.colnames:
             tb_sorted.rename_column('Source_Name', 'Source_id')
-        keep_cols = ['Source_id', 'RA', 'DEC','Total_flux','Peak_flux', 'Major', 'Minor', 'PA', 'DC_Maj', 'DC_Min', 'DC_PA', 'Isl_rms', 'Resolved']
+        keep_cols = ['Source_id', 'RA', 'DEC','Total_flux','Peak_flux', 'Majax', 'Minax', 'PA', 'DC_Maj', 'DC_Min', 'DC_PA', 'Isl_rms', 'Resolved']
         if 'LGZ_Size' in colnames:
             keep_cols = keep_cols + ['LGZ_Size', 'LGZ_Width', 'LGZ_PA']
 
@@ -184,7 +177,7 @@ def my_lbcs_catalogue( ms_input, Radius=1.5, outfile='' ):
         String from the list (map) of the target MSs
     Radius : float (default = 1.5)
         Radius for the LOTSS cone search in degrees
-    
+
     """
     ## first check if the file already exists
     print( outfile )
@@ -205,7 +198,7 @@ def my_lbcs_catalogue( ms_input, Radius=1.5, outfile='' ):
                 response = requests.get(url, stream=True,verify=True,timeout=60)
                 if response.status_code!=200:
                     print(response.headers)
-                    raise RuntimeError('Code was %i' % response.status_code)
+                    raise RuntimeError(f'Code was {response.status_code}')
             except requests.exceptions.ConnectionError:
                 print('Connection error! sleeping 30 seconds before retry...')
                 sleep(30)
@@ -227,7 +220,7 @@ def my_lbcs_catalogue( ms_input, Radius=1.5, outfile='' ):
         print( len( tb ) )
 
         #### Workaround to sort and pick good calibrator info from tb array ###########
-        
+
         if not len(tb) > 1:
             logging.critical('There are no LBCS sources within the given radius. Check your source is within the LBCS footprint and increase the search radius. Exiting...')
             return
@@ -235,7 +228,7 @@ def my_lbcs_catalogue( ms_input, Radius=1.5, outfile='' ):
             ## calculate the total FT goodness
             ft_total = []
             for xx in range(len(tb)):
-                    ft_total.append( sum_digits( tb[xx]['FT_Goodness'] ) )
+                ft_total.append( sum_digits( tb[xx]['FT_Goodness'] ) )
             ft_col = Column( ft_total, name='FT_total' )
             tb.add_column( ft_col )
 
@@ -262,12 +255,12 @@ def find_close_objs(lo, lbcs, tolerance=5.):
     lotss_coords = SkyCoord( lo['RA'], lo['DEC'], frame='icrs', unit='deg' )
     lbcs_coords = SkyCoord( lbcs['RA'], lbcs['DEC'], frame='icrs', unit='deg' )
 
-    ## search radius 
+    ## search radius
     search_rad = 5. / 60. / 60. * u.deg
 
     ## loop through the lbcs coordinates -- this will be much faster than looping through lotss
     lotss_idx = []
-    lbcs_idx = []    
+    lbcs_idx = []
     for xx in range(len(lbcs)):
         seps = lbcs_coords[xx].separation(lotss_coords)
         match_idx = np.where( seps < search_rad )[0]
@@ -288,7 +281,7 @@ def find_close_objs(lo, lbcs, tolerance=5.):
                 if not isinstance(m_idx,int):
                     m_idx = m_idx[0]
                 lbcs_idx.append(xx)
-                lotss_idx.append(m_idx) 
+                lotss_idx.append(m_idx)
     lbcs_matches = lbcs[lbcs_idx]
     lotss_matches = lo[lotss_idx]
 
@@ -302,7 +295,7 @@ def find_close_objs(lo, lbcs, tolerance=5.):
         for src_id in src_ids:
             idx = np.where( combined['Source_id'] == src_id )[0]
             if len(idx) > 1:
-                ## multiple matches found.  Count P's first and then break ties with Goodness_FT 
+                ## multiple matches found.  Count P's first and then break ties with Goodness_FT
                 num_P = []
                 total_ft = []
                 for yy in range( len( idx ) ):
@@ -317,7 +310,7 @@ def find_close_objs(lo, lbcs, tolerance=5.):
                     if len( best_idx ) == 1:
                         good_idx.append(idx[best_idx][0])  ## idx[best_idx][0] is a number
                     if len( best_idx ) > 1:
-                        currentmax = 0.0 
+                        currentmax = 0.0
                         for i in range(0,len(best_idx)):
                             if total_ft[best_idx[i]] > currentmax:
                                 currentmax = total_ft[best_idx[i]]
@@ -335,36 +328,84 @@ def find_close_objs(lo, lbcs, tolerance=5.):
     ## rename RA columns
     result.rename_column('RA_1','RA_LBCS')
     result.rename_column('DEC_1','DEC_LBCS')
-    result.rename_column('RA_2','RA_LOTSS')
-    result.rename_column('DEC_2','DEC_LOTSS')
+    result.rename_column('RA_2','RA')
+    result.rename_column('DEC_2','DEC')
 
     return result
 
+def is_resolved(Sint, Speak, rms):
+    """ Determines if a source is resolved or unresolved.
+    The calculation is presented in Shimwell et al. 2018 of the LOFAR DR1 paper splash.
+
+    Args:
+        Sint (float or ndarray): integrated flux density.
+        Speak (float or ndarray): peak flux density.
+        rms (float or ndarray): local rms around the source.
+    Returns:
+        resolved (bool or ndarray): True if the source is resolved, False if not.
+    """
+    resolved = ((Sint / Speak) > 1.25 + 3.1 * (Speak / rms) ** (-0.53))
+    return resolved
+
+def remove_multiples_position( mycat, racol='RA', decol='DEC' ):
+    radecstrings = []
+    for i in np.arange(0,len(mycat)):
+        radecstrings.append(str(mycat[i][racol]) + str(mycat[i][decol]) )
+    radecstrings = np.asarray(radecstrings)
+    if len( np.unique( radecstrings ) ) != len( mycat ):
+        radecs = np.unique( radecstrings )
+        good_idx = []
+        for radec in radecs:
+            idx = np.where( radecstrings == radec )[0]
+            if len(idx) > 1:
+                ## multiple matches found.  Count P's first and then break ties with Goodness_FT
+                num_P = []
+                total_ft = []
+                for yy in range( len( idx ) ):
+                    tmp = mycat[idx[yy]]['Goodness']
+                    num_P.append( count_p( tmp ) )
+                    tmp = mycat[idx[yy]]['FT_Goodness']
+                    total_ft.append( sum_digits( tmp ) )
+                ## check that the total_ft values are non-zero before looking for a best cal
+                if np.max( total_ft ) > 0:
+                    ## pick the one with the highest number of P's -- if tie, use total_ft
+                    best_idx = np.where( num_P == np.max( num_P ) )[0]  ## this is an array
+                    if len( best_idx ) == 1:
+                        good_idx.append(idx[best_idx][0])  ## idx[best_idx][0] is a number
+                    if len( best_idx ) > 1:
+                        currentmax = 0.0
+                        for i in range(0,len(best_idx)):
+                            if total_ft[best_idx[i]] > currentmax:
+                                currentmax = total_ft[best_idx[i]]
+                                ft_idx = i
+                        good_idx.append( idx[best_idx[ft_idx]] )
+                else:
+                    print( 'Duplicate sources have total_ft = 0, removing from results.' )
+            else:
+                good_idx.append(idx[0])
+    else:
+        print( 'All LBCS sources are unique' )
+
+    return( mycat[good_idx] )
+
 def plugin_main( args, **kwargs ):
 
     MSname = args[0]
-    #mapfile_in = kwargs['mapfile_in']
     lotss_radius = kwargs['lotss_radius']
     lbcs_radius  = kwargs['lbcs_radius']
+    im_radius = float(kwargs['im_radius'])
+    image_limit_Jy = float(kwargs['image_limit_Jy'])
     bright_limit_Jy = float(kwargs['bright_limit_Jy'])
     lotss_result_file = kwargs['lotss_result_file']
     lotss_catalogue = kwargs['lotss_catalogue']
     lbcs_catalogue = kwargs['lbcs_catalogue']
     delay_cals_file = kwargs['delay_cals_file']
-    subtract_file = kwargs['subtract_file']
     match_tolerance = float(kwargs['match_tolerance'])
-    subtract_limit = float(kwargs['subtract_limit'])
-    image_limit_Jy = float(kwargs['image_limit_Jy'])
     fail_lotss_ok = kwargs['continue_no_lotss'].lower().capitalize()
 
-    #mslist = DataMap.load(mapfile_in)
-    #MSname = mslist[0].file
-    # For testing
-    #MSname = kwargs['MSname']
- 
     ## first check for a valid delay_calibrator file
     if os.path.isfile(delay_cals_file):
-        print( 'Delay calibrators file {:s} exists! returning.'.format(delay_cals_file) )
+        print(f'Delay calibrators file {delay_cals_file} exists! returning.')
         return
 
     ## look for or download LBCS
@@ -381,13 +422,13 @@ def plugin_main( args, **kwargs ):
         return
     if len(lotss_catalogue) == 0 and not fail_lotss_ok:
         logging.error('LoTSS coverage does not exist, and contine_without_lotss is set to False.')
-        return 
+        return
 
     ## if the LoTSS catalogue is empty, write out the delay cals only and stop
     if len(lotss_catalogue) == 0:
-        print('Target field not in LoTSS coverage yet! Only writing {:s} and best_{:s} based on LBCS'.format(delay_cals_file,delay_cals_file))
-        lbcs_catalogue.write(delay_cals_file, format='csv')
-        ## pick the best calibrator based on LBCS information only
+        print('Target field not in LoTSS coverage yet! Only writing {:s} based on LBCS'.format(delay_cals_file))
+
+        ## Add the radius from phase centre to the catalogue
         RATar, DECTar = grab_coo_MS(input2strlist_nomapfile(MSname)[0])
         ptg_coords = SkyCoord( RATar, DECTar, frame='icrs', unit='deg' )
 
@@ -396,30 +437,26 @@ def plugin_main( args, **kwargs ):
         seps = Column( separations.deg, name='Radius' )
         lbcs_catalogue.add_column( seps )
 
-        ## highest FT_total; if tie use closest to phase centre
-        best_idx = np.where( lbcs_catalogue['FT_total'] == np.max( lbcs_catalogue['FT_total'] ) )[0]
-        if len( best_idx ) > 1:
-            tmp = seps[best_idx]
-            new_best_idx = np.where( tmp == np.min(tmp) )[0]
-            best_idx = best_idx[new_best_idx]
-
-        ## rename some columns 
-        lbcs_catalogue.rename_column('RA','RA_LOTSS')
-        lbcs_catalogue.rename_column('DEC','DEC_LOTSS')
+        ## rename the source_id column
         lbcs_catalogue.rename_column('Observation','Source_id')
+
         ## add in some dummy data
-        Total_flux = Column( np.ones(len(lbcs_catalogue)), name='Total_flux' )
+        Total_flux = Column( np.ones(len(lbcs_catalogue))*1e3, name='Total_flux', unit='mJy' )
         lbcs_catalogue.add_column( Total_flux )
-        LGZ_Size = Column( np.ones( len(lbcs_catalogue) )*20., name='LGZ_Size' ) ## set to a default of 20 arcsec
+        LGZ_Size = Column( np.ones( len(lbcs_catalogue) )*20., name='LGZ_Size', unit='arcsec' ) ## set to a default of 20 arcsec
         lbcs_catalogue.add_column( LGZ_Size )
 
-        best_result = lbcs_catalogue[best_idx]
-        best_file = delay_cals_file.replace('delay_','best_delay_')
-        best_result.write( best_file, format='csv' )
-        print( 'Writing best delay calibrator information to file {:s}'.format(best_file) )
+        ## remove duplicate sources if necessary
+        lbcs_catalogue = remove_multiples_position( lbcs_catalogue )
+
+        ## order based on radius from the phase centre
+        lbcs_catalogue.sort('Radius')
+
+        ## write the catalogue
+        lbcs_catalogue.write(delay_cals_file, format='csv')
         return
 
-    ## else continue 
+    ## else continue
     result = find_close_objs( lotss_catalogue, lbcs_catalogue, tolerance=match_tolerance )
 
     ## check if there are any matches
@@ -427,127 +464,62 @@ def plugin_main( args, **kwargs ):
         logging.error('LoTSS and LBCS coverage exists, but no matches found. This indicates something went wrong, please check your catalogues.')
         return
     else:
-        # find the best delay calibrator
+        # add radius to the catalogue
         RATar, DECTar = grab_coo_MS(input2strlist_nomapfile(MSname)[0])
         ptg_coords = SkyCoord( RATar, DECTar, frame='icrs', unit='deg' )
 
-        src_coords = SkyCoord( result['RA_LBCS'], result['DEC_LBCS'], frame='icrs', unit='deg' )
+        src_coords = SkyCoord( result['RA'], result['DEC'], frame='icrs', unit='deg' )
         separations = src_coords.separation(ptg_coords )
         seps = Column( separations.deg, name='Radius' )
         result.add_column( seps )
 
-        radsq = np.array( result['Radius'] )**2.
-        fluxjy = np.array( result['Total_flux'] )*1e-3
-        ftsq = np.array( result['FT_total'] )**2.
-
-        mystat = [ radsq[xx]/fluxjy[xx]/ftsq[xx] if ftsq[xx] > 0 else 100 for xx in range(len(radsq)) ]
-        mycol = Column( mystat, name='Select_stat' )
-        result.add_column( mycol )
-
-        ## pick the best one
-        best_idx = np.where( mystat == np.min( mystat ) )[0]
-        tmp = result['Source_id'][best_idx].data
-        tmp2 = result['Select_stat'][best_idx].data
-        print( 'Best delay calibrator candidate is {:s} with a stastic of {:f}'.format(str(tmp[0]),tmp2[0]) )
+        ## order by radius from the phase centre
+        result.sort( 'Radius' )
 
         ## Write catalogues
         ## 1 - delay calibrators -- from lbcs_catalogue
         result.write( delay_cals_file, format='csv' )
-        print('Writing delay calibrator candidate file {:s}'.format(delay_cals_file))
-
-        ## best delay calibrator
-        best_result = result[best_idx]
-        best_file = delay_cals_file.replace('delay_','best_delay_')
-        best_result.write( best_file, format='csv' )
-        print( 'Writing best delay calibrator information to file {:s}'.format(best_file) )
-
-        ## convert Jy to milliJy
-        subtract_index = np.where( result['Total_flux'] > subtract_limit*1e3 )[0]
-        subtract_cals = result[subtract_index]
-
-        ## also include bright sources from the lotss catalogue
-        ## convert Jy to milliJy
-        bright_index = np.where( lotss_catalogue['Total_flux'] >= bright_limit_Jy*1e3 )[0]
-        lotss_bright = lotss_catalogue[bright_index]
-        ## lotss catalogue has units, redefine a table that doesn't
-
-        subtract_sources = vstack( [subtract_cals, lotss_bright], join_type='outer' )
-
-        unique_srcs = np.unique( subtract_sources['Source_id'] )
-        if len( unique_srcs ) != len( subtract_sources ):
-            ## remove duplicates, keep LBCS information
-            ## this is untested and will probably fail ...
-            good_idx = []
-            for src_id in unique_srcs:
-                idx = np.where( subtract_sources['Source_id'] == src_id )[0]
-                if len( idx ) > 1:
-                    tmp = subtract_sources[idx]
-                    lbcs_idx = np.where( tmp['RA_LBCS'] > 0 )[0]
-                    good_idx.append( idx[lbcs_idx][0] )
-                else:
-                    good_idx.append( idx[0] )
-            subtract_sources = subtract_sources[good_idx]
+        print(f'Writing delay calibrator candidate file {delay_cals_file}')
 
-        subtract_sources.write( subtract_file, format='csv' )
-
-        ## sources to image -- first remove things that are already in the delay_cals_file and subtract_file
-        good_index = [ x for x, src_id in enumerate( lotss_catalogue['Source_id'] ) if src_id not in result['Source_id'] and src_id not in subtract_sources['Source_id'] ]
+        ## sources to image -- first remove things that are already in the delay_cals_file
+        good_index = [ x for x, src_id in enumerate( lotss_catalogue['Source_id'] ) if src_id not in result['Source_id'] ]
 
         tmp_cat = lotss_catalogue[good_index]
 
         ## make a flux cut
         image_index = np.where( tmp_cat['Peak_flux'] >= image_limit_Jy*1e3 )[0]
-        sources_to_image = tmp_cat[image_index]
+        flux_cut_sources = tmp_cat[image_index]
 
-        ## find unresolved
-        nsrcs = float( len( sources_to_image ) )
-        print("There are "+str(nsrcs)+" sources above "+str(image_limit_Jy)+" mJy.")
-        try:
-            unresolved_index = np.where( sources_to_image['Resolved'] == 'U' )[0]
-        except:
-            unresolved_index = np.where( is_resolved(sources_to_image['Total_flux'],  sources_to_image['Peak_flux'], sources_to_image['Isl_rms'] ) )[0]
-        if nsrcs==0:
-            print("Did not find any unresolved objects.")
-        else:
-            perc_unres = len( unresolved_index ) / nsrcs * 100.
-            print('Percentage of sources which are unresolved: '+str( perc_unres ))
+        ## make a radius cut
+        src_coords = SkyCoord( flux_cut_sources['RA'], flux_cut_sources['DEC'], frame='icrs', unit='deg' )
+        separations = src_coords.separation( ptg_coords )
+        seps = Column( separations.deg, name='Radius' )
+        flux_cut_sources.add_column( seps )
+        good_idx = np.where( flux_cut_sources['Radius'] <= im_radius )[0]
+        sources_to_image = flux_cut_sources[good_idx]
 
+        nsrcs = float( len( sources_to_image ) )
+        print(f"There are {nsrcs} sources above {image_limit_Jy} mJy within {im_radius} degrees of the phase centre.")
         sources_to_image.write( lotss_result_file, format='csv' )
 
     return
 
-def is_resolved(Sint, Speak, rms):
-    """ Determines if a source is resolved or unresolved.
-    The calculation is presented in Shimwell et al. 2018 of the LOFAR DR1 paper splash.
-    
-    Args:
-        Sint (float or ndarray): integrated flux density.
-        Speak (float or ndarray): peak flux density.
-        rms (float or ndarray): local rms around the source.
-    Returns:
-        resolved (bool or ndarray): True if the source is resolved, False if not.
-    """
-    resolved = ((Sint / Speak) > 1.25 + 3.1 * (Speak / rms) ** (-0.53))
-    return resolved
-
 if __name__ == "__main__":
 
     parser = argparse.ArgumentParser()
-    parser.add_argument( '--lotss_radius', dest='lotss_radius', type=float, help='Radius to search LoTSS', default=5. )
-    parser.add_argument( '--lbcs_radius', dest='lbcs_radius', type=float, help='Radius to search LBCS', default=5. )
+    parser.add_argument( '--lotss_radius', dest='lotss_radius', type=float, help='Radius to search LoTSS', default=1.5 )
+    parser.add_argument( '--lbcs_radius', dest='lbcs_radius', type=float, help='Radius to search LBCS', default=1.5 )
+    parser.add_argument( '--im_radius', dest='im_radius', type=float, help='Radius in which to image', default=1.24 )
     parser.add_argument( '--lotss_catalogue', dest='lotss_catalogue', type=str, help='input file for LoTSS catalogue [will be downloaded if does not exist]', default='lotss_catalogue.csv' )
     parser.add_argument( '--lbcs_catalogue', dest='lbcs_catalogue', type=str, help='input file for LBCS catalogue [will be downloaded if does not exist]', default='lbcs_catalogue.csv' )
     parser.add_argument( '--lotss_result_file', dest='lotss_result_file', type=str, help='output file of sources to image', default='image_catalogue.csv' )
     parser.add_argument( '--delay_cals_file', dest='delay_cals_file', type=str, help='output file of delay calibrators', default='delay_calibrators.csv' )
-    parser.add_argument( '--subtract_file', dest='subtract_file', type=str, help='output file of sources to subtract', default='subtract_sources.csv' )
     parser.add_argument( '--match_tolerance', dest='match_tolerance', type=float, help='radius for matching LBCS to LoTSS [arcsec]', default=5. )
-    parser.add_argument( '--subtract_limit', dest='subtract_limit', type=float, help='Flux limit for sources to subtract [Jy]', default=0.5 )
     parser.add_argument( '--bright_limit_Jy', dest='bright_limit_Jy', type=float, help='Flux limit for bright sources [Jy]', default=5.0 )
-    parser.add_argument( '--image_limit_Jy', dest='image_limit_Jy', type=float, help='Flux limit for which sources to image [Jy]', default = 0.05 )
+    parser.add_argument( '--image_limit_Jy', dest='image_limit_Jy', type=float, help='Flux limit for which sources to image [Jy]', default = 0.01 )
     parser.add_argument( '--continue_no_lotss', dest='continue_no_lotss', action='store_true', help='Continue with the pipeline if no lotss cross-matches can be found?', default = False )
     parser.add_argument( 'MSname', type=str, nargs='+', help='Measurement set name (to get pointing center)' )
 
     args = parser.parse_args()
 
-    plugin_main( args.MSname, MSname=args.MSname, lotss_radius=args.lotss_radius, lbcs_radius=args.lbcs_radius, bright_limit_Jy=args.bright_limit_Jy, lotss_catalogue=args.lotss_catalogue, lbcs_catalogue=args.lbcs_catalogue, lotss_result_file=args.lotss_result_file, delay_cals_file=args.delay_cals_file, subtract_file=args.subtract_file, match_tolerance=args.match_tolerance, subtract_limit=args.subtract_limit, image_limit_Jy=args.image_limit_Jy, continue_no_lotss = str(args.continue_no_lotss) )
-
+    plugin_main( args.MSname, MSname=args.MSname, lotss_radius=args.lotss_radius, lbcs_radius=args.lbcs_radius, im_radius=args.im_radius, bright_limit_Jy=args.bright_limit_Jy, lotss_catalogue=args.lotss_catalogue, lbcs_catalogue=args.lbcs_catalogue, lotss_result_file=args.lotss_result_file, delay_cals_file=args.delay_cals_file, match_tolerance=args.match_tolerance, image_limit_Jy=args.image_limit_Jy, continue_no_lotss = str(args.continue_no_lotss) )
diff --git a/scripts/skynet.py b/scripts/skynet.py
index f51906d328da0613256ac7f56317e38109250b06..f6d102ff5e5792171f1ce1a2899fe2a46b4c8c39 100755
--- a/scripts/skynet.py
+++ b/scripts/skynet.py
@@ -4,7 +4,7 @@ import subprocess
 import shutil
 import sys
 import glob
-import regex as re
+import re
 from astropy.coordinates import SkyCoord
 from astropy.table import Table
 
diff --git a/steps/Ateamclipper.cwl b/steps/Ateamclipper.cwl
index 8ef84aafcf5f150c2cffc0bf4637dab5b58f16b4..e03c3e6455ed5703d94b4bb46ad06394c34365a7 100755
--- a/steps/Ateamclipper.cwl
+++ b/steps/Ateamclipper.cwl
@@ -1,6 +1,6 @@
 class: CommandLineTool
 cwlVersion: v1.2
-id: check_ateam_separation
+id: Ateamclipper
 label: Ateamclipper
 
 baseCommand:
@@ -36,8 +36,13 @@ hints:
     listing:
       - entry: $(inputs.msin)
         writable: true
+  - class: InplaceUpdateRequirement
+    inplaceUpdate: true
   - class: DockerRequirement
     dockerPull: vlbi-cwl
   - class: InlineJavascriptRequirement
+  - class: ResourceRequirement
+    coresMin: 12
 
 stdout: Ateamclipper.log
+stderr: Ateamclipper_err.log
diff --git a/steps/aoflagger.cwl b/steps/aoflagger.cwl
index 91146ff671b5f186cf51ee978ad120ffec31e700..42fe271b46e8112a84afb469288d775d4aca16bc 100644
--- a/steps/aoflagger.cwl
+++ b/steps/aoflagger.cwl
@@ -25,7 +25,7 @@ inputs:
       doc: Input data Column
     - id: memoryperc
       type: int?
-      default: 15
+      default: 10
       inputBinding:
         position: 0
         prefix: aoflagger.memoryperc=
@@ -72,6 +72,8 @@ requirements:
       - entry: $(inputs.msin)
         writable: true
   - class: ShellCommandRequirement
+  - class: ResourceRequirement
+    coresMin: 6
 
 hints:
   DockerRequirement:
diff --git a/steps/delay_solve.cwl b/steps/delay_solve.cwl
index 235a3bac7be8438ed4afb1071ccfcb932f5dd663..8a6c0b32ec9a9ca62ec9e794a1832f4d010a3733 100644
--- a/steps/delay_solve.cwl
+++ b/steps/delay_solve.cwl
@@ -26,6 +26,10 @@ outputs:
       type: File[]
       outputBinding:
         glob: merged_selfcal*.h5
+    - id: images
+      type: File[]
+      outputBinding:
+        glob: image*.png
     - id: logfile
       type: File[]
       outputBinding:
diff --git a/steps/download_cats.cwl b/steps/download_cats.cwl
index 6c1c91736d3b6f9a1f031d3a1ce7772019996ad5..dfbbc6f438994138ba3bfe2a0d89ad17a5a39dad 100644
--- a/steps/download_cats.cwl
+++ b/steps/download_cats.cwl
@@ -7,28 +7,28 @@ baseCommand:
     - downloadCats.py
 
 inputs:
-    - id: msin # mapfile_in = kwargs['mapfile_in']
+    - id: msin
       type: Directory[]
       inputBinding:
         position: 0
     - id: lotss_radius
       type: float?
-      default: 2.0
+      default: 1.5
     - id: lbcs_radius
       type: float?
-      default: 2.0
+      default: 1.5
+    - id: im_radius
+      type: float?
+      default: 1.24
     - id: bright_limit_Jy
       type: float
       default: 5.0
     - id: match_tolerance
       type: float?
       default: 5.0
-    - id: subtract_limit
-      type: float?
-      default: 0.5
     - id: image_limit_Jy
       type: float?
-      default: 0.05
+      default: 0.01
     - id: continue_no_lotss
       type: boolean?
       default: true
@@ -47,7 +47,7 @@ inputs:
       default: "image_catalogue.csv"
     - id: delay_catalogue_name
       type: string?
-      default: "delay_catalogue.csv"
+      default: "delay_calibrators.csv"
     - id: subtract_catalogue_name
       type: string?
       default: "subtract_sources.csv"
@@ -65,18 +65,18 @@ outputs:
     #  type: File
     #  outputBinding:
     #    glob: $(inputs.lbcs_skymodel_name)
-    - id: best_delay_catalogue
-      type: File
-      outputBinding:
-        glob: best_delay_*.csv
-    #- id: image_catalogue
+    #- id: best_delay_catalogue
     #  type: File
     #  outputBinding:
-    #    glob: $(inputs.image_catalogue_name)
-    #- id: delay_catalogue
+    #    glob: best_delay_*.csv
+    #- id: image_catalogue
     #  type: File
     #  outputBinding:
-    #    glob: $(inputs.delay_catalogue_name)
+    #    glob: $(inputs.image_catalogue_name)
+    - id: delay_catalogue
+      type: File
+      outputBinding:
+        glob: $(inputs.delay_catalogue_name)
     #- id: subtract_catalogue
     #  type: File
     #  outputBinding:
diff --git a/steps/dp3_concat.cwl b/steps/dp3_concat.cwl
index 850d0b9e721f935306ac6b0a2029b253dd76f9e5..601b4eeee836e08bc72592e8e3f8317eb39fb1fa 100644
--- a/steps/dp3_concat.cwl
+++ b/steps/dp3_concat.cwl
@@ -89,6 +89,8 @@ hints:
     listing:
       - entry: $(inputs.msin)
         writable: true
+  - class: ResourceRequirement
+    coresMin: 6
 
 stdout: dp3_concat.log
 stderr: dp3_concat_err.log
diff --git a/steps/dp3_prep_target.cwl b/steps/dp3_prep_target.cwl
index 401ef3790e67c21f80cef8e0f3597196b7726dea..54d05253407178659e423e10ea4d3262c3a661f6 100644
--- a/steps/dp3_prep_target.cwl
+++ b/steps/dp3_prep_target.cwl
@@ -70,6 +70,8 @@ arguments:
 
 requirements:
   - class: InlineJavascriptRequirement
+  - class: ResourceRequirement
+    coresMin: 6
 #  - class: InitialWorkDirRequirement
 #    listing:
 #      - entry: $(inputs.msin)
diff --git a/steps/predict.cwl b/steps/predict.cwl
index 799e2535af8a7302cca636761d17331585008e6c..1859008f6adb3857f706c89545ca3e4eeca8268e 100644
--- a/steps/predict.cwl
+++ b/steps/predict.cwl
@@ -109,6 +109,8 @@ outputs:
 hints:
   - class: DockerRequirement
     dockerPull: vlbi-cwl:latest
+  - class: ResourceRequirement
+    coresMin: 6
 
 stdout: predict_ateam.log
 stderr: predict_ateam_err.log
diff --git a/workflows/delay-calibration.cwl b/workflows/delay-calibration.cwl
index ec1a84d4c420e1b15ae9ddfead6eb44de93ff7e1..3ad15abe282383211ec928443565f40ca17be89f 100644
--- a/workflows/delay-calibration.cwl
+++ b/workflows/delay-calibration.cwl
@@ -50,7 +50,7 @@ steps:
           source: phasesol
       out:
         - id: parset
-        - id: best_delay_cats
+        - id: delay_calibrators
         - id: logdir
         - id: msout
         - id: initial_flags
@@ -85,7 +85,7 @@ steps:
         - id: msin
           source: sort-concatenate-flag/msout
         - id: delay_calibrator
-          source: setup/best_delay_cats
+          source: setup/delay_calibrators
         - id: configfile
           source: configfile
         - id: selfcal
@@ -102,6 +102,7 @@ steps:
       out:
         - id: msout
         - id: solutions
+        - id: pictures
         - id: phaseup_flags
         - id: logdir
 #        - id: summary_file
@@ -129,13 +130,21 @@ outputs:
     type: Directory[]
 
   - id: delay_cat
-    outputSource: setup/best_delay_cats
+    outputSource: setup/delay_calibrators
     type: File
 
   - id: logs
     outputSource: store_logs/dir
     type: Directory
 
+  - id: pictures
+    outputSource: phaseup/pictures
+    type: File[]
+
+  - id: solutions
+    outputSource: phaseup/solutions
+    type: File[]
+
 #  - id: summary_file
 #    outputSource: phaseup-concat.cwl
 #    type: File
diff --git a/workflows/phaseup-concat.cwl b/workflows/phaseup-concat.cwl
index 42ac0d1ac837f85a4f953e859caf7c9e1b69f49c..3f4d0315bcd8d9e00a16a33a659f061fe5e12eba 100644
--- a/workflows/phaseup-concat.cwl
+++ b/workflows/phaseup-concat.cwl
@@ -204,8 +204,6 @@ steps:
     in:
       - id: msin
         source: delay_cal_model/msout
-        #source: phaseup_concatenate/msout
-        #valueFrom: $(self[0])
       - id: configfile
         source: configfile
       - id: selfcal
@@ -214,6 +212,7 @@ steps:
         source: h5merger
     out:
       - id: h5parm
+      - id: images
       - id: logfile
     run: ../steps/delay_solve.cwl
     label: delay_solve
@@ -262,7 +261,7 @@ steps:
         linkMerge: merge_flattened
         source:
           - prep_delay/logfile
-          - concat_logfiles_phaseup/output 
+          - concat_logfiles_phaseup/output
           - sort_concatenate/logfile
           - concat_logfiles_phaseup/output
           - delay_cal_model/logfile
@@ -288,6 +287,9 @@ outputs:
   - id: logdir
     outputSource: save_logfiles/dir
     type: Directory
+  - id: pictures
+    type: File[]
+    outputSource: delay_solve/images
 #  - id: summary_file
 #    type: File
 #    outputSource: summary/summary_file
diff --git a/workflows/setup.cwl b/workflows/setup.cwl
index 2c99b4432bc7d60023172b1bd64552af472a94eb..5d0df7bdf40e55aae5622f1f4ef83c1dbf248cd8 100644
--- a/workflows/setup.cwl
+++ b/workflows/setup.cwl
@@ -46,7 +46,7 @@ steps:
         - id: msin
           source: msin
       out:
-        - id: best_delay_catalogue
+        - id: delay_catalogue
         - id: logfile
       run: ../steps/download_cats.cwl
       label: download_cats
@@ -154,8 +154,8 @@ outputs:
     - id: logdir
       outputSource: save_logfiles/dir
       type: Directory
-    - id: best_delay_cats
-      outputSource: download_cats/best_delay_catalogue
+    - id: delay_calibrators
+      outputSource: download_cats/delay_catalogue
       type: File
     - id: parset
       outputSource: dp3_make_parset/parset