diff --git a/lofar_station_client/RCUs.py b/lofar_station_client/RCUs.py
new file mode 100644
index 0000000000000000000000000000000000000000..2637079d1228c6da62fe47f36d3904f47344e179
--- /dev/null
+++ b/lofar_station_client/RCUs.py
@@ -0,0 +1,6 @@
+RCU2L = [0, 1, 2, 3, 4, 5]
+RCU2H = [8, 9]
+# RCU2H = [8]
+RCU2Hpwr = [8]  # Should be even numbers
+RCU2Hctl = [9]  # Should be even numbers
+# RCU2Hctl = [rcu + 1 for rcu in RCU2Hpwr]  # This is then odd numbers
diff --git a/lofar_station_client/dts_outside.py b/lofar_station_client/dts_outside.py
new file mode 100644
index 0000000000000000000000000000000000000000..1e1df5ec55c56821eb94fb81866c0f8b42fa2c84
--- /dev/null
+++ b/lofar_station_client/dts_outside.py
@@ -0,0 +1,1161 @@
+# dts_outside.py: Module with functions to restart dts outside and
+# put it in a default configuration
+#
+# -*- coding: utf-8 -*-
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may
+# not use this file except in compliance with the License. You may obtain
+# a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations
+# under the License.
+
+
+"""  module dts_outside
+
+This module contains functions to:
+- restart the Dwingeloo Test Station (Outside)
+- configure the Dwingeloo Test Station (Outside) to a defined default
+
+Boudewijn Hut
+"""
+
+import time
+import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+import processing_lib as plib
+import RCUs
+import os
+
+
+# define some constants for the setup:
+n_fgpa = 4  # number of fpgas
+
+aps_location_labels = ['Slot%02d' % nr for nr in range(32)]
+# positions of RCU boards. (0,0) is bottom, left
+aps_rcu_pos_x = [11, 10, 9, 8, 7, 6] +\
+                [11, 10, 9, 8, 7]*2 +\
+                [4, 3, 2, 1, 0]*2 +\
+                [5, 4, 3, 2, 1, 0]
+aps_rcu_pos_y = [2]*6 + [1]*5 + [0]*5 + [2]*5 + [1]*5 + [0]*6
+# positions of the other modules in the APS
+aps_mod_pos_x = [0, 1, 3, 4, 5]
+aps_mod_pos_y = [0]*5
+
+# if the names are correctly provided by the RCU2s,
+# then update plot_subband_statistics() (see annotation there)
+# and remove these two lines:
+RCU2L_TAGS = ['8393812', '8416938', '8469028', '8386883', '8374523']
+RCU2H_TAGS = ['8514859', '8507272']
+
+plot_dir = 'inspection_plots'
+station_name = 'DTS-Outside'
+
+folders = [plot_dir]
+for folder in folders:
+    if not os.path.exists(folder):
+        os.makedirs(folder)
+
+
+def set_to_default_config(devices, skip_boot=False, silent=False):
+    '''
+    Restart system and set to default configuration
+
+    Input arguments:
+    devices = struct of software devices of m&c of the station
+              name: object
+
+    Output arguments:
+    foundError = True if an error ocurred, False otherwise
+    '''
+    foundError = False
+    # initialisation
+    if not skip_boot:
+        if not silent:
+            plib.log('Start initialisation')
+        foundError = foundError | initialise(devices)
+    else:
+        plib.log('SKIPPING INITIALISATION (skip_init was set to True)')
+    if not silent:
+        plib.log('Check initialisation')
+    foundError = foundError | check_initialise(devices)
+    # configuring
+    if not silent:
+        plib.log('Start setting to default configuration')
+    if not silent:
+        plib.log('Start configuring Low Receivers')
+    foundError = foundError | set_rcul_to_default_config(devices)
+    if not silent:
+        plib.log('Start configuring High Receivers')
+    foundError = foundError | set_rcuh_to_default_config(devices)
+    if not silent:
+        plib.log('Start configuring Station Digital Processor')
+    foundError = foundError | set_sdp_to_default_config(devices)
+    if not silent:
+        plib.log('Start configuring Subband Statistics')
+    foundError = foundError | set_sdp_sst_to_default_config(devices)
+    if not silent:
+        plib.log('Start configuring Crosslet Statistics')
+    foundError = foundError | set_sdp_xst_to_default_config(devices)
+    if not silent:
+        plib.log('Start configuring Beamlet Statistics')
+    foundError = foundError | set_sdp_bst_to_default_config(devices)
+    if not silent:
+        plib.log('Check setting to default configuration')
+    foundError = foundError |\
+        check_set_to_default_configuration(devices,
+                                           readonly=False)  # only at first time
+    # plot statistics data
+    if not silent:
+        plib.log('Plot statistics data')
+    plot_statistics(devices)
+    return foundError
+
+
+def initialise(devices, timeout=90, DEBUG=False):
+    '''
+    Initialise the system
+
+    Input arguments:
+    timeout = timeout to flag an error in seconds
+
+    Output arguments:
+    foundError
+    '''
+    devices = devices_list_2_dict(devices)
+    boot = devices['boot']
+    recv = devices['recv']
+    sst = devices['sst']
+    unb2 = devices['unb2']
+    # initialise boot device
+    boot.put_property({"DeviceProxy_Time_Out": 60})
+    boot.off()
+    boot.initialise()
+    boot.on()
+    # reboot system
+    # increase timeouts for that
+    if DEBUG:
+        print('Time out settings prior to increase for reboot:')
+        print(recv.get_timeout_millis())
+        print(unb2.get_timeout_millis())
+        print(sst.get_timeout_millis())
+    recv.set_timeout_millis(30000)
+    unb2.set_timeout_millis(60000)
+    sst.set_timeout_millis(20000)
+    unb2.put_property({"UNB2_On_Off_timeout" : 20})
+    recv.put_property({"RCU_On_Off_timeout" : 60})
+    # and start reboot process
+    boot.reboot()
+    print("Rebooting now")
+    t_start = time.time()
+    dt = time.time() - t_start
+    # indicate progress
+    while boot.booting_R and (dt < timeout):
+        if dt % 5 < 1:  # print every 5 seconds
+            print(f"Initialisation at {boot.progress_R}%: {boot.status_R}")
+        time.sleep(1)
+        dt = time.time() - t_start
+    print(f'Initialisation took {dt} seconds')
+    # check for errors
+    foundError = False
+    if boot.booting_R:
+        foundError = True
+        print(f"Warning! Initialisation still ongoing after timeout ({timeout} sec)")
+    if boot.uninitialised_devices_R:
+        foundError = True
+        print(f"Warning! Did not initialise {boot.uninitialised_devices_R}.")
+    return foundError
+
+
+def check_initialise(devices):
+    '''
+    Check the initialisation of the system by an independent monitoring
+    point.
+
+    Note that this method should be used after calling initialise()
+
+    All checks that are directly related to the boot device are already
+    in the method initialise()
+    '''
+    foundError = False
+    foundError = foundError | check_firmware_loaded(devices)
+    foundError = foundError | check_clock_locked(devices)
+    return foundError
+
+
+def check_firmware_loaded(devices, DEBUG=True):
+    '''
+    Verify that the firmware has loaded in the fpgas
+    Returns True if not loaded, False otherwise
+    '''
+    devices = devices_list_2_dict(devices)
+    sdp = devices['sdp']
+    res = sdp.FPGA_boot_image_R
+    if DEBUG:
+        print('')
+        print('Checking sdp.FPGA_boot_image_R:')
+        print('Reply: %s' % res)
+    if not np.array_equal(res[0:n_fgpa], np.ones(n_fgpa)):
+        print('Firmware is not loaded!')
+        return True
+    else:
+        if DEBUG:
+            print('Ok - Firmware loaded')
+        return False
+
+
+def check_clock_locked(devices, DEBUG=True):
+    '''
+    Verify that the APSCT clock is locked
+    Returns True if not loaded, False otherwise
+    '''
+    devices = devices_list_2_dict(devices)
+    apsct = devices['apsct']
+    res = apsct.APSCT_PLL_200MHz_error_R
+    if DEBUG:
+        print('')
+        print('Checking apsct.APSCT_PLL_200MHz_error_R:')
+        print('Reply: %s' % res)
+    if res:
+        print("APSCT clock not locked!")
+        return True
+    else:
+        if DEBUG:
+            print('Ok - APSCT clock locked')
+        return False
+
+
+def set_rcul_to_default_config(devices):
+    '''
+    Set RCU Low to its default configuration
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    #
+    # Set RCUs Lows in default modes
+    #
+
+    recv.set_defaults()
+
+    # Set attenuator, antenna power etc.
+    recv.RCU_band_select_RW = set_ant_masked(recv, LBAMask, recv.RCU_band_select_R,
+                                             [[1]*3]*32)  # 1 = 30-80 MHz
+    recv.RCU_PWR_ANT_on_RW = set_ant_masked(recv, LBAMask, recv.RCU_PWR_ANT_on_R,
+                                            [[True]*3]*32)  # Off
+    recv.RCU_attenuator_dB_RW = set_ant_masked(recv, LBAMask, recv.RCU_attenuator_dB_R,
+                                               [[0]*3]*32)  # 0dB attenuator
+    wait_receiver_busy(recv)
+
+    #
+    foundError = False
+    return foundError
+
+
+def set_rcuh_to_default_config(devices):
+    '''
+    Set RCU High to its default configuration
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    #
+    # Setup HBA RCUs and HBA Tiles
+    #
+
+    # Set RCU in correct mode
+    rcu_modes = recv.RCU_band_select_R
+    rcu_modes[8] = [2, 2, 2]  # 2 = 110-180 MHz, 1 = 170-230 MHz, 4 = 210-270 MHz ??
+    rcu_modes[9] = [2, 2, 2]
+    recv.RCU_band_select_RW = set_ant_masked(recv, hba_mask, recv.RCU_band_select_R,
+                                             rcu_modes)
+    wait_receiver_busy(recv)
+
+    # Switch on the Antenna
+    recv.RCU_PWR_ANT_on_RW = set_ant_masked(recv, hba_PwrMask, recv.RCU_PWR_ANT_on_R,
+                                            [[False]*3]*32)  # Switch off
+    wait_receiver_busy(recv)
+    recv.HBAT_PWR_LNA_on_RW = set_hba_mask(recv, hba_CtlMask, recv.HBAT_PWR_LNA_on_R,
+                                           [[True]*32]*96)  # LNA default on
+    recv.HBAT_BF_delay_steps_RW = set_hba_mask(recv, hba_CtlMask,
+                                               recv.HBAT_BF_delay_steps_R,
+                                               [[0]*32]*96)  # Default
+    recv.HBAT_PWR_on_RW = set_hba_mask(recv, hba_CtlMask, recv.HBAT_PWR_on_R,
+                                       [[True]*32]*96)  # Default
+    wait_receiver_busy(recv)
+    recv.RCU_PWR_ANT_on_RW = set_ant_masked(recv, hba_PwrMask, recv.RCU_PWR_ANT_on_R,
+                                            [[True]*3]*32)  # Switch on
+    wait_receiver_busy(recv)
+    # default in the tile: after power-on, zero delay
+    # delays should be set by the tile beam device
+    # recv.HBAT_BF_delay_steps_RW=[[1]*32]*32
+    # wait_receiver_busy(recv)
+    # by default: equal delay settings for all elements --> Pointing to zenith
+
+    # TODO: read values and only update the ones that need to be changed.
+    # Could be a function that manages the masks.
+    print(recv.RCU_PWR_ANT_on_RW)
+    print()
+    print("done")
+
+    #
+    foundError = False
+    return foundError
+
+####
+# RCU specific functions and variables for DTS-Outside
+####
+
+# General RCU functions
+
+
+def wait_receiver_busy(recv):
+    time.sleep(0.5)
+    while recv.RECVTR_translator_busy_R:
+        time.sleep(0.1)
+
+
+def set_ant_masked(recv, mask, old_value, new_value):
+    recv.ANT_mask_RW = mask
+
+    for imask, maski in enumerate(mask):
+        for rcu_input_index in range(3):
+            if maski[rcu_input_index]:
+                old_value[imask][rcu_input_index] = new_value[imask][rcu_input_index]
+# this used to be the following one liner, but that line ot too long.
+# Here is it, but cut in pieces:
+# old_value[imask] =
+# [new_value[imask][j] if maski[j] else old_value[imask][j] for j in range(3)]
+    return old_value
+
+
+def set_hba_mask(recv, mask, old_value, new_value):
+    recv.ANT_mask_RW = mask
+    for imask, maski in enumerate(mask):
+        for i in range(3):
+            if maski[i]:
+                old_value[imask*3 + i] = new_value[imask*3 + i]
+            else:
+                old_value[imask*3 + i] = old_value[imask*3 + i]
+    return old_value
+
+
+# LBA masks
+RCU2LMask = [True if rcu_nr in RCUs.RCU2L else False for rcu_nr in range(32)]
+LBAMask = [[True]*3 if rcu_nr in RCUs.RCU2L else [False]*3 for rcu_nr in range(32)]
+
+# HBA masks
+# signal inputs connected to HBAT power:
+si_hbat_pwr = [si*3 + x for si in RCUs.RCU2Hpwr for x in range(3)]
+# signal inputs. HBAT control connect to first element:
+si_nbat_ctl = [si*3 + x for si in RCUs.RCU2Hctl for x in range(3)]
+
+rcu_indices = range(32)
+hba_rcu_mask_pwr = [True if ri in RCUs.RCU2Hpwr else False for ri in rcu_indices]
+hba_rcu_mask_ctrl = [True if ri in RCUs.RCU2Hctl else False for ri in rcu_indices]
+hba_PwrMask = [[True]*3 if ri in RCUs.RCU2Hpwr else [False]*3 for ri in rcu_indices]
+hba_CtlMask = [[True]*3 if ri in RCUs.RCU2Hctl else [False]*3 for ri in rcu_indices]
+hba_mask = [[True]*3 if ri in RCUs.RCU2H else [False]*3 for ri in rcu_indices]
+
+####
+# Till here RCU specific functions and variables for DTS-Outside
+####
+
+
+def set_sdp_to_default_config(devices):
+    '''
+    Set RCU High to its default configuration
+    '''
+    devices = devices_list_2_dict(devices)
+    sdp = devices['sdp']
+    #
+    # Set SDP in default modes
+    #
+
+    plib.log('Start configuring SDP to default')
+    sdp.set_defaults()
+    # should be part of sdp.set_defaults:
+    Next_ring = [False]*16
+    Next_ring[3] = True
+    sdp.FPGA_ring_use_cable_to_previous_rn_RW = [True] + [False]*15
+    sdp.FPGA_ring_use_cable_to_next_rn_RW = Next_ring
+    sdp.FPGA_ring_nof_nodes_RW = [n_fgpa]*16
+    sdp.FPGA_ring_node_offset_RW = [0]*16
+
+    #
+    foundError = False
+    return foundError
+
+
+def set_sdp_sst_to_default_config(devices):
+    '''
+    Set SSTs to default
+    '''
+    devices = devices_list_2_dict(devices)
+    sst = devices['sst']
+
+    sst.set_defaults()
+    # prepare for subband stati
+    sst.FPGA_sst_offload_weighted_subbands_RW = [True]*16  # should be in set_defaults()
+
+    #
+    foundError = False
+    return foundError
+
+
+def set_sdp_xst_to_default_config(devices):
+    '''
+    Set XSTs to default
+    '''
+    devices = devices_list_2_dict(devices)
+    xst = devices['xst']
+
+    # prepare for correlations
+    int_time = 1  # used for 'medium' measurement using statistics writer
+    subband_step_size = 7  # 0 is no stepping
+    nofXST = 7
+    subband_select = [subband_step_size, 0, 1, 2, 3, 4, 5, 6]
+
+    xst.set_defaults()
+    # should be part of set_defaults()
+    xst.FPGA_xst_processing_enable_RW = [False]*16
+    # this (line above) should be first, then configure (and enable in user script)
+    xst.fpga_xst_ring_nof_transport_hops_RW = [3]*16  # [((n_fgpa/2)+1)]*16
+
+    Crosslets = [0]*16
+    for fpga_nr in range(n_fgpa):
+        Crosslets[fpga_nr] = nofXST
+    xst.FPGA_xst_offload_nof_crosslets_RW = Crosslets
+
+    xst.FPGA_xst_subband_select_RW = [subband_select]*16
+    xst.FPGA_xst_integration_interval_RW = [int_time]*16
+
+    #
+    foundError = False
+    return foundError
+
+
+def set_sdp_bst_to_default_config(devices):
+    '''
+    Set XSTs to default
+    '''
+    devices = devices_list_2_dict(devices)
+    beamlet = devices['beamlet']
+
+    # prepare for beamformer
+    beamlet.set_defaults()
+
+    beamlet.FPGA_bf_ring_nof_transport_hops_RW = [[1]]*3 + [[0]]*13
+    # line above should be part of set_defaults(), specifically for DTS-Outside
+
+    beamletoutput_mask = [False]*16
+    beamletoutput_mask[3] = True
+    # beamlet.FPGA_beamlet_output_enable_RW=beamletoutput_mask
+    beamlet.FPGA_beamlet_output_hdr_eth_destination_mac_RW = ["40:a6:b7:2d:4f:68"]*16
+    beamlet.FPGA_beamlet_output_hdr_ip_destination_address_RW = ["192.168.0.249"]*16
+
+    # define weights
+    # dim: FPGA, input, subband
+    xx = beamlet.FPGA_bf_weights_xx_R.reshape([16, 6, 488])
+    yy = beamlet.FPGA_bf_weights_yy_R.reshape([16, 6, 488])
+    # make a beam for the HBA
+    xx[:][:] = 0
+    yy[:][:] = 0
+    for si in [0, 1]:  # range(8*3,9*3+3):
+        unb2i = si//12
+        si2 = si//2
+        si2 = si2 % 6
+        if si % 2 == 0:
+            xx[unb2i, si2] = 16384
+        else:
+            yy[unb2i, si2] = 16384
+
+    beamlet.FPGA_bf_weights_xx_RW = xx.reshape([16, 6*488])
+    beamlet.FPGA_bf_weights_yy_RW = yy.reshape([16, 6*488])
+
+    blt = beamlet.FPGA_beamlet_subband_select_R.reshape([16, 12, 488])
+    blt[:] = np.array(range(10, 10 + 488))[np.newaxis, np.newaxis, :]
+    beamlet.FPGA_beamlet_subband_select_RW = blt.reshape([16, 12*488])
+
+    #
+    foundError = False
+    return foundError
+
+
+def enable_outputs(devices, enable_sst=True, enable_xst=True, enable_bst=True,
+                   enable_beam_output=True):
+    '''
+    Enable the station outputs after everything is set.
+    '''
+    devices = devices_list_2_dict(devices)
+    if enable_beam_output:
+        sdp = devices['sdp']
+        sdp.FPGA_processing_enable_RW = [True]*16
+    if enable_sst:
+        sst = devices['sst']
+        sst.FPGA_sst_offload_enable_RW = [True]*16
+    if enable_xst:
+        xst = devices['xst']
+        xst.FPGA_xst_processing_enable_RW = [True]*16
+    if enable_bst:
+        beamlet = devices['beamlet']
+        beamlet.FPGA_beamlet_output_enable_RW = [i == 3 for i in range(16)]
+
+    #
+    foundError = False
+    return foundError
+
+
+def check_set_to_default_configuration(devices, readonly=True):
+    '''
+    Check the default configuration of the system to be set correctly.
+    '''
+    foundError = False
+    foundError = foundError | check_default_configuration_low_receivers(devices)
+    foundError = foundError | check_default_configuration_hbat(devices,
+                                                               readonly=readonly)
+    foundError = foundError | check_rcu_fpga_interface(devices, si_to_check=range(18))
+    foundError = foundError | check_set_sdp_xst_to_default_config(devices,
+                                                                  lane_mask=range(0, 6))
+
+    print_fpga_input_statistics(devices)
+
+    return foundError
+
+
+def check_default_configuration_low_receivers(devices, mask=LBAMask, DEBUG=True):
+    '''
+    Check that the LNAs are on
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    foundError = False
+
+    res = recv.RCU_PWR_ANT_on_R
+    if len(mask) > 0:
+        res = res[mask]
+    if DEBUG:
+        print('')
+        print('Checking recv.RCU_PWR_ANT_on_R:')
+        print('Reply%s: %s' % (' (masked)' if len(mask) > 0 else '', res))
+    if not (res).all():
+        print('One or more LBAs are powered off!')
+        foundError = True
+    else:
+        if DEBUG:
+            print('Ok - Low Receivers operational')
+    return foundError
+
+
+def check_default_configuration_hbat(devices, readonly=True):
+    '''
+    Check that the HBA Tiles are on
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    HBATsi = [rcu*3 + r_input for rcu in RCUs.RCU2Hpwr for r_input in range(3)]
+    if not readonly:
+        recv.RECVTR_monitor_rate_RW = 1  # for this checking
+        time.sleep(3)
+
+    foundError = False
+    foundError = foundError | check_rcu_vin(devices, si_mask=HBATsi)
+    foundError = foundError | check_hbat_vout(devices, si_mask=HBATsi)
+    foundError = foundError | check_hbat_iout(devices, si_mask=HBATsi)
+
+    if not readonly:
+        recv.RECVTR_monitor_rate_RW = 10
+
+    # TODO - what is checked for here?
+    print(recv.HBAT_PWR_LNA_on_R[8:10])
+    print(recv.ANT_mask_RW[8:10])
+
+    return foundError
+
+
+def check_rcu_vin(devices, si_mask=[], valid_range=[30, 50], DEBUG=True):
+    '''
+    Verify that the HBA Tile is powered on by comparing the V_in monitoring
+    point on the RCU against the thresholds for the given signal inputs.
+
+    If no si_mask is given, all returned values are compared.
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+
+    res = recv.RCU_PWR_ANT_VIN_R
+    if len(si_mask) > 0:
+        res = [res[si//3][si % 3] for si in si_mask]
+    if DEBUG:
+        print('')
+        print('Checking recv.RCU_PWR_ANT_VIN_R:')
+        print('Reply%s: %s' % (' (masked)' if len(si_mask) > 0 else '', res))
+    if not (np.all(valid_range[0] < np.array(res)) and
+            np.all(np.array(res) < valid_range[1])):
+        print("No power at input RCU2!")
+        return True
+    else:
+        if DEBUG:
+            print('Ok - Power at input RCU2')
+        return False
+
+
+def check_hbat_vout(devices, si_mask=[], valid_range=[30, 50], DEBUG=True):
+    '''
+    Verify that the HBA Tile is powered on by comparing the V_out monitoring
+    point on the RCU against the thresholds for the given signal inputs.
+
+    If no si_mask is given, all returned values are compared.
+    '''
+    # Vout should be 48 +- 1 V
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+
+    res = recv.RCU_PWR_ANT_VOUT_R
+    if len(si_mask) > 0:
+        res = [res[si//3][si % 3] for si in si_mask]
+    if DEBUG:
+        print('')
+        print('Checking recv.RCU_PWR_ANT_VOUT_R:')
+        print('Reply%s: %s' % (' (masked)' if len(si_mask) > 0 else '', res))
+    if not (np.all(valid_range[0] < np.array(res)) and
+            np.all(np.array(res) < valid_range[1])):
+        print("Not all HBA Tiles powered on!")
+        return True
+    else:
+        if DEBUG:
+            print('Ok - HBA Tiles powered on')
+        return False
+
+
+def check_hbat_iout(devices, si_mask=[], valid_range=[0.6, 0.8], DEBUG=True):
+    '''
+    Verify that the HBA Tile is powered on by comparing the I_out monitoring
+    point on the RCU against the thresholds for the given signal inputs.
+
+    If no si_mask is given, all returned values are compared.
+    '''
+    # Iout = 0.7 +- 0.1 Amp
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+
+    res = recv.RCU_PWR_ANT_IOUT_R
+    if len(si_mask) > 0:
+        res = [res[si//3][si % 3] for si in si_mask]
+    if DEBUG:
+        print('')
+        print('Checking recv.RCU_PWR_ANT_IOUT_R:')
+        print('Reply%s: %s' % (' (masked)' if len(si_mask) > 0 else '', res))
+    if not (np.all(valid_range[0] < np.array(res)) and
+            np.all(np.array(res) < valid_range[1])):
+        print("Not all HBA Tiles powered on (draw current)!")
+        return True
+    else:
+        if DEBUG:
+            print('Ok - HBA Tiles powered on (draw current)')
+        return False
+
+
+def check_rcu_fpga_interface(devices, si_to_check=range(18)):
+    '''
+    Function to check the RCU - FPGA interface
+
+    # if error occurs, try to reboot by resetting Uniboard (switch)
+    # if only one works than the APSCT board is not making contact with the backplane.
+    # (Can take up to 3 times to get it right..)
+    # 23/6/2022T15:42 errorcode 0x4h: xxxxx on all of them.
+    # Error descriptions can be found here:
+    # https://support.astron.nl/jira/browse/L2SDP-774
+    '''
+    devices = devices_list_2_dict(devices)
+    sdp = devices['sdp']
+    foundError = False
+    # TODO: the number of calls can be reduced
+    for si in si_to_check:
+        lock = sdp.FPGA_jesd204b_csr_dev_syncn_R[si//12][si % 12]
+        if lock != 1:
+            foundError = True
+            print(f"Transceiver SI: {si} is not locked {lock}")
+        err0 = sdp.FPGA_jesd204b_rx_err0_R[si//12][si % 12]
+        if err0 != 0:
+            foundError = True
+            print(f"Si {si} has JESD err0: 0x{err0:x}h")
+        err1 = sdp.FPGA_jesd204b_rx_err1_R[si//12][si % 12]
+        if ((err1 & 0x02FF) != 0):
+            foundError = True
+            print(f"Si {si} has JESD err1: 0x{err1:x}h")
+    return foundError
+
+
+def check_set_sdp_xst_to_default_config(devices, lane_mask=[], t1=None, dt=2.5,
+                                        DEBUG=True):
+    '''
+    Check sdp/xst configuration by checking the returned timestamps for:
+    1. all being the same
+    2. not epoch zero (at 1 Jan 1970 00:00:00)
+    3. within 5 seconds from the local clock timestamp (now +/- 2.5 sec)
+
+    Input arguments:
+    t1 = number of seconds since the epoch to compare against
+         Default: None --> now()
+    dt = delta in seconds around t1, in which the xst timestamps should fall.
+         Default: 2.5 --> p/m 2.5 sec
+    '''
+    devices = devices_list_2_dict(devices)
+    xst = devices['xst']
+    if t1 is None:
+        t1 = datetime.datetime.now().timestamp()
+    res = xst.xst_timestamp_R
+    if len(lane_mask) > 0:  # mask lanes
+        res = [res[lane] for lane in lane_mask]
+    if DEBUG:
+        print('')
+        print('Checking xst.xst_timestamp_R:')
+        print('Reply%s: %s' % (' (masked)' if len(lane_mask) > 0 else '', res))
+        print('Converted:')
+        print("\n".join([f"{time.asctime(time.gmtime(xst_time))}" for xst_time in res]))
+    if not (res == res[0]).all():
+        print("Not all xst timestamps are equal!")
+        return True
+    if res[0] == 0:
+        print("One or more xst timestamps are zero!")
+        return True
+    if abs(res[0] - t1) > dt:
+        print("One or more xst timestamps are outside the specified time window!")
+        return True
+    else:
+        if DEBUG:
+            print('Ok - xst configured correctly')
+        return False
+
+
+def print_fpga_input_statistics(devices):
+    '''
+    Print FPGA input statistics
+    '''
+    devices = devices_list_2_dict(devices)
+    sdp = devices['sdp']
+
+    # get data from device
+    dc_offset = sdp.FPGA_signal_input_mean_R
+    signal_input_rms = sdp.FPGA_signal_input_rms_R
+    lock_to_adc = sdp.FPGA_jesd204b_csr_dev_syncn_R
+    err0 = sdp.FPGA_jesd204b_rx_err0_R
+    err1 = sdp.FPGA_jesd204b_rx_err1_R
+    count = sdp.FPGA_jesd204b_csr_rbd_count_R
+    signal_index = np.reshape(np.arange(count.shape[0]*count.shape[1]), count.shape)
+    # get some info on the hardware connections
+    pn_index = plib.signal_index_2_processing_node_index(signal_index)
+# processing_node_index = (signal_index%192)//12
+    pn_input_index = plib.signal_index_2_processing_node_input_index(signal_index)
+# processing_node_input_index = (signal_index%192)%12
+    rcu_index = plib.signal_index_2_rcu_index(signal_index)
+# rcu_index = (signal_index//6)*2 + signal_index%2
+    rcu_input_index = plib.signal_index_2_rcu_input_index(signal_index)
+# rcu_input_index = (signal_index//2)%3
+    # print
+    print('+------------+-------------+----------------------------------------------+')
+    print('|   Input    | Processed by|             Statistics of input              |')
+    print('+---+---+----+------+------+---------+---------+-----+------+------+------+')
+    print('|RCU|RCU| SI | Node | Node | DC      | RMS     | Lock| Err0 | Err1 | Count|')
+    print('|idx|inp|    | Index| Input| Offset  |   (LSB) |     | (0x) | (0x) | (0x) |')
+    print('|   |idx|    |      | Index|   (LSB) |         |     |      |      |      |')
+    print('+---+---+----+------+------+---------+---------+-----+------+------+------+')
+    pfmt = '| %02d| %1d |%3d |  %02d  |  %2d  | %7.2f | %7.2f | %3s | %4x | %4x | %4x |'
+    for ri, rii, si, pni, pnii, dco, rms, lock, e0, e1, cnt in zip(
+      rcu_index.flatten(), rcu_input_index.flatten(), signal_index.flatten(),
+      pn_index.flatten(), pn_input_index.flatten(),
+      dc_offset.flatten(), signal_input_rms.flatten(), lock_to_adc.flatten(),
+      err0.flatten(), err1.flatten(), count.flatten()):
+        print(pfmt % (ri, rii, si, pni, pnii, dco, rms, lock, e0, e1, cnt))
+    print('+---+---+----+------+------+---------+---------+-----+------+------+------+')
+
+
+def report_setup_configuration(devices, save_fig=True, DEBUG=False):
+    '''
+    Read name and version data from the setup and print that as a report
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    apsct = devices['apsct']
+    apspu = devices['apspu']
+    unb2 = devices['unb2']
+    sdp = devices['sdp']
+    # Receiver / High Band Antenna Tile
+
+    # APS Hardware
+    # APS / Receiver Unit2 side
+    RCU2_PCB_ids = recv.RCU_PCB_ID_R.tolist()
+    RCU2_PCB_nrs = list(recv.RCU_PCB_number_R)
+    RCU2_PCB_vrs = list(recv.RCU_PCB_version_R)
+
+    # APS / UniBoard2 side
+    # TODO, FIXME - shouldn't this be a list with multiple entries?
+    # See also LOFAR2-10774 (APSCT)
+    APSCT_PCB_ids = [apsct.APSCT_PCB_ID_R]
+    APSCT_PCB_nrs = [apsct.APSCT_PCB_number_R]
+    APSCT_PCB_vrs = [apsct.APSCT_PCB_version_R]
+    # TODO, FIXME - shoudln't this be a list with multiple entries?
+    # See also LOFAR2-11434 (APSPU)
+    APSPU_PCB_ids = [apspu.APSPU_PCB_ID_R]
+    APSPU_PCB_nrs = [apspu.APSPU_PCB_number_R]
+    APSPU_PCB_vrs = [apspu.APSPU_PCB_version_R]
+    UNB2_PCB_ids = unb2.UNB2_PCB_ID_R.tolist()
+    UNB2_PCB_nrs = list(unb2.UNB2_PCB_number_R)
+    UNB2_PCB_vrs = list(sdp.FPGA_hardware_version_R)
+    # TODO, FIXME:
+    # FPGA_hardware_version_R returns hardware version info per FPGA. Hardware is shared
+    # per 4 fpgas, so why repeating values here?
+    UNB2_PCB_vrs = [UNB2_PCB_vrs[idx] for idx in [0, 4, 8, 12]]
+
+    # TODO, FIXME:
+    # fixing list lengths for a large station with multiple APS
+    # this makes sure that the list indexing below works
+    n_aps = 4
+    APSCT_PCB_ids = APSCT_PCB_ids + [0]*(n_aps-len(APSCT_PCB_ids))
+    APSCT_PCB_nrs = APSCT_PCB_nrs + ['']*(n_aps-len(APSCT_PCB_nrs))
+    APSCT_PCB_vrs = APSCT_PCB_vrs + ['']*(n_aps-len(APSCT_PCB_vrs))
+    APSPU_PCB_ids = APSPU_PCB_ids + [0]*(n_aps*2-len(APSPU_PCB_ids))
+    APSPU_PCB_nrs = APSPU_PCB_nrs + ['']*(n_aps*2-len(APSPU_PCB_nrs))
+    APSPU_PCB_vrs = APSPU_PCB_vrs + ['']*(n_aps*2-len(APSPU_PCB_vrs))
+    UNB2_PCB_ids = UNB2_PCB_ids + [0]*(n_aps*2-len(UNB2_PCB_ids))
+    UNB2_PCB_nrs = UNB2_PCB_nrs + ['']*(n_aps*2-len(UNB2_PCB_nrs))
+    UNB2_PCB_vrs = UNB2_PCB_vrs + ['']*(n_aps*2-len(UNB2_PCB_vrs))
+
+    if DEBUG:
+        print('APSPU:')
+        print('APSCT_PCB_ids: %s' % APSCT_PCB_ids)
+        print('APSCT_PCB_nrs: %s' % APSCT_PCB_nrs)
+        print('APSCT_PCB_vrs: %s' % APSCT_PCB_vrs)
+        print('APSPU:')
+        print('APSPU_PCB_ids: %s' % APSPU_PCB_ids)
+        print('APSPU_PCB_nrs: %s' % APSPU_PCB_nrs)
+        print('APSPU_PCB_vrs: %s' % APSPU_PCB_vrs)
+        print('UNB2:')
+        print('UNB2_PCB_ids: %s' % UNB2_PCB_ids)
+        print('UNB2_PCB_nrs: %s' % UNB2_PCB_nrs)
+        print('UNB2_PCB_vrs: %s' % UNB2_PCB_vrs)
+
+    APS_PCB_ids = [APSPU_PCB_ids[1], UNB2_PCB_ids[0],
+                   APSCT_PCB_ids[0], APSPU_PCB_ids[0],
+                   UNB2_PCB_ids[1]]
+    APS_PCB_nrs = [APSPU_PCB_nrs[1], UNB2_PCB_nrs[0],
+                   APSCT_PCB_nrs[0], APSPU_PCB_nrs[0],
+                   UNB2_PCB_nrs[1]]
+    APS_PCB_vers = [APSCT_PCB_vrs[1], UNB2_PCB_vrs[0],
+                    APSCT_PCB_vrs[0], APSCT_PCB_vrs[0],
+                    UNB2_PCB_vrs[1]]
+
+    if not save_fig:
+        fig_filenames = []
+    else:
+        fig_filenames = [f'{plot_dir}/{plib.get_timestamp("filename")}_aps_config.png',
+                         f'{plot_dir}/{plib.get_timestamp("filename")}_aps_config.pdf']
+    plot_aps_configuration(RCU2_PCB_ids, RCU2_PCB_nrs,
+                           RCU2_PCB_vrs, APS_PCB_ids, APS_PCB_nrs,
+                           APS_PCB_vers, filenames=fig_filenames)
+
+    # Services
+    # Firmware images
+    firmware_versions = list(sdp.FPGA_firmware_version_R)
+    # make sure there are versions for at least n_aps subracks
+    firmware_versions = firmware_versions + ['']*(n_aps*8-len(firmware_versions))
+    aps_services = {}
+    for apsIdx, fpgaIdx in zip(sorted([i for i in range(n_aps)]*8), range(n_aps*8)):
+        aps_services['aps%i-fpga%02d' % (apsIdx, fpgaIdx)] = firmware_versions[fpgaIdx]
+        aps_services['aps%i-unb2tr' % apsIdx] = ''  # TODO
+        aps_services['aps%i-recvtr' % apsIdx] = ''  # TODO
+        aps_services['aps%i-apscttr' % apsIdx] = ''  # TODO
+        aps_services['aps%i-apsputr' % apsIdx] = ''  # TODO
+    # Translators
+    sdptr_software_versions = [sdp.TR_software_version_R]
+    sdptr_software_versions = sdptr_software_versions +\
+        ['']*(2 - len(sdptr_software_versions))
+    station_services = {}
+    station_services['sdptr0'] = sdptr_software_versions[0]
+    station_services['sdptr1'] = sdptr_software_versions[1]
+    station_services['ccdtr'] = ''  # TODO, currently missing
+    # station control tango devices
+    for dev in devices:
+        station_services[devices[dev].name().lower()] = devices[dev].version_R
+    if not save_fig:
+        figs = []
+    else:
+        figs = [f'{plot_dir}/{plib.get_timestamp("filename")}_station_services.png',
+                f'{plot_dir}/{plib.get_timestamp("filename")}_station_services.pdf']
+    plot_services_configuration(aps_services, station_services, filenames=figs)
+
+
+def plot_aps_configuration(RCU_PCB_identifiers, RCU_PCB_numbers,
+                           RCU_PCB_versions, APS_PCB_ids,
+                           APS_PCB_nrs, APS_PCB_vers,
+                           filenames=[f'{plot_dir}/aps_config.png',
+                                      f'{plot_dir}/aps_config.pdf']):
+
+    fig, axs = plt.subplots(1, 2, figsize=(20, 8))
+
+    plot_aps_rcu2_configuration(axs[0], RCU_PCB_identifiers,
+                                RCU_PCB_numbers, RCU_PCB_versions)
+    plot_aps_unb2_configuration(axs[1], APS_PCB_ids,
+                                APS_PCB_nrs, APS_PCB_vers)
+    for filename in filenames:
+        fig.savefig(filename)
+    fig.show()
+
+
+def plot_aps_rcu2_configuration(ax, RCU_PCB_identifiers, RCU_PCB_numbers,
+                                RCU_PCB_versions,
+                                RCU_locations=aps_location_labels,
+                                RCU_pos_x=aps_rcu_pos_x,
+                                RCU_pos_y=aps_rcu_pos_y):
+    ax.plot(0, 0)
+    dx = dy = 0.9
+    ax.set_xlim(-1 + dx, np.max(RCU_pos_x) + 1)
+    ax.set_ylim(-1 + dy, np.max(RCU_pos_y) + 1)
+    for x, y, loc, PCB_id, PCB_nr, PCB_v in zip(RCU_pos_x, RCU_pos_y,
+                                                RCU_locations,
+                                                RCU_PCB_identifiers,
+                                                RCU_PCB_numbers,
+                                                RCU_PCB_versions):
+        ax.plot([x, x + dx, x + dx, x, x], [y, y, y + dy, y + dy, y],
+                color='black', alpha=0.2)
+        ax.text(x+dx/2, y+dy, loc, va='top', ha='center')
+        ax.text(x+dx/2, y+dy/2,
+                '%s\n%s\n%s' % (PCB_v, PCB_id, PCB_nr),
+                va='center', ha='center', rotation=90)
+        if PCB_id == 0 and PCB_v == '' and PCB_nr == '':
+            ax.plot([x, x + dx], [y, y + dy], color='black', alpha=0.2)
+            ax.plot([x + dx, x], [y, y + dy], color='black', alpha=0.2)
+    ax.set_title('Configuration of RCU2s')
+    ax.set_xticks([])
+    ax.set_yticks([])
+
+
+def plot_aps_unb2_configuration(ax, APS_PCB_identifiers, APS_PCB_numbers,
+                                APS_PCB_versions,
+                                mod_locations=aps_location_labels,
+                                mod_pos_x=aps_mod_pos_x,
+                                mod_pos_y=aps_mod_pos_y):
+    ax.plot(0, 0)
+    dx = 0.9
+    dy = 3.9
+    ax.set_xlim(-1 + dx, np.max(aps_mod_pos_x) + 2)
+    ax.set_ylim(-1 + dx, np.max(aps_mod_pos_y) + 4)
+    for x, y, loc, PCB_id, PCB_nr, PCB_v in zip(mod_pos_x, mod_pos_y,
+                                                mod_locations,
+                                                APS_PCB_identifiers,
+                                                APS_PCB_numbers,
+                                                APS_PCB_versions):
+        ax.plot([x, x + dx, x + dx, x, x], [y, y, y + dy, y + dy, y],
+                color='black', alpha=0.2)
+        ax.text(x+dx/2, y+dy, loc, va='top', ha='center')
+        ax.text(x+dx/2, y+dy/2,
+                '%s\n%s\n%s' % (PCB_v, PCB_id, PCB_nr),
+                va='center', ha='center', rotation=90)
+        if PCB_id == 0 and PCB_v == '' and PCB_nr == '':
+            ax.plot([x, x + dx], [y, y + dy], color='black', alpha=0.2)
+            ax.plot([x + dx, x], [y, y + dy], color='black', alpha=0.2)
+    ax.set_title('Configuration of UniBoard2 side')
+    ax.set_xticks([])
+    ax.set_yticks([])
+
+
+def plot_services_configuration(aps_services, station_services,
+                                filenames=[f'{plot_dir}/station_services.png',
+                                           f'{plot_dir}/station_services.pdf']):
+
+    fig, ax = plt.subplots(1, 1, figsize=(16, 15))
+
+    aps_service_dx = 0.9
+    aps_service_dy = 0.9
+    first_aps_service = sorted(aps_services)[0][5:]
+    y = y0 = len(aps_services)/2 + 2
+    for aps_service in sorted(aps_services):
+        x = 0 if (aps_service[3] in '01') else 1
+        y = y0 if (aps_service[5:] == first_aps_service and
+                   aps_service[3] in '02') else y - 1
+        y = y - 1 if (aps_service[5:] == first_aps_service and
+                      aps_service[3] in '13') else y
+        ax.plot([x, x+aps_service_dx, x+aps_service_dx, x, x],
+                [y, y, y+aps_service_dy, y+aps_service_dy, y],
+                color='black', alpha=0.2)
+        ax.text(x, y + aps_service_dy/2,
+                aps_service,
+                va='center', ha='left', color='black')
+        ax.text(x + aps_service_dx, y + aps_service_dy/2,
+                aps_services[aps_service],
+                va='center', ha='right', color='blue')
+        if aps_services[aps_service] == '':
+            ax.plot([x, x+aps_service_dx], [y, y+aps_service_dy],
+                    color='black', alpha=0.2)
+            ax.plot([x+aps_service_dx, x], [y, y+aps_service_dy],
+                    color='black', alpha=0.2)
+
+    stat_service_dx = 1.9
+    stat_service_dy = 0.9
+    x = 0
+    y = 0
+    for stat_service in sorted(station_services):
+        ax.plot([x, x+stat_service_dx, x+stat_service_dx, x, x],
+                [y, y, y+stat_service_dy, y+stat_service_dy, y],
+                color='black', alpha=0.2)
+        ax.text(x, y + stat_service_dy/2,
+                stat_service,
+                va='center', ha='left', color='black')
+        ax.text(x + stat_service_dx, y + stat_service_dy/2,
+                station_services[stat_service],
+                va='center', ha='right', color='blue')
+        if station_services[stat_service] == '':
+            ax.plot([x, x+stat_service_dx], [y, y+stat_service_dy],
+                    color='black', alpha=0.2)
+            ax.plot([x+stat_service_dx, x], [y, y+stat_service_dy],
+                    color='black', alpha=0.2)
+        y -= 1
+    ax.set_title('Services in setup')
+    ax.set_xticks([])
+    ax.set_yticks([])
+
+    for filename in filenames:
+        fig.savefig(filename)
+    fig.show()
+
+
+def plot_statistics(devices, save_fig=True, silent=False):
+    if not silent:
+        plib.log('Plot subband statistics')
+    if not save_fig:
+        figs = []
+    else:
+        figs = [f'{plot_dir}/{plib.get_timestamp("filename")}_subband_statistics.png',
+                f'{plot_dir}/{plib.get_timestamp("filename")}_subband_statistics.pdf']
+    plot_subband_statistics(devices,
+                            filenames=figs)
+    # TODO:
+    # if not silent:
+    #     plib.log('Plot crosslet statistics')
+    # TODO:
+    # if not silent:
+    #     plib.log('Plot beamlet statistics')
+
+
+def plot_subband_statistics(devices, xlim=[-5, 305], ylim=[-105, -30],
+                            filenames=[f'{plot_dir}/subband_statistics.png',
+                                       f'{plot_dir}/subband_statistics.pdf']):
+    '''
+    Make the subband statistics plot
+    '''
+    sst_data, band_names, frequency_axes = get_subband_statistics(devices)
+
+    fig, axs = plt.subplots(1, 1, figsize=(18, 6))
+    for signal_index in range(len(sst_data)):
+        if band_names[signal_index] is None:
+            continue
+        freq = frequency_axes[band_names[signal_index]]
+        # plot data in dB full scale
+        plot_data = 10*np.log10(sst_data[signal_index, :] + 1) - 128 - 6 * 4
+        # since we're lacking antenna names,
+        # show the RCU input info as label
+        rcu_index, rcu_input_index = plib.gsi_2_rcu_index_and_rcu_input_index(
+            signal_index)
+        label = 'RCU%02d, Inp%d' % (rcu_index, rcu_input_index)
+        plt.plot(freq/1e6, plot_data, label=label)
+    plt.grid()
+    plt.legend(bbox_to_anchor=(1, 1), loc="upper left")
+    plt.xlim(xlim)
+    plt.ylim(ylim)
+    plt.xlabel('Frequency (MHz)')
+    plt.ylabel('Power (dB full scale)')
+    plt.suptitle('Subband Statistics for all fully-connected receivers')
+    plt.title(f'{station_name}; {plib.get_timestamp()}')
+
+    for filename in filenames:
+        fig.savefig(filename)
+    fig.show()
+
+
+def get_subband_statistics(devices, DEBUG=True):
+    '''
+    Get subband statistics data from the setup, including proper axes and labels
+    '''
+    devices = devices_list_2_dict(devices)
+    recv = devices['recv']
+    sst = devices['sst']
+
+    # get band names for the current configuration
+    cur_band_names = plib.get_band_names(recv.RCU_band_select_R,
+                                         # TODO - when correct LCNs are provided
+                                         # in 'RCU_PCB_number_R', uncomment the
+                                         # following line:
+                                         # recv.RCU_PCB_number_R,
+                                         # and remove the following line:
+                                         recv.RCU_PCB_ID_R,
+                                         LB_tags=RCU2L_TAGS,
+                                         HB_tags=RCU2H_TAGS)
+
+    # replace the aliases of LB1 and LB2 into LB: (same frequency axis)
+    for alias in plib.LB_aliases:
+        while alias in cur_band_names:
+            cur_band_names[cur_band_names.index(alias)] = 'LB'
+    # get frequency axes
+
+    # TODO - get frequency axes correct (indexing is incorrect/inconsequent)
+    key1 = 'HB2'
+    key2 = 'HB1'
+    tmp_key = 'magic_band'
+    while key1 in cur_band_names:
+        cur_band_names[cur_band_names.index(key1)] = tmp_key
+    while key2 in cur_band_names:
+        cur_band_names[cur_band_names.index(key1)] = key1
+    while tmp_key in cur_band_names:
+        cur_band_names[cur_band_names.index(tmp_key)] = key2
+
+    frequency_axes = {}
+    for band_name in plib.get_valid_band_names():
+        if band_name in cur_band_names:
+            frequency_axes[band_name] = plib.get_frequency_for_band_name(band_name)
+    # get sst data
+    sst_data = sst.sst_r
+    # make sure header data is of correct size
+    cur_band_names = cur_band_names + [None]*(len(sst_data) - len(cur_band_names))
+    return sst_data, cur_band_names, frequency_axes
+
+
+def devices_list_2_dict(devicesList,
+                        deviceKeys=['boot', 'unb2', 'recv', 'sdp', 'sst',
+                                    'bst', 'xst', 'digitalbeam', 'tilebeam',
+                                    'beamlet', 'apsct', 'apspu']):
+    '''
+    Convert a list of devices to a dictionary with known keys for internal
+    use.
+
+    Devices are selected based on substring presence in the devices name
+
+    Input arguments:
+    devicesList = list of Tango devices
+    deviceKeys  = list of keys, should be a substring of the device names
+
+    Output arguments:
+    devicesDict = dict of devices with known keys
+    '''
+    devicesDict = {}
+    if isinstance(devicesList, dict):
+        # if by accident devicesList is already a dict,
+        # then simply return the dict
+        return devicesList
+    for device in devicesList:
+        for deviceKey in deviceKeys:
+            if deviceKey in device.name().lower():
+                devicesDict[deviceKey] = device
+    return devicesDict
+
+
+def main():
+    return False
+
+
+if __name__ == '__main__':
+    exit(main())
diff --git a/lofar_station_client/processing_lib.py b/lofar_station_client/processing_lib.py
new file mode 100644
index 0000000000000000000000000000000000000000..6461e2af9ed8e2b2ce9b4174d43d06e628477d6c
--- /dev/null
+++ b/lofar_station_client/processing_lib.py
@@ -0,0 +1,414 @@
+# processing_lib.py: Module with generic functions for processing
+# LOFAR2 station data
+#
+# -*- coding: utf-8 -*-
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may
+# not use this file except in compliance with the License. You may obtain
+# a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations
+# under the License.
+
+
+"""  module processing_lib
+
+This module contains generic functions for processing LOFAR2 station data.
+
+"""
+
+# import h5py
+# from os import listdir
+# from os.path import isfile, join
+import datetime
+import numpy as np
+# from mpl_toolkits.mplot3d import axes3d
+# import matplotlib.pyplot as plt
+
+
+cst_n_sub = 512  # number of subbands as output of subband generation in firmware
+cst_fs = 100e6   # sampling frequency in Hz
+bandlims = {'LB': [0, 100e6], 'HB1': [100e6, 200e6], 'HB2': [200e6, 300e6]}
+LB_aliases = ['LB1', 'LB2']
+
+
+def signal_index_2_processing_node_index(signal_index):
+    '''
+    Convert signal index to processing node index.
+    Can be used to determine on which processing node the signal is processed.
+    '''
+    return (signal_index % 192)//12
+
+
+def signal_index_2_processing_node_input_index(signal_index):
+    '''
+    Convert signal index to processing node input index.
+    Can be used to determine wich input to the processing node is used for the signal.
+    '''
+    return (signal_index % 192) % 12
+
+
+def signal_index_2_uniboard_index(signal_index):
+    '''
+    Convert signal index to UniBoard index.
+    Can be used to determine on which UniBoard the signal is processed.
+    '''
+    return signal_index//48
+
+
+def signal_index_2_rcu_index(signal_index):
+    '''q
+    Convert signal index to RCU index.
+    Can be used to determine which RCU is used to provide the signal.
+    '''
+    return (signal_index//6)*2 + signal_index % 2
+
+
+def signal_index_2_rcu_input_index(signal_index):
+    '''
+    Convert signal index to RCU input index.
+    Can be used to determine which input to the RCU is used for the signal.
+    '''
+    return (signal_index//2) % 3
+
+
+def signal_index_2_aps_index(signal_index):
+    '''
+    Convert signal index to APS index.
+    Can be used to determine in which Antenna Processing Subrack this signal is
+    processed.
+    '''
+    return signal_index//96
+
+
+def signal_input_2_rcu_index_and_rcu_input_index(signal_index):
+    '''
+    Given the gsi, get the RCU index and the RCU input index in one call.
+    '''
+    return (signal_index_2_rcu_index(signal_index),
+            signal_index_2_rcu_input_index(signal_index))
+
+
+# the inverse:
+def rcu_index_2_signal_index(rcu_index):
+    '''
+    Convert RCU index to signal indices.
+    '''
+    return rcu_index_and_rcu_input_index_2_signal_index(rcu_index)
+
+
+def rcu_index_and_rcu_input_index_2_signal_index(rcu_index,
+                                                 rcu_input_index=np.array([0, 1, 2])):
+    '''
+    Convert RCU index and RCU input index to signal index.
+    '''
+    return rcu_index + rcu_input_index * 2 + (rcu_index//2) * 4
+
+
+# backward compatible methods:
+def gsi_2_rcu_input_index(gsi):
+    print('Warning! Do not use gsi_2_rcu_input_index().')
+    print('Use signal_index_2_rcu_input_index() instead')
+    return signal_index_2_rcu_input_index(gsi)
+
+
+def gsi_2_rcu_index(gsi):
+    print('Warning! Do not use gsi_2_rcu_index().')
+    print('Use signal_index_2_rcu_index() instead')
+    return signal_index_2_rcu_index(gsi)
+
+
+def gsi_2_rcu_index_and_rcu_input_index(gsi):
+    print('Warning! Do not use gsi_2_rcu_index_and_rcu_input_index().')
+    print('Use signal_input_2_rcu_index_and_rcu_input_index() instead')
+    return signal_input_2_rcu_index_and_rcu_input_index(gsi)
+
+
+def rcu_index_2_gsi(rcu_index):
+    print('Warning! Do not use rcu_index_2_gsi().')
+    print('Use rcu_index_2_signal_index() instead')
+    return rcu_index_2_signal_index(rcu_index)
+
+
+def rcu_index_and_rcu_input_index_2_gsi(rcu_index, rcu_input_index=np.array([0, 1, 2])):
+    print('Warning! Do not use rcu_index_and_rcu_input_index_2_gsi().')
+    print('Use rcu_index_and_rcu_input_index_2_signal_index() instead')
+    return rcu_index_and_rcu_input_index_2_signal_index(rcu_index, rcu_input_index)
+
+
+def test_rcu_index_2_gsi():
+    test_input = []
+    test_output = []
+    # test 1
+    test_input.append(np.array([0]))
+    test_output.append(np.array([0, 2, 4]))
+    # test 2
+    test_input.append(np.array([1]))
+    test_output.append(np.array([1, 3, 5]))
+    for input_rcu_index, expected_output in zip(test_input, test_output):
+        actual_output = rcu_index_2_gsi(input_rcu_index)
+        if not np.all(expected_output == actual_output):
+            print(expected_output)
+            print(actual_output)
+            raise Exception('Tests fails!')
+
+
+def test_rcu_index_and_rcu_input_index_2_gsi():
+    test_input = []
+    test_output = []
+    # test 1
+    test_input.append([np.array([0]), np.array([0])])
+    test_output.append(np.array([0]))
+    # test 2
+    test_input.append([np.array([0]), np.array([0, 1, 2])])
+    test_output.append(np.array([0, 2, 4]))
+    for inputs, expected_output in zip(test_input, test_output):
+        input_rcu_index = inputs[0]
+        input_rcu_input_index = inputs[1]
+        actual_output = rcu_index_and_rcu_input_index_2_gsi(input_rcu_index,
+                                                            input_rcu_input_index)
+        if not np.all(expected_output == actual_output):
+            print(expected_output)
+            print(actual_output)
+            raise Exception('Tests fails!')
+
+
+def test_gsi_2_rcu_index_and_rcu_input_index():
+    test_input = []
+    test_output = []
+    # test 1
+    test_input.append(np.array([0]))
+    test_output.append((np.array([0]), np.array([0])))
+    # test 2
+    test_input.append(np.array([1]))
+    test_output.append((np.array([1]), np.array([0])))
+    for input_signal_index, expected_output in zip(test_input, test_output):
+        actual_output = gsi_2_rcu_index_and_rcu_input_index(input_signal_index)
+        print(actual_output)
+        print(expected_output)
+        if not np.all(expected_output == actual_output):
+            print(expected_output)
+            print(actual_output)
+            raise Exception('Tests fails!')
+
+
+def get_timestamp(f=None):
+    '''
+    Get the timestamp in standard format
+
+    Input argument:
+    f = format specifier.
+        iso format (default)
+        'filename': filename without spaces and special characters,
+                    format: yyyymmddThhmmss
+
+    Output argument:
+    timestamp = timestamp in correct format
+    '''
+    timestamp = datetime.datetime.isoformat(datetime.datetime.now())
+    if f == 'filename':
+        timestamp = datetime.datetime.now().strftime("%Y%m%dT%H%M%S")
+    return timestamp
+
+
+def log(string):
+    print(f'{get_timestamp()} : {string}')
+
+
+def get_valid_band_names():
+    band_names = []
+    for band_name in bandlims:
+        band_names.append(band_name)
+    return band_names
+
+
+def get_band_name_and_subband_index_for_frequency(freqs):
+    '''
+    Get band name and subband index for one or more frequencies.
+
+    Input arguments:
+    freqs = one or more frequencies in Hz
+    bandlims = dict with band names and their limits (default: LOFAR2 bands)
+
+    Output arguments:
+    band_names = list of band names (see get_valid_band_names)
+    sbi        = one or more indices of the subbands for which the frequencies
+                 should be returned
+    '''
+    # make sure freqs is a list
+    if not isinstance(freqs, list):
+        freqs = [freqs]
+    # get band name and subband index
+    freqs = np.array(freqs)
+    band_names = []
+    for freq in freqs:
+        band_name = get_band_name_for_frequency(freq)
+        band_names.append(band_name)
+    # get indices for band_name 'LB'
+    sbis = np.round(freqs*cst_n_sub/cst_fs)
+    # and update for the other bands:
+    for idx in np.where(np.array(band_names) == 'HB1')[0]:
+        sbis[idx] = 2*cst_n_sub - sbis[idx]
+    for idx in np.where(np.array(band_names) == 'HB2')[0]:
+        sbis[idx] = sbis[idx] - 2*cst_n_sub
+    return band_names, sbis
+
+
+def get_band_name_for_frequency(freq, bandlims=bandlims):
+    '''
+    Get band name for frequency
+
+    Input arguments:
+    freq = one frequency in Hz, float
+    bandlims = dict with band names and their limits (default: LOFAR2 bands)
+
+    Output arguments:
+    band_name = name of analog band (see get_valid_band_names),
+                empty string ('') if unsuccessful
+    '''
+    for band_name in bandlims:
+        if np.min(bandlims[band_name]) <= freq and freq <= np.max(bandlims[band_name]):
+            return band_name
+    return ''
+
+
+def get_frequency_for_band_name(band_name):
+    '''
+    Get all frequencies for given band name
+
+    Input arguments:
+    band_name = name of analog band (see get_valid_band_names)
+
+    Output arguments:
+    freqs = frequencies in Hz for the selected band
+    '''
+    return get_frequency_for_band_name_and_subband_index(band_name, sbi=np.arange(512))
+
+
+def get_frequency_for_band_name_and_subband_index(band_name, sbi):
+    '''
+    Get frequency for band name and subband(s)
+
+    Input arguments:
+    band_name = name of analog band (see get_valid_band_names)
+    sbi       = one or more indices of the subbands for which the frequencies
+                should be returned
+
+    Output arguments:
+    freqs = frequencies in Hz for the selected band and subband
+    '''
+    # check input
+    valid_band_names = get_valid_band_names()
+    if band_name not in valid_band_names and band_name not in LB_aliases:
+        raise Exception('Unknown band_name "%s". Valid band_names are: %s' %
+                        (band_name, valid_band_names))
+    # generate frequencies
+    freqs = sbi*cst_fs/cst_n_sub
+    if band_name == 'HB1':
+        return 200e6 - freqs
+    elif band_name == 'HB2':
+        return 200e6 + freqs
+    # else: # LB
+    return freqs
+
+
+def get_band_names(RCU_band_select_R, RCU_PCB_version_R,
+                   LB_tags=['RCU2L'], HB_tags=['RCU2H']):
+    '''
+    Get the frequency band names per receiver, based on their band selection
+    and RCU PCB name.
+
+    PCB name is used to determine the band (Low Band 'LB' or High Band 'HB',
+    None otherwise). The band selection is added as a postfix.
+    RCU PCB version can also be used, but then the default tags will not be sufficient
+
+    Input Arguments:
+    RCU_band_select_R = as returned by recv.RCU_band_select_R
+    RCU_PCB_version_R = as returned by recv.RCU_PCB_version_R
+    LB_tags           = substrings in RCU_PCB_version_R fields to indicate Low Band
+    HB_tags           = substrings in RCU_PCB_version_R fields to indicate High Band
+
+    Note:
+    RCU_PCB_version_R can also hold IDs as returned by recv.RCU_PCB_ID_R
+    In that case, the LB_tags and HB_tags should be passed on by the user
+
+    Output Arguments:
+    band_name = list of strings indicating the band name per receiver.
+                containts None if no band name could be determined
+
+    '''
+    # if IDs are passed on, convert to list of strings
+    if isinstance(RCU_PCB_version_R, np.ndarray):
+        RCU_PCB_version_R = [str(x) for x in RCU_PCB_version_R]
+    #
+    signal_indices = np.arange(RCU_band_select_R.shape[0] * RCU_band_select_R.shape[1])
+    band_name_per_signal_index = [None]*len(signal_indices)
+
+    for rcu_index in range(len(RCU_band_select_R)):
+        for rcu_input_index in range(len(RCU_band_select_R[rcu_index])):
+            signal_index = rcu_index_and_rcu_input_index_2_gsi(rcu_index,
+                                                               rcu_input_index)
+            this_band = None
+            for tag in LB_tags:
+                if tag in RCU_PCB_version_R[rcu_index]:
+                    this_band = 'LB%1s' % RCU_band_select_R[rcu_index][rcu_input_index]
+                    break
+            if this_band is None:  # continue with high band
+                for tag in HB_tags:
+                    if tag in RCU_PCB_version_R[rcu_index]:
+                        this_band = 'HB%1s' % (
+                            RCU_band_select_R[rcu_index][rcu_input_index])
+                        break
+            band_name_per_signal_index[signal_index] = this_band
+    return band_name_per_signal_index
+
+
+# some basic test methods, for internal use:
+def test_get_band_name_and_subband_index_for_frequency():
+    input_frequencies = [138.1e6, 137.9e6, 0.1, 0, 195312.5, 234e6]
+    expected_band_names = ['HB1', 'HB1', 'LB', 'LB', 'LB', 'HB2']
+    expected_subband_index = np.array([317., 318., 0., 0., 1., 174.])
+    actual_band_names, actual_subband_index = \
+        get_band_name_and_subband_index_for_frequency(input_frequencies)
+    if not expected_band_names == actual_band_names:
+        print(expected_band_names)
+        print(actual_band_names)
+        raise Exception('Tests fails!')
+    if not (expected_subband_index == actual_subband_index).all():
+        print(expected_subband_index)
+        print(actual_subband_index)
+        raise Exception('Tests fails!')
+
+
+def test_get_frequency_for_band_name_and_subband_index():
+    for input_band_name, input_subband_index, expected_frequencies in zip(
+      ['HB1', 'HB1', 'LB', 'LB', 'LB', 'HB2', 'LB1'],
+      [195., 194., 0., 0., 1., 174., 1.],
+      [161914062.5, 162109375.0, 0., 0., 195312.5, 233984375, 195312.5]):
+        actual_frequencies = get_frequency_for_band_name_and_subband_index(
+            input_band_name, input_subband_index)
+        if not expected_frequencies == actual_frequencies:
+            print(expected_frequencies)
+            print(actual_frequencies)
+            raise Exception('Tests fails!')
+
+
+def test():
+    test_rcu_index_2_gsi()
+    test_rcu_index_and_rcu_input_index_2_gsi()
+    test_gsi_2_rcu_index_and_rcu_input_index()
+    test_get_band_name_and_subband_index_for_frequency()
+    test_get_frequency_for_band_name_and_subband_index()
+
+
+def main():
+    return False
+
+
+if __name__ == '__main__':
+    exit(main())