diff --git a/CHANGELOG.md b/CHANGELOG.md
index ea5d18c..162ecbc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,6 +40,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Use spectral-cube common beam calculations for common beams (#385).
 - Mosaic weights are now generated after convolution and regridding, and will undo any primary beam corrections (#396).
 - Don't skip making large mosaics (but keep the log messages to inform users this is happening) (#405).
+- Begin removal of analysisUtils (#402).
 
 ### Fixed
 
diff --git a/phangsPipeline/casaCubeRoutines.py b/phangsPipeline/casaCubeRoutines.py
index bd99f3a..37b41ee 100644
--- a/phangsPipeline/casaCubeRoutines.py
+++ b/phangsPipeline/casaCubeRoutines.py
@@ -7,7 +7,6 @@
 import logging
 import os
 
-import analysisUtils as au
 import numpy as np
 import scipy.ndimage as nd
 from astropy.io import fits
@@ -45,7 +44,7 @@ def check_getchunk_putchunk_memory_issue(
 
     Args:
         infile: input image file.
-        myia: optional, if not given then we will run `myia = au.createCasaTool(casaStuff.iatool)` then `myia.open(infile)`.
+        myia: optional, if not given then we will run `myia = casaStuff.iatool()` then `myia.open(infile)`.
         return_data: if True then we will return the data if `getchunk()` runs properly withut a memory issue.
         return_mask: if True then we will return the mask if `getchunk(getmask=True)` runs properly withut a memory issue.
         return_shape: if True then we will return the shape of the cube data.
@@ -64,7 +63,7 @@ def check_getchunk_putchunk_memory_issue(
     cube_mask = None
     has_opened_file = False
     if myia is None:
-        myia = au.createCasaTool(casaStuff.iatool)
+        myia = casaStuff.iatool()
         myia.open(infile)
         has_opened_file = True
     if not myia.isopen():
@@ -177,14 +176,14 @@ def get_mask(infile, huge_cube_workaround=True):
     #    os.system('rm -rf ' + infile + '.temp.fits')
     #
     #else:
-    #    myia = au.createCasaTool(casaStuff.iatool)
+    #    myia = casaStuff.iatool()
     #    myia.open(infile+'.temp_deg')
     #    mask = myia.getchunk(getmask=True)
     #    myia.close()
     #
     #os.system('rm -rf ' + infile + '.temp_deg')
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     myia.open(infile)
 
@@ -237,7 +236,7 @@ def copy_mask(infile, outfile, huge_cube_workaround=True):
 
     # use putregion to update pixel mask
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     myia.open(outfile)
 
@@ -259,7 +258,7 @@ def copy_mask(infile, outfile, huge_cube_workaround=True):
     else: # getchunk was unsuccessful, has memory issue
         # putregion channel by channel
         logger.debug('putregion channel by channel for known memory issue')
-        myrg = au.createCasaTool(casaStuff.rgtool)
+        myrg = casaStuff.rgtool()
         nx = out_shape[0]
         ny = out_shape[1]
         #logger.debug('copy_mask mask.shape {}'.format(mask.shape))
@@ -324,13 +323,13 @@ def multiply_cube_by_value(infile, value, brightness_unit, huge_cube_workaround=
     #                         overwrite=True)
     #    os.system('rm -rf ' + infile + '.fits')
     #else:
-    #    myia = au.createCasaTool(casaStuff.iatool)
+    #    myia = casaStuff.iatool()
     #    myia.open(infile)
     #    vals = myia.getchunk()
     #    vals *= value
     #    myia.putchunk(vals)
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     myia.open(infile)
 
@@ -627,7 +626,7 @@ def trim_rind(
         target_file = infile
 
     # Figure out the extent of the image inside the cube
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(target_file)
     mask = myia.getchunk(getmask=True)
     elt = nd.generate_binary_structure(2,1)
diff --git a/phangsPipeline/casaFeatherRoutines.py b/phangsPipeline/casaFeatherRoutines.py
index 5b79feb..d804227 100644
--- a/phangsPipeline/casaFeatherRoutines.py
+++ b/phangsPipeline/casaFeatherRoutines.py
@@ -6,7 +6,6 @@
 import logging
 import os
 
-import analysisUtils as au
 import numpy as np
 
 from . import casaCubeRoutines as ccr
@@ -115,12 +114,12 @@ def prep_sd_for_feather(
             overwrite=overwrite)
 
     # Get out and match axis orders from interferometric images
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(interf_file)
     interf_cs = myia.coordsys()
     myia.close()
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(current_outfile)
     sd_cs = myia.coordsys()
     myia.close()
@@ -235,7 +234,7 @@ def feather_two_cubes(
         os.system('cp -rf '+interf_file+' '+current_interf_file)
         os.system('cp -rf '+sd_file+' '+current_sd_file)
 
-        myia = au.createCasaTool(casaStuff.iatool)
+        myia = casaStuff.iatool()
 
         # As noted in [https://casa.nrao.edu/docs/casaref/image.putchunk.html], 
         # "If all the pixels didn't easily fit in memory, you would iterate through 
@@ -390,7 +389,7 @@ def feather_two_cubes(
 
     # Finally, if we don't have a 4th (Stokes I) axis, add this
     # here to avoid feather crashing
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     myia.open(current_interf_file)
     interf_shape = myia.shape()
diff --git a/phangsPipeline/casaImagingRoutines.py b/phangsPipeline/casaImagingRoutines.py
index 00cafda..06bf37f 100644
--- a/phangsPipeline/casaImagingRoutines.py
+++ b/phangsPipeline/casaImagingRoutines.py
@@ -7,11 +7,11 @@
 import logging
 import os
 
-import analysisUtils as au
 import numpy as np
 
 from . import casaMaskingRoutines as cmr
 from . import casaStuff
+from .casaVisRoutines import get_spw_info, pick_cell_and_im_size
 from .clean_call import CleanCall
 
 logger = logging.getLogger(__name__)
@@ -50,65 +50,21 @@ def estimate_cell_and_imsize(
     else:
         pblevel = clean_call.get_param('pblimit')
 
-    # These are the CASA-preferred sizes for fast FFTs
-
-    valid_sizes = []
-    for ii in range(10):
-        for kk in range(3):
-            for jj in range(3):
-                valid_sizes.append(2 ** (ii + 1) * 5 ** jj * 3 ** kk)
-    valid_sizes = sorted(valid_sizes)
-    valid_sizes = np.array(valid_sizes)
-
-    # Cell size implied by baseline distribution from analysis
-    # utilities.
-
-    au_cellsize, au_imsize, _ = au.pickCellSize(infile,
-                                                imsize=True,
-                                                npix=oversamp,
-                                                intent='',
-                                                pblevel=pblevel,
-                                                )
-    xextent = au_cellsize * au_imsize[0] * 1.2
-    yextent = au_cellsize * au_imsize[1] * 1.2
-
-    # Make the cell size a nice round number
-
-    if au_cellsize < 0.01:
-        cell_size = au_cellsize
-    elif 0.01 <= au_cellsize < 0.05:
-        cell_size = np.floor(au_cellsize / 0.005) * 0.005
-    elif 0.05 <= au_cellsize < 0.1:
-        cell_size = np.floor(au_cellsize / 0.01) * 0.01
-    elif 0.1 <= au_cellsize < 0.5:
-        cell_size = np.floor(au_cellsize / 0.05) * 0.05
-    elif 0.5 <= au_cellsize < 1.0:
-        cell_size = np.floor(au_cellsize / 0.1) * 0.1
-    elif 1.0 <= au_cellsize < 2.0:
-        cell_size = np.floor(au_cellsize / 0.25) * 0.25
-    elif 2.0 <= au_cellsize < 5.0:
-        cell_size = np.floor(au_cellsize / 0.5) * 0.5
-    else:
-        cell_size = np.floor(au_cellsize / 1.0) * 0.5
-
-    # Now make the image size a good number for the FFT
-
-    need_cells_x = xextent / cell_size
-    need_cells_y = yextent / cell_size
-
-    cells_x = np.min(valid_sizes[valid_sizes > need_cells_x])
-    cells_y = np.min(valid_sizes[valid_sizes > need_cells_y])
+    # Cell size implied by baseline distribution
+    cell_size, image_size = pick_cell_and_im_size(
+        infile,
+        npix=oversamp,
+        pblevel=pblevel,
+    )
 
     # If requested, force the mosaic to be square. This avoids
     # pathologies in CASA versions 5.1 and 5.3.
-
     if force_square:
-        if cells_y < cells_x:
-            cells_y = cells_x
-        if cells_x < cells_y:
-            cells_x = cells_y
+        if image_size[0] < image_size[1]:
+            image_size[1] = image_size[0]
+        if image_size[0] > image_size[1]:
+            image_size[0] = image_size[1]
 
-    image_size = [int(cells_x), int(cells_y)]
     cell_size_string = str(cell_size) + 'arcsec'
 
     return cell_size_string, image_size
@@ -786,18 +742,17 @@ def clean_loop(
     working_call.set_param('calcpsf', False)
 
     # Note the number of channels, which is used in setting the number
-    # of iterations that we give to an individual clean call.
-
-    vm = au.ValueMapping(working_call.get_param('vis'))
-
-    # If we have a number of channels defined, then use that
+    # of iterations that we give to an individual clean call. If we
+    # have a number of channels defined, then use that
     nchan = working_call.get_param('nchan')
 
     # Otherwise, pull the total number of channels from the spwInfo
     if nchan is None:
         nchan = -1
     if nchan <= 0:
-        nchan = vm.spwInfo[0]['numChannels']
+        vis = working_call.get_param('vis')
+        spw_info = get_spw_info(vis)
+        nchan = spw_info[0]['numChannels']
 
     # Create a text record of progress through successive clean calls.
 
@@ -1016,7 +971,7 @@ def clean_loop(
             pb_image = (_imagename + '.pb') if _imagename else None
             if pb_image and os.path.isdir(pb_image):
                 try:
-                    myia = au.createCasaTool(casaStuff.iatool)
+                    myia = casaStuff.iatool()
                     myia.open(pb_image)
                     pb_data = myia.getchunk()
                     myia.close()
@@ -1115,7 +1070,7 @@ def calc_residual_statistics(
         logger.error('Error! The input file "' + mask_name + '" was not found!')
         return
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     myia.open(mask_name)
     mask = myia.getchunk()
diff --git a/phangsPipeline/casaLegacySingleDishRoutines.py b/phangsPipeline/casaLegacySingleDishRoutines.py
index 53d2548..60f3711 100644
--- a/phangsPipeline/casaLegacySingleDishRoutines.py
+++ b/phangsPipeline/casaLegacySingleDishRoutines.py
@@ -49,7 +49,13 @@
 import numpy as np
 
 from . import casaStuff
-from .utilsSingleDish import getTPSampling, get_first_arr_val
+from .casaVisRoutines import get_science_spws
+from .utilsSingleDish import (
+    getTPSampling,
+    get_first_arr_val,
+    get_sourcename,
+    read_source_coordinates,
+)
 
 es = au.stuffForScienceDataReduction()
 
@@ -86,19 +92,11 @@ def checkpipeline():
         logger.info("> Setting the variable 'pipeline' to False")
         return False
 
-# Creating CASA tools
-#def createCasaTool(mytool):
-#
-#    if (type(casac.Quantity) != type):  # casa 4.x
-#        myt = mytool()
-#    else:  # casa 3.x
-#        myt = mytool.create()
-#    return(myt)
 
 # Retrieve name of the column
 def getDataColumnName(inputms):
 
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(inputms)
     colnames = mytb.colnames()
     if 'FLOAT_DATA' in colnames:
@@ -112,7 +110,7 @@ def getDataColumnForSDBaseline(vis):
     """
     Returns the names of the corrected data columns (corrected) in a measurement set.
     """
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(vis)
     names = copy.copy(mytb.colnames())
     mytb.close()
@@ -178,8 +176,8 @@ def scaleAutocorr(vis, scale=1., antenna='', spw='', field='', scan=''):
         logger.warning("No table.dat. This does not appear to be an MS.")
         return
 
-    mymsmd = au.createCasaTool(casaStuff.msmdtool)
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mymsmd = casaStuff.msmdtool()
+    mytb = casaStuff.tbtool()
 
     conditions = ["ANTENNA1==ANTENNA2"]
 
@@ -242,7 +240,7 @@ def scaleAutocorr(vis, scale=1., antenna='', spw='', field='', scan=''):
 # Create vector with antenna names
 def read_ants_names(filename):
 
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(filename + '/ANTENNA')
     vec_ants = mytb.getcol('NAME')
     mytb.close()
@@ -266,7 +264,7 @@ def get_tsysmap(tsysmap,spws_scie,spws_tsys,freq_rep_scie,freq_rep_tsys):
 def read_spw(filename,source):
 
     # Tsys spws (index)
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(filename + '/SYSCAL')
 
     spwstsys  = mytb.getcol('SPECTRAL_WINDOW_ID')
@@ -285,7 +283,7 @@ def read_spw(filename,source):
     names          = mytb.getcol('NAME')
     rest_freq_scie = [rest_freq_scie[i] for i in range(len(spws_scie)) if "FULL_RES" in names[spws_scie[i]]]
     spws_scie      = [spw for spw in spws_scie if "FULL_RES" in names[spw]]
-    spws_scie      = au.getScienceSpws(filename)
+    spws_scie      = get_science_spws(filename)
     spws_scie      = spws_scie.split(",")
     spws_scie = [int(i) for i in spws_scie]
 
@@ -315,7 +313,7 @@ def read_spw(filename,source):
 # Get information of the source velocity
 def read_vel_source(filename,source):
 
-    mytb  = au.createCasaTool(casaStuff.tbtool)
+    mytb  = casaStuff.tbtool()
     mytb.open(filename + '/SOURCE')
     names = mytb.getcol('NAME')
     numli = mytb.getcol('NUM_LINES')
@@ -456,7 +454,7 @@ def convert_vel2chan_line(filename_in,freq_rest,vel_line,spw_line,coords,date):
     freq2_topo = au.lsrkToTopo(freq2, date, ra, dec)
 
     if fsuffix == '.asap':
-        mytb  = au.createCasaTool(casaStuff.tbtool)
+        mytb  = casaStuff.tbtool()
         mytb.open(filename_in)
         nchan = mytb.getkeyword('nChan')
         if_eq = mytb.getcol('FREQ_ID',startrow=1,nrow=1)
@@ -470,7 +468,7 @@ def convert_vel2chan_line(filename_in,freq_rest,vel_line,spw_line,coords,date):
         mytb.close()
 
     else:
-        mytb  = au.createCasaTool(casaStuff.tbtool)
+        mytb  = casaStuff.tbtool()
         mytb.open(filename_in+os.sep+'SPECTRAL_WINDOW')
         spwids = np.arange(mytb.nrows())
         spw_line = -1
@@ -581,73 +579,6 @@ def extract_jyperk(filename, spw_line, pipeline):
         jyperk = kw['jyperk']
         return jyperk
 
-# Read source coordinates
-def read_source_coordinates(filename,source):
-
-    coord_source = au.getRADecForSource(filename,source)
-    RA_h  = (coord_source.split(' ')[0]).split(':')[0]
-    RA_m  = (coord_source.split(' ')[0]).split(':')[1]
-    RA_s  = (coord_source.split(' ')[0]).split(':')[2]
-    DEC_d = (coord_source.split(' ')[1]).split(':')[0]
-    DEC_m = (coord_source.split(' ')[1]).split(':')[1]
-    DEC_s = (coord_source.split(' ')[1]).split(':')[2]
-    coord = "J2000  "+str(RA_h)+"h"+str(RA_m)+"m"+str(RA_s[0:6])+" "+str(DEC_d)+"d"+str(DEC_m)+"m"+str(DEC_s)
-    return coord
-
-# Get source name
-def get_sourcename(filename,
-                   source='all',
-                   ):
-
-    mytb   = au.createCasaTool(casaStuff.msmdtool)
-    mytb.open(filename)
-
-    # Get source field names
-    fieldnames = mytb.fieldnames()
-
-    # If we're looking for all intents, then just take the (first) science target
-    if source == 'all':
-        field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')[0]
-
-    # Else, key in on the source name. Do this in a regex-style
-    else:
-
-        # For field selection, this can be done with forward slashes,
-        # so make sure to strip them
-        if source.startswith("/"):
-            source = source.strip("/")
-
-        # Get science fieldnames as names
-        fieldnames_asname = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE', asnames=True)
-
-        regex_matches = [re.fullmatch(source, x) for x in fieldnames_asname]
-        regex_idxs = [x is not None for x in regex_matches]
-        found_field = any(regex_idxs)
-
-        # If we've found a match, update the source name with the first
-        # fieldname
-        if found_field:
-            source = fieldnames_asname[regex_idxs][0]
-
-        if not found_field:
-            logger.warning(f"source {source} not in field names: {', '.join(fieldnames)}")
-            return None
-
-        field = mytb.fieldsforname(source)
-
-        # If we don't find anything, get the science targets
-        if len(field) == 0:
-            field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')
-
-        # Take the first result
-        field = field[0]
-
-    source_name = fieldnames[field]
-
-    mytb.close()
-
-    return source_name
-
 # Create string of spws to apply the Tsys
 def str_spw_apply_tsys(spws_info):
     #science spws
@@ -1491,7 +1422,7 @@ def imaging(source, name_line, phcenter, vel_source, source_vel_kms, vwidth_kms,
         infile = ms_concat_filename
 
     # Read frequency
-    mymsmd = au.createCasaTool(casaStuff.msmdtool)
+    mymsmd = casaStuff.msmdtool()
     mymsmd.open(infile)
     freq = mymsmd.meanfreq(0)
     mymsmd.close()
@@ -1585,7 +1516,7 @@ def imaging(source, name_line, phcenter, vel_source, source_vel_kms, vwidth_kms,
     #ia.open('ALMA_TP.'+source+'.'+name_line+suffix+'.image')
     #ia.setrestoringbeam(major = str(sfbeam)+'arcsec', minor = str(sfbeam)+'arcsec', pa = '0deg')
     #ia.done()
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(outimage) # 'ALMA_TP.'+source+'.'+name_line+suffix+'.image'
     myia.setrestoringbeam(major = str(sfbeam)+'arcsec', minor = str(sfbeam)+'arcsec', pa = '0deg')
     mask = myia.getchunk(getmask=True)
@@ -1933,4 +1864,3 @@ def run_ALMA_TP_tools(
     os.chdir(ori_path)
 
 
-
diff --git a/phangsPipeline/casaMaskingRoutines.py b/phangsPipeline/casaMaskingRoutines.py
index d59c996..828ad1d 100644
--- a/phangsPipeline/casaMaskingRoutines.py
+++ b/phangsPipeline/casaMaskingRoutines.py
@@ -6,7 +6,6 @@
 import logging
 import os
 
-import analysisUtils as au
 import numpy as np
 import scipy.ndimage as ndimage
 from astropy.io import fits
@@ -162,7 +161,7 @@ def noise_for_cube(
             logger.error('maskfile specified but not found - ' + maskfile)
             return (None)
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(infile)
 
     has_memory_issue, cube_shape = ccr.check_getchunk_putchunk_memory_issue(
@@ -191,7 +190,7 @@ def noise_for_cube(
 
         myia_mask = None
         if maskfile is not None:
-            myia_mask = au.createCasaTool(casaStuff.iatool)
+            myia_mask = casaStuff.iatool()
             myia_mask.open(maskfile)
 
         per_channel_noise = []
@@ -303,7 +302,7 @@ def read_cube(infile, huge_cube_workaround=True):
         # Remove intermediate fits file
         os.system('rm -rf ' + infile + '.fits')
     else:
-        myia = au.createCasaTool(casaStuff.iatool)
+        myia = casaStuff.iatool()
         myia.open(infile)
         cube = myia.getchunk()
         myia.close()
@@ -344,7 +343,7 @@ def write_mask(infile, outfile, mask, huge_cube_workaround=True):
         # Remove the intermediate fits file
         os.system('rm -rf ' + outfile + '.fits')
     else:
-        myia = au.createCasaTool(casaStuff.iatool)
+        myia = casaStuff.iatool()
         myia.open(outfile)
         myia.putchunk(mask)
         myia.close()
@@ -478,7 +477,7 @@ def apply_additional_mask(
     beam based masks by setting the PB file to new_mask_file and the
     pb_limit as new_thresh.
     """
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(new_mask_file)
     new_mask = myia.getchunk()
     myia.close()
@@ -515,8 +514,8 @@ def import_and_align_mask(
                          overwrite=True)
 
     # Prepare analysis utility tool
-    myia = au.createCasaTool(casaStuff.iatool)
-    myim = au.createCasaTool(casaStuff.imtool)
+    myia = casaStuff.iatool()
+    myim = casaStuff.imtool()
 
     # Read mask data
     # myia.open(out_file+'.temp_copy')
diff --git a/phangsPipeline/casaMosaicRoutines.py b/phangsPipeline/casaMosaicRoutines.py
index 2d90f2e..d103e3b 100644
--- a/phangsPipeline/casaMosaicRoutines.py
+++ b/phangsPipeline/casaMosaicRoutines.py
@@ -7,7 +7,6 @@
 import logging
 import os
 
-import analysisUtils as au
 import numpy as np
 
 from . import casaMaskingRoutines as cma
@@ -760,7 +759,7 @@ def generate_weight_file(
     # Copy the template and open it
     os.system("cp -r "+template+" "+outfile)
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(outfile)
 
     # Check for memory issues
@@ -1004,7 +1003,7 @@ def mosaic_aligned_data(
 
     # We might need to add a 4th (Stokes) axis
     # to prevent crashing
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
 
     for this_infile in infile_list:
 
@@ -1142,8 +1141,8 @@ def mosaic_aligned_data(
         os.system("cp -r "+local_imlist[0]+" "+sum_file)
         os.system("cp -r "+local_imlist[0]+" "+weight_file)
         
-        myia_sum = au.createCasaTool(casaStuff.iatool)
-        myia_weight = au.createCasaTool(casaStuff.iatool)
+        myia_sum = casaStuff.iatool()
+        myia_weight = casaStuff.iatool()
         myia_sum.open(sum_file)
         myia_weight.open(weight_file)
         
@@ -1171,7 +1170,7 @@ def mosaic_aligned_data(
                     casaStuff.immath(imagename=temp_chan_images, mode='evalexpr', 
                                     expr=lel_exp_sum, outfile=temp_sum)
                     
-                    myia_temp = au.createCasaTool(casaStuff.iatool)
+                    myia_temp = casaStuff.iatool()
                     myia_temp.open(temp_sum)
                     sum_slice = myia_temp.getchunk()
                     myia_temp.close()
@@ -1237,7 +1236,7 @@ def mosaic_aligned_data(
                         casaStuff.immath(imagename=temp_chan_images, mode='evalexpr',
                                         expr=lel_exp_sum, outfile=temp_sum)
                         
-                        myia_temp = au.createCasaTool(casaStuff.iatool)
+                        myia_temp = casaStuff.iatool()
                         myia_temp.open(temp_sum)
                         sum_slice = myia_temp.getchunk()
                         myia_temp.close()
@@ -1267,7 +1266,7 @@ def mosaic_aligned_data(
 
     # Just to be safe, reset the masks on the two images.
 
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(sum_file)
     myia.set(pixelmask=1)
     myia.close()
@@ -1298,8 +1297,8 @@ def mosaic_aligned_data(
         os.system("cp -r "+sum_file+" "+temp_file)
         os.system("cp -r "+weight_file+" "+local_maskfile)
         
-        myia_temp = au.createCasaTool(casaStuff.iatool)
-        myia_mask = au.createCasaTool(casaStuff.iatool)
+        myia_temp = casaStuff.iatool()
+        myia_mask = casaStuff.iatool()
         myia_temp.open(temp_file)
         myia_mask.open(local_maskfile)
 
@@ -1368,7 +1367,7 @@ def mosaic_aligned_data(
             myia_mask.close()
 
     # Add in potential Stokes axis to the mask file
-    myia = au.createCasaTool(casaStuff.iatool)
+    myia = casaStuff.iatool()
     myia.open(local_maskfile)
     in_shape = myia.shape()
     if len(in_shape) != 4:
diff --git a/phangsPipeline/casaRoutineTests.py b/phangsPipeline/casaRoutineTests.py
index 0f149c1..8cfd67e 100644
--- a/phangsPipeline/casaRoutineTests.py
+++ b/phangsPipeline/casaRoutineTests.py
@@ -8,7 +8,6 @@
 
 import logging
 
-import analysisUtils as au
 import numpy as np
 
 from . import casaMaskingRoutines as cma
diff --git a/phangsPipeline/casaSingleDishALMAWrapper.py b/phangsPipeline/casaSingleDishALMAWrapper.py
index d5ea31e..3a2f40b 100644
--- a/phangsPipeline/casaSingleDishALMAWrapper.py
+++ b/phangsPipeline/casaSingleDishALMAWrapper.py
@@ -8,9 +8,8 @@
 import numpy as np
 from astropy.table import Table
 
-from . import casaLegacySingleDishRoutines as csdr
 from . import casaStuff
-from .utilsSingleDish import getTPSampling
+from .utilsSingleDish import getTPSampling, get_sourcename, read_source_coordinates
 
 logger = logging.getLogger(__name__)
 
@@ -87,7 +86,7 @@ def SDImaging(filename,
     imsize = int(round(maxsize/cell)*1.5)
 
     if phcenter == '':
-        coord_phase = csdr.read_source_coordinates(filename, source)
+        coord_phase = read_source_coordinates(filename, source)
     else:
         coord_phase = phcenter
 
@@ -374,7 +373,7 @@ def runALMAPipeline(path_galaxy,
     this_vis = ms_concat_filename
 
     # read the source name directly from the ms
-    source = csdr.get_sourcename(this_vis, source=in_source)
+    source = get_sourcename(this_vis, source=in_source)
 
     for this_product in product_dict:
 
diff --git a/phangsPipeline/casaVisRoutines.py b/phangsPipeline/casaVisRoutines.py
index 744e5e9..9242f8e 100644
--- a/phangsPipeline/casaVisRoutines.py
+++ b/phangsPipeline/casaVisRoutines.py
@@ -7,7 +7,6 @@
 import os
 import shutil
 
-import analysisUtils as au
 import astropy.constants as const
 import astropy.units as u
 import numpy as np
@@ -234,7 +233,7 @@ def split_science_targets(
     # Verify the column to use. If present, we use the corrected
     # column. If not, then we use the data column.
 
-    mytb = au.createCasaTool(casaStuff.tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(infile, nomodify = True)
     colnames = mytb.colnames()
     if 'CORRECTED_DATA' in colnames:
@@ -571,14 +570,13 @@ def find_spws_for_line(
     spw_lowest_ghz = None
     spw_highest_ghz = None
 
-    logger.debug("... vm = au.ValueMapping(infile) ...")
-    vm = au.ValueMapping(infile)
-    logger.debug("... vm = au.ValueMapping(infile) done")
+    spw_info = get_spw_info(infile)
+    scans_for_spw = get_scans_for_spw(infile)
 
-    for this_spw in vm.spwInfo.keys():
+    for this_spw in spw_info.keys():
 
-        spw_high_ghz = np.max(vm.spwInfo[this_spw]['edgeChannels'])/1e9
-        spw_low_ghz = np.min(vm.spwInfo[this_spw]['edgeChannels'])/1e9
+        spw_high_ghz = np.max(spw_info[this_spw]['edgeChannels'])/1e9
+        spw_low_ghz = np.min(spw_info[this_spw]['edgeChannels'])/1e9
         logger.debug(
             "... spw: %s, freq: %.6f - %.6f GHz" %
             (this_spw, spw_low_ghz, spw_high_ghz))
@@ -590,12 +588,12 @@ def find_spws_for_line(
             continue
 
         if max_chanwidth_ghz is not None:
-            spw_chanwidth_ghz = abs(vm.spwInfo[this_spw]['chanWidth'])/1e9
+            spw_chanwidth_ghz = abs(spw_info[this_spw]['chanWidth'])/1e9
             if spw_chanwidth_ghz > max_chanwidth_ghz:
                 continue
 
         if require_data:
-            if len(vm.scansForSpw[this_spw]) == 0:
+            if len(scans_for_spw[this_spw]) == 0:
                 continue
 
         if require_full_line_coverage and not (
@@ -647,8 +645,7 @@ def find_spws_for_science(
         infile=None, require_data=False, exit_on_error=True, as_list=False):
     """
     List all spectral windows that we judge likely to be used for
-    science. Mostly wraps analysisUtils rather than reinventing the
-    wheel.
+    science.
     """
 
     # Check inputs
@@ -667,25 +664,26 @@ def find_spws_for_science(
             'Error! The input uv data measurement set "'+infile +
             '"does not exist!')
 
-    # Call the analysisUtil version.
-
-    spw_string = au.getScienceSpws(
-        infile, intent='OBSERVE_TARGET*')
+    # Get science SPWs
+    spw_string = get_science_spws(
+        vis=infile,
+        intent='OBSERVE_TARGET*',
+    )
     if spw_string is None or len(spw_string) == 0:
-        spw_string = au.getScienceSpws(
-            infile, intent='OBSERVE_TARGET#ON_SOURCE')
+        spw_string = get_science_spws(
+            vis=infile,
+            intent='OBSERVE_TARGET#ON_SOURCE',
+        )
 
     spw_list = []
     for this_spw_string in spw_string.split(','):
         spw_list.append(int(this_spw_string))
 
-    # Shouldn't get here, I think, because of the analysisUtils logic
-
     if require_data:
-        vm = au.ValueMapping(infile)
+        scans_for_spw = get_scans_for_spw(infile)
 
         for spw in spw_list:
-            if len(vm.scansForSpw[spw]) == 0:
+            if len(scans_for_spw[spw]) == 0:
                 spw_list.remove(spw)
     # Return
 
@@ -740,18 +738,18 @@ def spw_string_for_freq_ranges(
     if not isinstance(just_spw, list):
         just_spw = [just_spw]
 
-    vm = au.ValueMapping(infile)
+    spw_info = get_spw_info(infile)
 
     # Loop over spectral windows
     spw_flagging_string = ''
     first_string = True
-    for this_spw in vm.spwInfo:
+    for this_spw in spw_info:
 
         if len(just_spw) > 0:
             if this_spw not in just_spw:
                 continue
 
-        freq_axis = vm.spwInfo[this_spw]['chanFreqs']
+        freq_axis = spw_info[this_spw]['chanFreqs']
         half_chan = abs(freq_axis[1]-freq_axis[0])*0.5
         chan_axis = np.arange(len(freq_axis))
         mask_axis = np.zeros_like(chan_axis, dtype='bool')
@@ -843,7 +841,7 @@ def compute_common_chanwidth(
             this_infile, line, vsys_kms=vsys_kms, vwidth_kms=vwidth_kms,
             require_full_line_coverage=require_full_line_coverage)
 
-        chan_widths_hz = au.getChanWidths(this_infile, spw_list_string)
+        chan_widths_hz = get_chan_widths(this_infile, spw_list_string)
 
         # Convert to km/s and return
         for this_chan_width_hz in chan_widths_hz:
@@ -951,10 +949,10 @@ def batch_extract_line(
 
                 # This zaps the whole table:
                 if os.path.exists(this_outfile+os.sep+'POINTING'):
-                    au.clearPointingTable(this_outfile)
+                    clear_pointing_table(this_outfile)
                 else:
                     copy_pointing = False
-                    #logger.debug('Warning! Failed to run au.clearPointingTable(%r)'%(this_outfile))
+                    #logger.debug('Warning! Failed to run clear_pointing_table(%r)'%(this_outfile))
 
     # Allow a small tolerance in the channel width
     if allow_freqtol_chanfrac:
@@ -1045,7 +1043,7 @@ def suggest_extraction_scheme(
 
     for this_infile in infile_list:
 
-        vm = au.ValueMapping(this_infile)
+        spw_info = get_spw_info(this_infile)
 
         spw_list = find_spws_for_line(
             this_infile, restfreq_ghz=restfreq_ghz,
@@ -1058,7 +1056,7 @@ def suggest_extraction_scheme(
 
         for this_spw in spw_list:
 
-            chan_width_ghz = np.abs(vm.spwInfo[this_spw]['chanWidth'])/1e9
+            chan_width_ghz = np.abs(spw_info[this_spw]['chanWidth'])/1e9
 
             # Using RADIO convention:
             chan_width_kms = chan_width_ghz/restfreq_ghz * sol_kms
@@ -1072,7 +1070,7 @@ def suggest_extraction_scheme(
             # Figure out the binfactor
             this_binfactor = int(np.floor(target_chan_kms/chan_width_kms))
             # clamp to nchan if the binfactor exceeds the number of channels in the spw
-            nchan_spw = vm.spwInfo[this_spw]['numChannels']
+            nchan_spw = spw_info[this_spw]['numChannels']
             if this_binfactor > nchan_spw:
                 this_binfactor = nchan_spw
 
@@ -1536,7 +1534,7 @@ def build_mstransform_call(
     # Determine the column to use
 
     if datacolumn is None:
-        mytb = au.createCasaTool(casaStuff.tbtool)
+        mytb = casaStuff.tbtool()
         mytb.open(infile, nomodify = True)
         colnames = mytb.colnames()
         if 'CORRECTED_DATA' in colnames:
@@ -1596,7 +1594,7 @@ def build_mstransform_call(
             vlow_kms=vstart_kms, vhigh_kms=vstart_kms+vwidth_kms)
         # line_freq_ghz = (line_low_ghz+line_high_ghz)*0.5
 
-        max_chan_ghz = np.max(np.abs(au.getChanWidths(infile, spw)))/1e9
+        max_chan_ghz = np.max(np.abs(get_chan_widths(infile, spw)))/1e9
 
         # Using RADIO velocity convention:
         current_chan_kms = max_chan_ghz/restfreq_ghz*sol_kms
@@ -1740,7 +1738,7 @@ def reweight_data(
     # Determine column to use
 
     if datacolumn is None:
-        mytb = au.createCasaTool(casaStuff.tbtool)
+        mytb = casaStuff.tbtool()
         mytb.open(infile, nomodify = True)
         colnames = mytb.colnames()
         if 'CORRECTED_DATA' in colnames:
@@ -1757,15 +1755,15 @@ def reweight_data(
     exclude_str = ''
     if (edge_chans is not None) or (edge_kms is not None):
 
-        vm = au.ValueMapping(infile)
+        spw_info = get_spw_info(infile)
 
         first = True
-        for this_spw in vm.spwInfo.keys():
+        for this_spw in spw_info.keys():
 
             if edge_kms is not None:
-                spw_high_ghz = np.max(vm.spwInfo[this_spw]['edgeChannels'])/1e9
-                spw_low_ghz = np.min(vm.spwInfo[this_spw]['edgeChannels'])/1e9
-                spw_chanwidth_ghz = abs(vm.spwInfo[this_spw]['chanWidth'])/1e9
+                spw_high_ghz = np.max(spw_info[this_spw]['edgeChannels'])/1e9
+                spw_low_ghz = np.min(spw_info[this_spw]['edgeChannels'])/1e9
+                spw_chanwidth_ghz = abs(spw_info[this_spw]['chanWidth'])/1e9
 
                 mean_freq_ghz = 0.5*(spw_high_ghz+spw_low_ghz)
 
@@ -1787,7 +1785,7 @@ def reweight_data(
 
                 edge_chans = int(np.ceil(edge_kms / mean_chanwidth_kms))
 
-            nchan = vm.spwInfo[this_spw]['numChannels']
+            nchan = spw_info[this_spw]['numChannels']
 
             if edge_chans*2 > nchan:
                 logger.warning(
@@ -1999,7 +1997,7 @@ def batch_extract_continuum(
                     # all SPWs except the first one.
 
                     if clear_pointing:
-                        au.clearPointingTable(this_outfile)
+                        clear_pointing_table(this_outfile)
 
     # Concatenate and combine the output data sets
 
@@ -2107,7 +2105,7 @@ def extract_continuum(
         if not os.path.isdir(outfile + '.touch'):
             os.mkdir(outfile + '.touch')
 
-        mytb = au.createCasaTool(casaStuff.tbtool)
+        mytb = casaStuff.tbtool()
         mytb.open(infile + '.temp_copy', nomodify=True)
         colnames = mytb.colnames()
         if 'CORRECTED_DATA' in colnames:
@@ -2316,9 +2314,9 @@ def noise_spectrum(
         return None
 
     # Note the number of channels in SPW 0
+    spw_info = get_spw_info(vis)
 
-    vm = au.ValueMapping(vis)
-    nchan = vm.spwInfo[0]['numChannels']
+    nchan = spw_info[0]['numChannels']
     spec = np.zeros(nchan)
     for ii in range(nchan):
         if start_chan is not None:
@@ -2402,89 +2400,12 @@ def estimate_mrs(
 
         tb_subset.close()
 
-        tb_subset = tb.query(" && ".join(subset_conditions))
-        total_rows = tb_subset.nrows()
-        
-        uv_distance_m = []
-
-        # Loop over in chunks to reduce memory cost
-        for startrow in range(0, total_rows, chunk_size):
-            nrow = min(chunk_size, total_rows - startrow)
-
-            uvw = np.asarray(
-                tb_subset.getcol("UVW",
-                                   startrow=startrow,
-                                   nrow=nrow,
-                                   rowincr=1,
-                                   )
-            )
-            antenna1 = np.asarray(
-                tb_subset.getcol(
-                    "ANTENNA1",
-                    startrow=startrow,
-                    nrow=nrow,
-                    rowincr=1,
-                )
-            )
-            antenna2 = np.asarray(
-                tb_subset.getcol(
-                    "ANTENNA2",
-                    startrow=startrow,
-                    nrow=nrow,
-                    rowincr=1,
-                )
-            )
-            ddid = np.asarray(
-                tb_subset.getcol(
-                    "DATA_DESC_ID",
-                    startrow=startrow,
-                    nrow=nrow,
-                    rowincr=1,
-                )
-            )
-            flag_row = np.asarray(
-                tb_subset.getcol(
-                    "FLAG_ROW",
-                    startrow=startrow,
-                    nrow=nrow,
-                    rowincr=1,
-                ),
-                dtype=bool,
-            )
-            flags = np.asarray(
-                tb_subset.getcol(
-                    "FLAG",
-                    startrow=startrow,
-                    nrow=nrow,
-                    rowincr=1,
-                ),
-                dtype=bool,
-            )
-
-            # CASA normally returns UVW as (3, nrow).
-            if uvw.shape[0] != 3 and uvw.shape[-1] == 3:
-                uvw = uvw.T
-
-            uvdm = np.hypot(uvw[0], uvw[1])
-
-            # FLAG normally has dimensions (ncorr, nchan, nrow).
-            # Keep a row if at least one correlation/channel is unflagged.
-            flag_axes = tuple(range(flags.ndim - 1))
-            completely_flagged = np.all(flags, axis=flag_axes)
-
-            valid = (
-                ~flag_row
-                & ~completely_flagged
-                & (antenna1 != antenna2)
-                & np.isfinite(uvdm)
-                & (uvdm > 0)
-                & (ddid >= 0)
-            )
-
-            # Only take valid values
-            uv_distance_m.extend(uvdm[valid])
-
-        uv_distance_m = np.asarray(uv_distance_m)
+        # Get UV distances in meters
+        uv_distance_m = get_uv_distance_m(
+            tb=tb,
+            subset_conditions=subset_conditions,
+            chunk_size=chunk_size,
+        )
 
         # Now we loop over, take the 5th percentile baseline
         baseline_percentiles = []
@@ -2516,6 +2437,218 @@ def estimate_mrs(
 
     return result
 
+def estimate_synthesised_beam(
+    vis: str,
+    baseline_percentile: float = 80.0,
+    beam_factor: float = 0.574,
+    use_first_field: bool = True,
+    chunk_size: int = 100_000,
+) -> float:
+    """
+    Estimate the synthesised beam for a measurement set.
+
+    This function is specifically designed to sidestep the biases that can arise
+    from concatenating measurement sets. It calculates some minimum baseline percentile
+    from each unique observation ID, and then uses the minimum of these to for the beam size.
+    By default, this function uses 0.574 * 80th baseline percentile, to match the analysisUtils default.
+
+    Args:
+        vis (str): Path to measurement set
+        baseline_percentile (float): The percentile of the baseline distribution to use.
+            Defaults to 80.
+        beam_factor (float): Factor to convert from baseline percentile to synthesised beam.
+            Defaults to 0.574.
+        use_first_field (bool): Whether to only use the first field.
+            Defaults to True.
+        chunk_size (int): Number of rows to process in each chunk, to limit memory usage.
+            Defaults to 100,000.
+
+    Returns:
+        float: The estimated synthesised beam in arcseconds
+    """
+
+    # Get the representative frequency
+    rep_freq = get_representative_freq(vis)
+    logger.debug(f"Representative frequency: {rep_freq}")
+
+    # Get a list of the observation IDs
+    obs_ids = get_obs_ids(vis=vis)
+    logger.debug(f"Found {len(obs_ids)} observation IDs")
+
+    tb = casaStuff.tbtool()
+    tb.open(vis)
+
+    # Set up a dictionary to hold all the calculations
+    result = {
+        "representative_frequency": rep_freq.to(u.GHz).value,
+        "synthesised_beam_per_obs_id": {},
+    }
+
+    for obs_id in obs_ids:
+
+        # Downselect on observation ID
+        subset_conditions = [
+            f"OBSERVATION_ID == {obs_id}",
+        ]
+        tb_subset = tb.query(" && ".join(subset_conditions))
+
+        # If we're only using first field, get the first field ID for further downselecting
+        if use_first_field:
+            field_id = tb_subset.getcell("FIELD_ID", 0)
+            subset_conditions.append(f"FIELD_ID == {field_id}")
+
+        tb_subset.close()
+
+        # Get UV distances in meters
+        uv_distance_m = get_uv_distance_m(
+            tb=tb,
+            subset_conditions=subset_conditions,
+            chunk_size=chunk_size,
+        )
+        baseline_for_beam = np.nanpercentile(uv_distance_m, baseline_percentile) * u.m
+
+        # Calculate the synthesised beam in arcsec
+        synthesised_beam = (
+            beam_factor
+            * const.c.to(u.m * u.Hz)
+            / (rep_freq.to(u.Hz) * baseline_for_beam.to(u.m))
+            * u.rad
+        )
+        synthesised_beam = synthesised_beam.to(u.arcsec).value
+
+        result["synthesised_beam_per_obs_id"][obs_id] = synthesised_beam
+
+        logger.debug(f"Calculated synthesised beam of {synthesised_beam} for obs ID {obs_id}")
+
+    tb.close()
+
+    # The synthesised beam is the minimum of the synthesised beams from each observation ID
+    synthesised_beam = np.nanmin([r[-1] for r in result["synthesised_beam_per_obs_id"].items()])
+
+    return synthesised_beam
+
+
+def pick_cell_and_im_size(
+    vis: str,
+    npix: float = 5.0,
+    baseline_percentile: float = 80.0,
+    cellstring: bool = False,
+    roundcell: int = 2,
+    pblevel: float = 0.2,
+    beam_size_factor: float = 0.574,
+    beam_size_use_first_field: bool = True,
+    beam_size_chunk_size: int = 100_000,
+) -> tuple[float | str, list[int]]:
+    """Pick a cell size and image size for a measurement set.
+
+    Args:
+        vis (str): Path to measurement set.
+        npix (float): Number of pixels across the synthesised beam.
+            Defaults to 5.0.
+        baseline_percentile (float): The percentile of the baseline distribution to use for estimating the synthesised
+            beam. Defaults to 80.0.
+        cellstring (bool): Whether to return the cell size as a string with units.
+            Defaults to False.
+        roundcell (int): Number of decimal places to round the cell size to.
+            Defaults to 2.
+        pblevel (float): The primary beam level for calculating image size.
+            Defaults to 0.2.
+        beam_size_factor (float): Factor to convert from baseline percentile to synthesised beam.
+            Defaults to 0.574.
+        beam_size_use_first_field (bool): Whether to only use the first field for estimating the synthesised beam.
+            Defaults to True.
+        beam_size_chunk_size (int): Number of rows to process in each chunk when estimating the synthesised beam.
+            Defaults to 100,000.
+
+    Returns:
+        float or str, list: Tuple of the cell size (potentially as a string) and an associated list of imsize.
+    """
+
+    # Get cell size from synthesised beam and oversample factor
+    synthesised_beam = estimate_synthesised_beam(
+        vis=vis,
+        baseline_percentile=baseline_percentile,
+        beam_factor=beam_size_factor,
+        use_first_field=beam_size_use_first_field,
+        chunk_size=beam_size_chunk_size,
+    )
+    cell_size = synthesised_beam / npix
+
+    # If selected, round to requested number of significant figures
+    if roundcell > 0:
+        cell_size = float(f"{cell_size:.{roundcell}g}")
+
+    # Convert to a cell size with units
+    cell_size_unit = cell_size * u.arcsec
+        
+    if cellstring:
+        cell_size = f"{cell_size}arcsec"
+
+    # Get dish diameter for image size calculation. Use max to get minimum FOV
+    tb = casaStuff.table()
+    tb.open(os.path.join(f"{vis}/ANTENNA"))
+    dish_dia = tb.getcol("DISH_DIAMETER").max()
+    dish_dia = dish_dia * u.m
+    tb.close()
+
+    # Get representative frequency
+    rep_freq = get_representative_freq(vis)
+
+    # Get the source name from the MS metadata
+    msmd = casaStuff.msmdtool()
+    msmd.open(vis)
+
+    # Get a list of all science field IDs, turn into a source name by taking the first
+    # in the MS
+    intents = msmd.intents()
+    field_ids = msmd.fieldsforintent(intents[0])
+    sourcename = msmd.namesforfields(field_ids)[0]
+    logger.debug(f"Using source name {sourcename} for image size calculation")
+    msmd.close()
+    
+    tb = casaStuff.table()
+    tb.open(f"{vis}/FIELD")
+    tb_subset = tb.query(f"NAME == '{sourcename}'")
+    phase_dirs = tb_subset.getcol("PHASE_DIR")
+
+    tb_subset.close()
+    tb.close()
+
+    # Extract RA and Dec arrays (in radians)
+    ra_fields = phase_dirs[0, 0, :] * u.rad
+    dec_fields = phase_dirs[1, 0, :] * u.rad
+    
+    # Get primary beam FWHM
+    fwhm = (1.14 * 1.22 * const.c.to(u.m * u.Hz) / rep_freq.to(u.Hz) / dish_dia) * u.rad
+    
+    # Apply the PB cutoff
+    radius_pb = (fwhm / 2.0) * np.sqrt(np.log(pblevel) / np.log(0.5))
+
+    ra_mins = ra_fields - radius_pb
+    ra_maxs = ra_fields + radius_pb
+    dec_mins = dec_fields - radius_pb
+    dec_maxs = dec_fields + radius_pb
+
+    ra_extent = ra_maxs.max() - ra_mins.min()
+    dec_extent = dec_maxs.max() - dec_mins.min()
+
+    x_extent = ra_extent / cell_size_unit
+    y_extent = dec_extent / cell_size_unit
+
+    # Round this up to a good FFT number for the FFT,
+    # cast to int
+    x_extent = int(1.2 * x_extent)
+    y_extent = int(1.2 * y_extent)
+
+    su = casaStuff.synthesisutils()
+    x_extent = int(su.getOptimumSize(x_extent))
+    y_extent = int(su.getOptimumSize(y_extent))
+    su.done()
+
+    im_size = [x_extent, y_extent]
+
+    return cell_size, im_size
+
 def get_representative_freq(
         vis: str,
 ) -> u.Quantity:
@@ -2592,3 +2725,439 @@ def get_obs_ids(vis):
     obs_ids = list(range(nrows))
 
     return obs_ids
+
+
+def get_chan_widths(
+    vis: str,
+    spw: str | int | list = "",
+    velocity: bool = False,
+) -> np.ndarray:
+    """Get channel widths for a given measurement set and SPW selection.
+
+    Args:
+        vis (str): Path to measurement set.
+        spw (str | int | list, optional): SPW selection. Can be a string, integer, or list of integers.
+            Defaults to "" (all SPWs).
+        velocity (bool, optional): If True, return channel widths in km/s.
+            Defaults to False.
+
+    Returns:
+        np.ndarray: Array of channel widths for the selected SPWs. If velocity is True,
+            the widths are in velocity units; otherwise, they are in Hz.
+    """
+
+    if not os.path.exists(vis):
+        raise FileNotFoundError(f"Measurement set '{vis}' not found.")
+
+    msmd = casaStuff.msmdtool()
+    msmd.open(vis)
+
+    # If an empty string or list is passed, selected all SPWs
+    if spw in ["", []]:
+        spw = list(range(msmd.nspw()))
+
+    # If we still have a string, turn to int
+    if isinstance(spw, str):
+        spw = [int(i) for i in spw.split(",")]
+
+    # Make sure spw is a list
+    if not isinstance(spw, list):
+        spw = [spw]
+
+    chan_widths = []
+    mean_freqs = []
+    for s in spw:
+        chan_width = np.median(msmd.chanwidths(s))
+        mean_freq = msmd.meanfreq(s)
+
+        chan_widths.append(chan_width)
+        mean_freqs.append(mean_freq)
+
+    chan_widths = chan_widths * u.Hz
+
+    # Convert to velocity units if requested
+    if velocity:
+        mean_freqs = mean_freqs * u.Hz
+        chan_widths *= const.c.to(u.m * u.Hz) / mean_freqs
+        chan_widths = chan_widths.to(u.km / u.s).value
+    else:
+        chan_widths = chan_widths.to(u.Hz).value
+
+    msmd.close()
+
+    return chan_widths
+
+def get_uv_distance_m(
+    tb: casaStuff.tbtool,
+    subset_conditions: list,
+    chunk_size: int = 100_000,
+) -> np.ndarray:
+    """Get the UV distance in meters for a given table and subset conditions.
+
+    Args:
+        tb (casaStuff.tbtool): The CASA table tool object.
+        subset_conditions (list): List of conditions to subset the table.
+        chunk_size (int): Number of rows to process in each chunk, to limit memory usage.
+        
+    Returns:
+        np.ndarray: Array of UV distances in meters.
+    """
+
+    tb_subset = tb.query(" && ".join(subset_conditions))
+    total_rows = tb_subset.nrows()
+
+    uv_distance_m = []
+
+    # Loop over in chunks to reduce memory cost
+    for startrow in range(0, total_rows, chunk_size):
+        nrow = min(chunk_size, total_rows - startrow)
+
+        uvw = np.asarray(
+            tb_subset.getcol(
+                "UVW",
+                startrow=startrow,
+                nrow=nrow,
+                rowincr=1,
+            )
+        )
+        antenna1 = np.asarray(
+            tb_subset.getcol(
+                "ANTENNA1",
+                startrow=startrow,
+                nrow=nrow,
+                rowincr=1,
+            )
+        )
+        antenna2 = np.asarray(
+            tb_subset.getcol(
+                "ANTENNA2",
+                startrow=startrow,
+                nrow=nrow,
+                rowincr=1,
+            )
+        )
+        ddid = np.asarray(
+            tb_subset.getcol(
+                "DATA_DESC_ID",
+                startrow=startrow,
+                nrow=nrow,
+                rowincr=1,
+            )
+        )
+        flag_row = np.asarray(
+            tb_subset.getcol(
+                "FLAG_ROW",
+                startrow=startrow,
+                nrow=nrow,
+                rowincr=1,
+            ),
+            dtype=bool,
+        )
+
+        # Sometimes trying to pull out the FLAG column fails.
+        # In that case, just make a dummy array of all False flags
+        try:
+            flags = np.asarray(
+                tb_subset.getcol(
+                    "FLAG",
+                    startrow=startrow,
+                    nrow=nrow,
+                    rowincr=1,
+                ),
+                dtype=bool,
+            )
+        except RuntimeError:
+            flags = np.zeros((1, 1, nrow), dtype=bool)
+
+        # CASA normally returns UVW as (3, nrow).
+        if uvw.shape[0] != 3 and uvw.shape[-1] == 3:
+            uvw = uvw.T
+
+        uvdm = np.hypot(uvw[0], uvw[1])
+
+        # FLAG normally has dimensions (ncorr, nchan, nrow).
+        # Keep a row if at least one correlation/channel is unflagged.
+        flag_axes = tuple(range(flags.ndim - 1))
+        completely_flagged = np.all(flags, axis=flag_axes)
+
+        valid = (
+            ~flag_row
+            & ~completely_flagged
+            & (antenna1 != antenna2)
+            & np.isfinite(uvdm)
+            & (uvdm > 0)
+            & (ddid >= 0)
+        )
+
+        # Only take valid values
+        uv_distance_m.extend(uvdm[valid])
+
+    tb_subset.close()
+
+    uv_distance_m = np.asarray(uv_distance_m)
+
+    return uv_distance_m
+
+def get_spw_info(
+    vis: str,
+    ignore_wvr: bool = True,
+) -> dict:
+    """Get SPW information from a measurement set.
+
+    Args:
+        vis (str): Path to measurement set.
+        ignore_wvr (bool): Whether to ignore WVR SPWs. Defaults to True.
+
+    Returns:
+        dict: Dictionary containing SPW information, including bandwidth, channel frequencies,
+            channel width, edge channels, sideband, mean frequency, and number of channels.
+    """
+
+    spw_info = {}
+
+    mytb = casaStuff.tbtool()
+    mytb.open(f"{vis}/SPECTRAL_WINDOW")
+
+    # Keep track of number of rows in the table
+    nrows = range(mytb.nrows())
+
+    for i in nrows:
+
+        # Get out useful info
+        bandwidth = mytb.getcell("TOTAL_BANDWIDTH", i)
+        chan_freqs = mytb.getcell("CHAN_FREQ", i)
+        min_freq = min(chan_freqs)
+        max_freq = max(chan_freqs)
+        mean_freq = chan_freqs.mean()
+        num_channels = chan_freqs.shape[0]
+        chan_width = mytb.getcell("CHAN_WIDTH", i)[0]
+        net_sideband = mytb.getcell("NET_SIDEBAND", i)
+
+        # Put this all into a dictionary
+        spw_info[i] = {}
+        spw_info[i]["bandwidth"] = bandwidth
+        spw_info[i]["chanFreqs"] = chan_freqs
+        spw_info[i]["chanWidth"] = chan_width
+        spw_info[i]["edgeChannels"] = [min_freq, max_freq]
+        if net_sideband == 2:
+            spw_info[i]["sideband"] = 1
+        else:
+            spw_info[i]["sideband"] = -1
+        spw_info[i]["meanFreq"] = mean_freq
+        spw_info[i]["numChannels"] = num_channels
+
+        # If ignoring WVR, then remove and log this
+        if ignore_wvr and (num_channels == 4):
+            logger.debug(f"Ignoring spectral window {i} because it is WVR related")
+            spw_info.pop(i)
+
+    mytb.close()
+
+    return spw_info
+
+def get_scans_for_spw(
+        vis: str,
+) -> dict:
+    """Get the unique scan numbers for each spectral window in a measurement set.
+
+    Args:
+        vis (str): Path to measurement set.
+
+    Returns:
+        dict: Dictionary mapping SPW IDs to unique scan numbers.
+    """
+
+    # Get out SPWs for each data description ID
+    tb = casaStuff.tbtool()
+    tb.open(f"{vis}/DATA_DESCRIPTION")
+    spw_for_data_desc_id = tb.getcol("SPECTRAL_WINDOW_ID")
+    tb.close()
+
+    tb = casaStuff.tbtool()
+    tb.open(vis)
+
+    data_desc_id = tb.getcol('DATA_DESC_ID')
+    scans = tb.getcol('SCAN_NUMBER')
+
+    scans_for_spw = {}
+    for i in spw_for_data_desc_id:
+        spw = spw_for_data_desc_id[i]
+        indices = np.where(data_desc_id == i)
+        scans_for_spw[spw] = np.unique(scans[indices])
+        
+    tb.close()
+
+    return scans_for_spw
+
+def get_science_spws(
+    vis: str,
+    intent: str = "OBSERVE_TARGET#ON_SOURCE",
+    return_string: bool = True,
+    return_list_of_strings: bool = False,
+    return_freq_ranges: bool = False,
+    tdm: bool = True,
+    fdm: bool = True,
+    sqld: bool = False,
+    chavg: bool = False,
+) -> str | list[int | str] | dict:
+    """ Get science SPWs from a measurement set based.
+
+    Will select SPWs based on specified intent. For ALMA data,
+    it also can ignore channel-averaged and SQLD SPWs.
+
+    Args:
+        vis (str): Path to measurement set.
+        intent (str): The intent to filter by.
+            Defaults to "OBSERVE_TARGET#ON_SOURCE".
+        return_string (bool): If True, return a comma-separated string of SPWs.
+            Defaults to True.
+        return_list_of_strings (bool): If True, return a list of SPWs as strings.
+            Defaults to False.
+        return_freq_ranges (bool): If True, return a dictionary with frequency ranges.
+            Defaults to False.
+        tdm (bool): If True, include TDM SPWs.
+            Defaults to True.
+        fdm (bool): If True, include FDM SPWs.
+            Defaults to True.
+        sqld (bool): If True, include SQLD SPWs.
+            Defaults to False.
+        chavg (bool): If True, include channel-averaged SPWs.
+            Defaults to False.
+
+    Returns:
+        str: Comma-separated string of science SPWs if return_string is True.
+        list[int]: List of science SPWs as integers if return_list_of_strings is False.
+        list[str]: List of science SPWs as strings if return_list_of_strings is True.
+        dict: Dictionary of science SPWs with frequency ranges if return_freq_ranges is True.
+    """
+
+    if return_string and return_list_of_strings:
+        raise ValueError("You can only specify one of: return_string, return_list_of_strings")
+
+    msmd = casaStuff.msmdtool()
+    msmd.open(vis)
+
+    all_intents = msmd.intents()
+
+    if intent not in all_intents and intent != "":
+        for i in all_intents:
+            if i.find(intent) >= 0:
+                intent = i
+                logger.debug(f"Translated intent to {i}")
+                break
+
+    # Minimum match OBSERVE_TARGET to OBSERVE_TARGET#UNSPECIFIED
+    value = [i.find(intent.replace("*", "")) for i in all_intents]
+
+    # If any intent gives a match, the mean value of the location list will be > -1
+    if np.mean(value) == -1 and intent != "":
+        logger.warning(f"{intent} not found in this dataset. Available intents: {all_intents}")
+
+        if return_string:
+            science_spws = ""
+        else:
+            science_spws = []
+
+    else:
+        # If we don't have an intent, match on wildcards
+        if intent == "":
+            intent = "*"
+
+        # Get SPWs for intent, cast back to int
+        spws = msmd.spwsforintent(intent)
+        spws = [int(i) for i in spws]
+
+        # Get observatory name, if ALMA data then perform further filtering
+        observatory_name = get_observatory_name(vis)
+
+        if observatory_name.find("ALMA") >= 0 or observatory_name.find("OSF") >= 0:
+
+            logger.debug(
+                f"Calling almaspws with:\n  chavg={chavg}\n  tdm={tdm}\n  fdm={fdm}\n  sqld={sqld}"
+            )
+            alma_spws = msmd.almaspws(
+                chavg=chavg,
+                tdm=tdm,
+                fdm=fdm,
+                sqld=sqld,
+            )
+            alma_spws = [int(i) for i in alma_spws]
+            
+            if chavg and not sqld:
+                sqld_spws = msmd.almaspws(sqld=True)
+                sqld_spws = [int(i) for i in sqld_spws]
+
+                logger.debug(f"removing SQLD ({sqld_spws}) from list")
+                alma_spws = list(set(alma_spws) - set(sqld_spws))
+            if len(spws) == 0 or len(alma_spws) == 0:
+                science_spws = []
+            else:
+                science_spws = [int(i) for i in np.intersect1d(spws, alma_spws)]
+        else:
+            science_spws = spws
+
+        science_spws_dict = {}
+        for spw in science_spws:
+            science_spws_dict[spw] = sorted([msmd.chanfreqs(spw)[0], msmd.chanfreqs(spw)[-1]])
+
+        msmd.close()
+
+        if return_freq_ranges:
+            science_spws = science_spws_dict
+        if return_string:
+            science_spws = ",".join(str(i) for i in science_spws)
+        elif return_list_of_strings:
+            science_spws = list([str(i) for i in science_spws])
+        else:
+            science_spws = list(science_spws)
+
+    return science_spws
+
+def get_observatory_name(
+        vis: str,
+) -> str:
+    """Get observatory name for given measurement set.
+
+    Args:
+        vis (str): The path to the measurement set.
+
+    Returns:
+        str: The observatory name, if found; otherwise, an empty string.
+    """
+
+    tb = casaStuff.tbtool()
+    tb.open(f"{vis}/OBSERVATION")
+
+    try:
+        observatory_name = tb.getcell("TELESCOPE_NAME")
+    except RuntimeError:
+        observatory_name = ""
+
+    tb.close()
+
+    return observatory_name
+
+def clear_pointing_table(
+    vis: str,
+):
+    """Removes all rows from the POINTING table of a measurement set.
+
+    Args:
+        vis (str): Path to the measurement set.
+
+    Returns:
+        bool: True if the POINTING table was cleared successfully.
+    """
+
+    tb = casaStuff.tbtool()
+    tb.open(
+        f"{vis}/POINTING",
+        nomodify=False,
+    )
+    row_numbers = tb.rownumbers()
+    tb.removerows(row_numbers)
+    tb.close()
+
+    logger.info(f"Cleared the POINTING table for {vis}")
+
+    return True
diff --git a/phangsPipeline/clean_call.py b/phangsPipeline/clean_call.py
index 12870db..23b66c0 100644
--- a/phangsPipeline/clean_call.py
+++ b/phangsPipeline/clean_call.py
@@ -6,11 +6,10 @@
 import os
 import re
 
-import analysisUtils as au
 import numpy as np
 
 from .casaStuff import imhead
-from .casaVisRoutines import estimate_mrs
+from .casaVisRoutines import estimate_mrs, estimate_synthesised_beam
 
 logger = logging.getLogger(__name__)
 
@@ -276,7 +275,7 @@ def estimate_synthesised_beam(
             logger.info("Estimating beam from visibilities")
             
             vis = self.get_param("vis")
-            beam = au.estimateSynthesizedBeam(vis)
+            beam = estimate_synthesised_beam(vis)
 
         return beam
 
diff --git a/phangsPipeline/handlerImaging.py b/phangsPipeline/handlerImaging.py
index e0084b6..7ba4ebf 100644
--- a/phangsPipeline/handlerImaging.py
+++ b/phangsPipeline/handlerImaging.py
@@ -57,9 +57,6 @@
 
 if casa_enabled:
 
-    # Analysis utilities
-    import analysisUtils as au
-
     from .clean_call import CleanCall, CleanCallFunctionDecorator
 
     from . import casaImagingRoutines as imr
diff --git a/phangsPipeline/handlerImagingChunked.py b/phangsPipeline/handlerImagingChunked.py
index 2ab62ee..183c9de 100644
--- a/phangsPipeline/handlerImagingChunked.py
+++ b/phangsPipeline/handlerImagingChunked.py
@@ -29,8 +29,6 @@
 
 if casa_enabled:
 
-    import analysisUtils as au
-
     from . import casaImagingRoutines as imr
     from . import casaMaskingRoutines as msr
     from . import casaStuff
@@ -539,7 +537,7 @@ def return_valid_chunks(self, chunk_num=None):
         #     chunks_iter = self.return_valid_chunks(chunk_num=chunk_num)
 
         #     # TODO: move to helper function.
-        #     mytb = au.createCasaTool(casaStuff.tbtool)
+        #     mytb = casaStuff.tbtool()
         #     mytb.open(self.full_vis_file, nomodify = True)
         #     colnames = mytb.colnames()
         #     if 'CORRECTED_DATA' in colnames:
@@ -590,8 +588,8 @@ def task_split_cube_to_chunks(self, imagename, overwrite=False, chunk_num=None,
 
             chunks_iter = self.return_valid_chunks(chunk_num=chunk_num)
 
-            myia = au.createCasaTool(casaStuff.iatool)
-            myrg = au.createCasaTool(casaStuff.rgtool)
+            myia = casaStuff.iatool()
+            myrg = casaStuff.rgtool()
 
             myia.open(imagename)
 
diff --git a/phangsPipeline/handlerTestImaging.py b/phangsPipeline/handlerTestImaging.py
index 4d6ef7b..0d7c7b6 100644
--- a/phangsPipeline/handlerTestImaging.py
+++ b/phangsPipeline/handlerTestImaging.py
@@ -55,9 +55,6 @@
 
 if casa_enabled:
 
-    # Analysis utilities
-    import analysisUtils as au
-
     from .clean_call import CleanCall
 
     from . import utilsLines as lines
@@ -257,7 +254,7 @@ def _get_beam_properties(
                 logger.warning(f"PSF file not found: {psf_file}")
                 return None
 
-            myia = au.createCasaTool(casaStuff.iatool)
+            myia = casaStuff.iatool()
             myia.open(psf_file)
 
             beam_info = {}
diff --git a/phangsPipeline/utilsFieldSelection.py b/phangsPipeline/utilsFieldSelection.py
index e4b3ee1..7ac6081 100644
--- a/phangsPipeline/utilsFieldSelection.py
+++ b/phangsPipeline/utilsFieldSelection.py
@@ -5,12 +5,12 @@
 import re
 import shutil
 
-import analysisUtils as aU
 import numpy as np
 
 from . import casaStuff
+from .casaVisRoutines import get_science_spws
 
-tb = aU.createCasaTool(casaStuff.tbtool)
+tb = casaStuff.tbtool()
 split = casaStuff.split
 
 # 
@@ -108,8 +108,14 @@ def extract_field_selections(
     # 
     # compute pb
     if pb_fwhm is None:
-        spw_dict = aU.getScienceSpws(vis, intent='OBSERVE_TARGET#ON_SOURCE', returnFreqRanges=True)
-        #print('spw_dict.values()', spw_dict.values())
+        spw_dict = get_science_spws(
+            vis,
+            intent='OBSERVE_TARGET#ON_SOURCE',
+            return_freq_ranges=True,
+        )
+        if not isinstance(spw_dict, dict):
+            raise TypeError('Expecting spw_dict to be a dictionary')
+
         min_freq = np.min(np.array(list(spw_dict.values())).ravel())
         pb_fwhm = 1.14*1.22*(3e8/min_freq)*3600*180/(ant_diam*3.1415926)
     match_radius = pb_factor * pb_fwhm / 2.0 # arcsec
diff --git a/phangsPipeline/utilsSingleDish.py b/phangsPipeline/utilsSingleDish.py
index b4be142..5aea694 100644
--- a/phangsPipeline/utilsSingleDish.py
+++ b/phangsPipeline/utilsSingleDish.py
@@ -1,5 +1,7 @@
 import copy
+import logging
 import os
+import re
 
 import analysisUtils as au
 import astropy.units as u
@@ -8,19 +10,17 @@
 import pylab as pb
 from analysisUtils import mjdSecondsListToDateTime, mjdsecToUT
 
+from . import casaStuff
 from .casaStuff import (
     casa_version,
-    tbtool,
-    msmdtool,
-    metool,
-    qatool,
-    iatool,
     imhead,
     importfits,
     imregrid,
     imtrans,
 )
 
+logger = logging.getLogger(__name__)
+
 casaVersion = "{0}.{1}.{2}".format(*casa_version)
 
 def get_first_arr_val(arr):
@@ -69,7 +69,7 @@ def prepare_sd_image(
     os.system('mv -f ' + sd_image_file + '_reorder ' + sd_image_file)
 
     # Regrid this to the input measurement set to avoid any weirdness with overlap.
-    mytb = au.createCasaTool(tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(clean_call.get_param('vis') + '/SPECTRAL_WINDOW')
     freq = mytb.getcol('CHAN_FREQ')
     n_chan = mytb.getcol('NUM_CHAN')[0]
@@ -108,7 +108,7 @@ def prepare_sd_image(
     n_chan = cube_info['shape'][-1]
     n_pol = cube_info['shape'][-2]
 
-    myia = au.createCasaTool(iatool)
+    myia = casaStuff.iatool()
     myia.open(sd_image_file)
     restoring_beam = myia.restoringbeam()
     myia.setrestoringbeam(remove=True)
@@ -137,7 +137,7 @@ def get_dish_diameter(
     freq = hdr["crval4"] * u.Unit(hdr["cunit4"])
 
     # Get restoring beam in radians
-    myia = au.createCasaTool(iatool)
+    myia = casaStuff.iatool()
     myia.open(infile=sdimage)
 
     # Beam size is the same for all channels, so just take the first
@@ -308,7 +308,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False,
     if field == '':
         field = None
 
-    mytb = au.createCasaTool(tbtool)
+    mytb = casaStuff.tbtool()
     mytb.open(vis)
     allObsIDs = mytb.getcol('OBSERVATION_ID')
     obsIDs = np.unique(allObsIDs)
@@ -329,7 +329,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False,
         print("This function requires casa >= 4.1.0.")
         return
 
-    mymsmd = au.createCasaTool(msmdtool)
+    mymsmd = casaStuff.msmdtool()
     mymsmd.open(vis)
     timeranges = {}
     timecenters = []
@@ -653,7 +653,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False,
     if (coordSys.find('AZEL') >= 0 and convert):
         coordSys = 'J2000'
         print("Converting %d coordinates from AZELGEO to J2000..." % (len(x)))
-        my_metool = au.createCasaTool(metool)
+        my_metool = casaStuff.metool()
         for i in range(len(x)):
             if ((i+1) % 10000 == 0): print("%d/%d" % (i+1,len(x)))
             ra, dec = au.computeRADecFromAzElMJD([xrad[i],yrad[i]], mjd=times[i]/86400.,
@@ -1175,7 +1175,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False,
             sample = i*labelIncrement
             pb.text(x[sample],y[sample],str(sample),size=8)
 
-        myqa = au.createCasaTool(qatool)
+        myqa = casaStuff.qatool()
         if (coordSys.find('AZEL') >= 0):
             azimString = myqa.formxxx('%frad'%(rightAscension), format='deg', prec=5)
             elevString = myqa.formxxx('%frad'%(declination), format='deg', prec=5)
@@ -1252,3 +1252,70 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False,
     else:
         return xSampling, ySampling, largestDimension
     # end of getTPSampling
+
+# Get source name
+def get_sourcename(filename,
+                   source='all',
+                   ):
+
+    mytb   = casaStuff.msmdtool()
+    mytb.open(filename)
+
+    # Get source field names
+    fieldnames = mytb.fieldnames()
+
+    # If we're looking for all intents, then just take the (first) science target
+    if source == 'all':
+        field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')[0]
+
+    # Else, key in on the source name. Do this in a regex-style
+    else:
+
+        # For field selection, this can be done with forward slashes,
+        # so make sure to strip them
+        if source.startswith("/"):
+            source = source.strip("/")
+
+        # Get science fieldnames as names
+        fieldnames_asname = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE', asnames=True)
+
+        regex_matches = [re.fullmatch(source, x) for x in fieldnames_asname]
+        regex_idxs = [x is not None for x in regex_matches]
+        found_field = any(regex_idxs)
+
+        # If we've found a match, update the source name with the first
+        # fieldname
+        if found_field:
+            source = fieldnames_asname[regex_idxs][0]
+
+        if not found_field:
+            logger.warning(f"source {source} not in field names: {', '.join(fieldnames)}")
+            return None
+
+        field = mytb.fieldsforname(source)
+
+        # If we don't find anything, get the science targets
+        if len(field) == 0:
+            field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')
+
+        # Take the first result
+        field = field[0]
+
+    source_name = fieldnames[field]
+
+    mytb.close()
+
+    return source_name
+
+# Read source coordinates
+def read_source_coordinates(filename,source):
+
+    coord_source = au.getRADecForSource(filename,source)
+    RA_h  = (coord_source.split(' ')[0]).split(':')[0]
+    RA_m  = (coord_source.split(' ')[0]).split(':')[1]
+    RA_s  = (coord_source.split(' ')[0]).split(':')[2]
+    DEC_d = (coord_source.split(' ')[1]).split(':')[0]
+    DEC_m = (coord_source.split(' ')[1]).split(':')[1]
+    DEC_s = (coord_source.split(' ')[1]).split(':')[2]
+    coord = "J2000  "+str(RA_h)+"h"+str(RA_m)+"m"+str(RA_s[0:6])+" "+str(DEC_d)+"d"+str(DEC_m)+"m"+str(DEC_s)
+    return coord
diff --git a/phangsPipeline/utilsTestImagingPlots.py b/phangsPipeline/utilsTestImagingPlots.py
index ba28b4f..36b694c 100644
--- a/phangsPipeline/utilsTestImagingPlots.py
+++ b/phangsPipeline/utilsTestImagingPlots.py
@@ -61,7 +61,6 @@
 casa_enabled = is_casa_installed()
 
 if casa_enabled:
-    import analysisUtils as au
     from . import casaStuff
 
 from . import utilsFilenames
@@ -107,7 +106,7 @@ def _read_image_data(image_file):
     # Try CASA image
     if casa_enabled and os.path.isdir(image_file):
         try:
-            myia = au.createCasaTool(casaStuff.iatool)
+            myia = casaStuff.iatool()
             myia.open(image_file)
             data = myia.getchunk()
             csys = myia.coordsys()
diff --git a/tests/test_casa_imports.py b/tests/test_casa_imports.py
index 3e2c336..8dd1caf 100644
--- a/tests/test_casa_imports.py
+++ b/tests/test_casa_imports.py
@@ -390,27 +390,20 @@ class TestAnalysisUtilsImports:
     au_funcs = [
         "Atmcal",
         "MAD",
-        "ValueMapping",
         "angularSeparationRadians",
-        "clearPointingTable",
         "computeAzElFromRADecMJD",
         "computeRADecFromAzElMJD",
-        "createCasaTool",
         "dateStringToMJD",
         "getAntennaPadXYZ",
-        "getBaselineStats",
-        "getChanWidths",
         "getObservationStart",
         "getObservationStartDate",
         "getOffSourceTimes",
         "getRADecForField",
         "getRADecForSource",
-        "getScienceSpws",
         "lsrkToTopo",
         "mjdSecondsListToDateTime",
         "mjdSecondsToMJDandUT",
         "parseTimerangeArgument",
-        "pickCellSize",
         "radec2rad",
         "stuffForScienceDataReduction",
         "timeOnSource",
@@ -473,22 +466,6 @@ class TestAnalysisUtilsArguments:
                 "usemstool",
             ],
         ),
-        (
-            "getScienceSpws",
-            [
-                "intent",
-                "returnFreqRanges",
-            ],
-        ),
-        (
-            "pickCellSize",
-            [
-                "imsize",
-                "intent",
-                "npix",
-                "pblevel",
-            ],
-        ),
     ]
 
     @pytest.mark.casa