Skip to content
Snippets Groups Projects
Select Git revision
  • a795a90d8c45cd71b80cf7d88e53a95a5bf61c58
  • main default protected
  • 65_async_query
  • 169_codemeta
  • issue/152_cleaning_up
  • adex-main protected
  • sdc380-aladin-cone-search
  • sdc-222_multi_archive_query_review
  • 69_add_diracIAM
  • dev-nico
  • dev-dirac
  • acceptance
  • dev-zooniverse
13 results

ZooniverseResults.js

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    apertif_arts_firmware_model.py 69.88 KiB
    ###############################################################################
    #
    # Copyright (C) 2018
    # ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
    # P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
    #
    # This program is free software: you can redistribute it and/or modify
    # it under the terms of the GNU General Public License as published by
    # the Free Software Foundation, either version 3 of the License, or
    # (at your option) any later version.
    #
    # This program is distributed in the hope that it will be useful,
    # but WITHOUT ANY WARRANTY; without even the implied warranty of
    # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    # GNU General Public License for more details.
    #
    # You should have received a copy of the GNU General Public License
    # along with this program.  If not, see <http://www.gnu.org/licenses/>.
    #
    ###############################################################################
    
    # Author: Eric Kooistra, 28 Oct 2018,
    
    """Model the Apertif BF, XC / Arts TAB, IAB firmware data path signal levels
    
    The Apertif / Arts data path processing is described in:
    
      ASTRON_SP_062_Detailed_Design_Arts_FPGA_Beamformer.pdf
    
    The purpose of the apertif_arts_firmware_model is to understand signal levels
    in the data path processing and to prove that the internal requantization in
    the firmware implementation is does not increase the system noise while
    providing the maximum dynamic range, given data widths at the external
    interfaces (i.e. W_adc = 8, W_beamlet = 6, W_vis = 32, W_volt = 8 and
    W_pow = 8).
     
    This apertif_arts_firmware_model.py model consists of the following sections:
    
    - User parameters
    - Firmware parameters
    - A static model to determine the expected internal signal levels as function
      of the WG amplitude or the ADC noise sigma level.
    - A dynamic model using FFT, BF and ()**2 to model the data path signal
      processing for a WG input signal or a system noise input signal. The FIR
      part of the F_sub and F_chan is not modelled.
      . WG input for test purposes to verify the linearity and dynamic range
        as function of the WG amplitude and the number of dishes.
      . Noise input for operational purposes to verify that the system noise 
        maintains represented by suffcient LSbits throughout the data path
        processing chain from ADC via beamlet and fine channel to correlator
        visibility, Arts SC beamlet TAB, Arts SC4 power TAB and IAB.
    - Plot results of the static and the dynamic model.
    
    The model matches the drawings of the firmware parameters and signal levels in:
    
      apertif_arts_firmware_model.vsd
      apertif_arts_firmware_model.pdf 
    
    More information and a log an plots are also reported in:
    
      ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf
      
    Parameters:
    * wgEnable
      When True then a WG input is used for the dynamic model, else a noise input.
    * wgSweep
      The WG makes a sweep through wgSweep number of subbands if wgSweep != 0. The
      wgSweep can be < 1 and negative. If wgSweep != 0 then an image plot is made
      of the subbands spectrum and the channels spectrum. The expected signal
      levels that are reported by the dynamic model are only relevant when wgSweep
      = 0.
    * quantize
      The quantize parameter provides rounding of internal signal levels, weights
      and gain factors. For the dynamic model clipping is not supported, because
      the static model already shows when the signals clip. In the static model the
      clipping can be controlled via the adjust parameter.
    
    Remark:
    - power = rms**2
    - mean powers = (rms voltages)**2
    - rms**2 = std**2 + mean**2 = std**2 when mean = 0
    - rms = std = ampl / sqrt(2) when mean is 0
      --> pi_apr.wgScale = sqrt(2)
      
    - power complex = power real + power imag
                    = (rms real)**2 + (rms imag)**2
                    = (ampl real/sqrt(2))**2 + (ampl imag/sqrt(2))**2
    - power real = power imag = power complex / 2, when signal is random or periodic
    - rms real = rms imag = rms complex / sqrt(2), when signal is random or periodic
      --> pi_apr.rcScale = sqrt(2) = sqrt(nof_complex)
    - rms complex = sqrt(power complex)
    - ampl real = ampl imag = sqrt(power complex / 2) * sqrt(2)
                            = sqrt(power complex)
                            = rms complex
                            = (rms real) * sqrt(2)
                            = (rms imag) * sqrt(2)
    
    A = 1.0
    N = 100000
    realSamples = A * np.random.randn(N)
    realPowers = realSamples * realSamples
    x = rms(realSamples)**2
    y = np.mean(realPowers)
    z = np.std(realPowers)
    x # = A * 1.0
    y # = A * 1.0
    z # = A * sqrt(2)
    complexSamples =  A * np.random.randn(N) +  A * np.random.randn(N) * 1j
    complexPowers = complexSamples * np.conjugate(complexSamples)
    complexPowers = complexPowers.real
    x  = rms(complexSamples)**2
    xr = rms(complexSamples.real)**2
    xi = rms(complexSamples.imag)**2
    y = np.mean(complexPowers)
    z = np.std(complexPowers)
    x # = A * 2.0
    xr
    xi
    y # = A * 2.0
    z # = A * 2.0
    
    Use prefix rmsComplexY for rms of complex signal Y
    Use prefix rmsY for rms of real signal Y or of Y.real or of Y.imag, so
        for complex signal Y then rmsY = rmsComplexY / pi_apr.rcScale,
        for real signal Y then rmsY = rmsY.
    
    Usage:
    # in iPython:
    import os
    filepath='E:\\svn_radiohdl\\trunk\\applications\\apertif\\matlab'
    os.chdir(filepath)
    run apertif_arts_firmware_model.py
    
    # on command line
    python apertif_arts_firmware_model.py
    
    """
    
    import argparse
    import math
    import numpy as np
    import common as cm
    import test_plot
    import pi_apertif_system as pi_apr
    
    ################################################################################    
    # General
    ################################################################################    
    
    def rms(arr):
        """Root mean square of values in arr
        Use this root-mean-square rms() instead of std() for CW signals that are at
        center of bin and result in static phasor with std() = 0. For noise
        signals with zero mean rms() is equal to std(). This rms() also works for
        complex input thanks to using np.abs().
        """
        return np.sqrt(np.mean(np.abs(arr)**2.0))
    
    nofSigma = 4    # Number of sigma to fit in full scale
    
    sc1_use_qua = True       # use SC1 TAB bf quantized and clipped output
    sc1_use_qua = False      # use SC1 TAB bf raw LSbits and wrapped MSbits output (as implemented in arts_unb1_sc1)
    
    ################################################################################    
    # Parse arguments to derive user settings
    ################################################################################    
    
    _parser = argparse.ArgumentParser('apertif_arts_firmware_model')
    _parser.add_argument('--app', default='all', type=str, help='Application to model: dish, apertif, arts_sc1, arts_sc4 or all')
    _parser.add_argument('--model', default='all', type=str, help='Model: static, dynamic or all')
    _parser.add_argument('--sky_sigma', default=4.0, type=float, help='Signal sigma level at ADC input')
    _parser.add_argument('--wg_enable', default=False, action='store_true', help='Model coherent ADC input using WG sinus or model incoherent ADC input using Gaussian white noise')
    _parser.add_argument('--wg_ampl', default=4.0 * 2**0.5, type=float, help='WG amplitude at ADC input')
    _parser.add_argument('--wg_sweep', default=0.0, type=float, help='Number of subbands fsub to sweep with WG in dynamic model')
    _parser.add_argument('--nof_periods', default=250, type=int, help='Number of channel periods in dynamic model')
    _parser.add_argument('--nof_dishes', default=1, type=int, help='Number of dishes')
    _parser.add_argument('--quantize', default=False, action='store_true', help='Use fixed point rounding or no rounding of weights, gains and internal signals')
    _parser.add_argument('--useplot', default=False, action='store_true', help='Use plots or skip plotting')
    _parser.add_argument('--cb_weight_adjust', default=1.0, type=float, help='Factor to adjust unit CB weight')
    _parser.add_argument('--tab_weight_adjust_sc1', default=1.0, type=float, help='Factor to adjust unit TAB weight for SC1')
    _parser.add_argument('--tab_weight_adjust_sc4', default=1.0, type=float, help='Factor to adjust unit TAB weight for SC4')
    _parser.add_argument('--iab_gain_adjust', default=1.0, type=float, help='Factor to adjust unit IAB gain for SC4')
    args = _parser.parse_args()
    
    if args.app=='all':
        apps = ['dish', 'apertif', 'arts_sc1', 'arts_sc4']
    else:
        apps = [args.app]
    if args.model=='all':
        models = ['static', 'dynamic']
    else:
        models = [args.model]
    
    skySigma  = args.sky_sigma
    #skySigma  = 4
    #skySigma  = 89.8        # 89.8 * 2**0.5 = 127 for full scale WG amplitude
    #skySigma  = 26.1        # 26.1 * 2**0.5 =  37
    #skySigma  = 11.3        # 11.3 * 2**0.5 =  16 for full scale 8 bit beamlet amplitude
    
    wgAmpl    = args.wg_ampl
    wgSigma   = wgAmpl / np.sqrt(2)
    
    wgEnable  = args.wg_enable
    wgSweep   = args.wg_sweep
    #wgSweep   = 5.0          # number of fsub to sweep over nofTchanx
    #wgSweep   = -5.0         # number of fsub to sweep over nofTchanx
    #wgSweep   = 2.0          # number of fsub to sweep over nofTchanx
    #wgSweep   = 1.0/32       # number of fsub to sweep over nofTchanx
    if wgSweep != 0:
        wgEnable = True   # force WG
        apps = ['dish']   # only support sweep for dish
    
    quantize  = args.quantize
    
    useplot   = args.useplot
    
    cbWeightAdjust     = args.cb_weight_adjust
    tabWeightAdjustSc1 = args.tab_weight_adjust_sc1
    tabWeightAdjustSc4 = args.tab_weight_adjust_sc4
    iabGainAdjust      = args.iab_gain_adjust
    #cbWeightAdjust     = 1.0/8    # use <= 1/8 to have no beamlet (8 bit) clipping
    #cbWeightAdjust     = 1.0/32   # use <= 1/32 to have no fine channel_x (9 bit) clipping
    #tabWeightAdjustSc1 = 1.0      # = 1.0 to use unit weight
    #tabWeightAdjustSc4 = 1.0/4    # use <= 1/4 and cbWeightAdjust = 1/32 to have no TAB Stokes I (8 bit) clipping
    #iabGainAdjust      = 1.0/16   # use <= 1/16 and cbWeightAdjust = 1/32 to have no IAB Stokes I (8 bit) clipping
    
    # Number of dishes
    # - not used for Apertif correlator, because it only models auto visibilities
    # - only used for Arts array beamformers:
    #   . use  1 to verify the expected TAB, IAB output signal level
    #   . use 12 is all dishes for Arts SC1 or for IAB
    #   . use 10 is equidistant dishes for Arts SC4
    #   . use  8 is equidistant and fixed dishes for Arts SC3
    # The actual number of dishes does not alter the TAB or IAB output level, 
    # because the number of dishes is compensated via the TAB weights in case of
    # SC1 when sc1_use_qua = True and SC4 or via the IAB gain in case of the IAB.
    # When sc1_use_qua = False then the TAB output does increase with the number
    # of dishes, because then the TAB weigth is fixed at 1. Still for SC1 the
    # maximum gain with N_dish = 12 dishes is log2(sqrt(12)) = 1.8 bits, so with
    # 6 bit beamlet input and 8 bit voltage TAB data this gain still fits, and
    # therefore it is fine to use sc1_use_qua = False for SC1.
    # Typical values:
    # . nofDishes = 16    # model only
    # . nofDishes = 12    # RT2-9,a-d, all for Apertif XC and Arts SC1
    # . nofDishes = 10    # RT2-9,a,b at equidistant for Arts SC4
    # . nofDishes = 8     # RT2-9 fixed at equidistant for commensal Arts SC3
    # . nofDishes = 1     # model default, because tabWeightSc1, tabWeightSc4 and
    #                     # iabGainReg nomalize for nofDishes. Therefore use 
    #                     # nofDishes = 1 to speed up simulation, except when modeling
    #                     # Arts SC1 with fixed tabWeightSc1 = 1 (when sc1_use_qua =
    #                     # False).
    nofDishes = args.nof_dishes
    
    # Time series
    if 'dynamic' in models:
        # use args.nof_periods for number of channel periods to simulate:
        if wgEnable:
            nofTchanx = args.nof_periods
            #nofTchanx = 10          # use >~10 for WG
            #nofTchanx = 100
        else:
            nofTchanx = args.nof_periods
            #nofTchanx = 100         # number of channel periods to simulate, use >~ 100 for noise
            #nofTchanx = 1000        # use >= 1000 for more accuracy
            #nofTchanx = 250        # use < 1000 to speed up simulation
        nofTStokes = nofTchanx
    
        # also use args.nof_periods for number of Ts, Tsub, Tchanx periods to plot
        # . to speed up plotting for Ts and Tsub
        # . to have same density of samples per plot
        nof_periods = args.nof_periods
    
    # Frequency selection
    if 'dynamic' in models:
        # Select frequency for WG, SST, BST, XC
        selSub           = 65
        selChanx         = 33          # N_chan_x / 2 = center of subband
        selChanxFraction = 1 * 1.0/nofTchanx   # yields 1 rotation in fine channel
        #selChanxFraction = 0.0         # center channel, yields constant in fine channel
        if wgSweep:
            selChanx         = 32      # N_chan_x / 2 = center of subband
            selChanxFraction = 0.0     # center channel, yields constant in fine channel
    
    ################################################################################
    # Derived gains
    
    # Aperif BF subband CB voltage weight
    cbWeight = pi_apr.cbUnitWeight * cbWeightAdjust         # = 8192 = 2**13 * 1.0
    if quantize:
        cbWeight = np.round(cbWeight)
    cbGain = cbWeight/pi_apr.cbMaxWeight   # normalized gain relative to maximum gain 1.0
    
    # Arts SC1 beamlet TAB voltage weight
    if sc1_use_qua:
        # Intended implementation using quantized and clipped BF data with 
        # TAB weights that are scaled by 1/nofDishes for coherent WG input
        # and by 1/sqrt(nofDishes) for incoherent noise input.
        W_tab_unit_scale_sc1 = -2                                   # = -2
        W_tab_weight_fraction_sc1 = pi_apr.W_tab_max_weight_sc1 + W_tab_unit_scale_sc1  # = 
        tabUnitWeightSc1 = 2.0**W_tab_weight_fraction_sc1           # = 2**13 = 8192
        if wgEnable:
            tabWeightSc1 = tabUnitWeightSc1 / nofDishes             # coherent TAB input, like WG
        else:
            tabWeightSc1 = tabUnitWeightSc1 / math.sqrt(nofDishes)  # incoherent TAB input, like sky noise
        tabWeightSc1 *= tabWeightAdjustSc1
    else:
        # Current implementation using LSbits and wrapped MSbits of raw BF data
        # and unit weight = 1 independent of nofDishes. The TAB weight
        # tabWeightSc1 = 1 is the smallest weight value, so effectively the range
        # of W_tab_weight_sc1 is not used.
        W_tab_unit_scale_sc1 = -pi_apr.W_tab_max_weight_sc1         # = -15
        W_tab_weight_fraction_sc1 = pi_apr.W_tab_max_weight_sc1 + W_tab_unit_scale_sc1 # = 0 = 15 + -15
        tabUnitWeightSc1 = 1.0                                      # = 1 = 2.0**W_tab_weight_fraction_sc1 = 2**0
        tabWeightSc1 = tabUnitWeightSc1                             # = 1 independent of nofDishes
    if quantize:
        tabWeightSc1 = np.round(tabWeightSc1)
    tabGainSc1 = tabWeightSc1/pi_apr.tabMaxWeightSc1  # Normalized gain relative to maximum gain 1.0
    
    # Arts SC4 channel TAB voltage weight
    if wgEnable:
        tabWeightSc4 = pi_apr.tabUnitWeightSc4 / nofDishes             # coherent TAB input, like WG
    else:
        tabWeightSc4 = pi_apr.tabUnitWeightSc4 / math.sqrt(nofDishes)  # incoherent TAB input, like sky noise
    tabWeightSc4 = tabWeightSc4 * tabWeightAdjustSc4
    if quantize:
        tabWeightSc4 = np.round(tabWeightSc4)
    tabGainSc4 = tabWeightSc4/pi_apr.tabMaxWeightSc4    # Normalized gain relative to maximum gain 1.0
    
    # Arts SC4 power IAB power gain
    iabGainReg = pi_apr.iabUnitGain * iabGainAdjust / nofDishes     # IAB input powers add coherently
    if quantize:
        iabGainReg = np.round(iabGainReg)
    iabGain = iabGainReg/pi_apr.iabMaxGain    # Normalized gain relative to maximum gain 1.0
    
    
    ################################################################################
    # Apertif / Arts settings
    # . BF
    N_fft         = pi_apr.N_FFT            # = 1024
    N_sub         = pi_apr.N_sub            # = 512
    N_int_sub     = pi_apr.N_int_x          # = 800000
    
    # . XC, SC4
    N_chan_x      = pi_apr.N_chan_x         # = 64
    N_chan_bf     = pi_apr.N_chan_bf        # = 4
    N_int_chan_x  = pi_apr.N_int_chan_x     # = 12500 = N_int_x / N_chan_x
    
    if 'static' in models:
        ################################################################################    
        # Static model of signal level as function of:
        # . WG input amplitude level
        # . Noise input sigma level
        ################################################################################    
        
        # Input axis
        inputResolution = 0.05
        
        # Dish
        A_wg = np.arange(inputResolution, 2**pi_apr.W_adc_max, inputResolution)    # 0:127
        S_wg = A_wg / pi_apr.wgScale
        S_wg_complex = S_wg
        S_wg_complex_phase = S_wg_complex / pi_apr.rcScale
        A_wg_complex_phase = S_wg_complex
        S_noise = np.arange(inputResolution, 2**pi_apr.W_adc_max / nofSigma, inputResolution)  # 0:31
        S_noise_complex = S_noise
        S_noise_complex_phase = S_noise / pi_apr.rcScale
        
        SST_wg_inputs = S_wg_complex / pi_apr.wgSstScale
        SST_wg_values = SST_wg_inputs**2 * N_int_sub
        SST_wg_dBs = 10*np.log10(SST_wg_values)
        SST_noise_inputs = S_noise_complex / pi_apr.noiseSstScale
        SST_noise_values = SST_noise_inputs**2 * N_int_sub
        SST_noise_dBs = 10*np.log10(SST_noise_values)
        
        BST_wg_inputs = S_wg_complex * cbWeightAdjust / pi_apr.wgBstScale
        BST_wg_values = BST_wg_inputs**2 * N_int_sub
        BST_wg_dBs = 10*np.log10(BST_wg_values)
        BST_noise_inputs = S_noise_complex * cbWeightAdjust / pi_apr.noiseBstScale
        BST_noise_values = BST_noise_inputs**2 * N_int_sub
        BST_noise_dBs = 10*np.log10(BST_noise_values)
        
        clip = 2.0**pi_apr.W_beamlet_max - 1
        beamlet_wg_complex_phase_values = S_wg_complex_phase * cbWeightAdjust / pi_apr.wgBeamletScale
        beamlet_wg_complex_phase_values = np.clip(beamlet_wg_complex_phase_values, -clip, clip)
        beamlet_noise_complex_phase_values = S_noise_complex_phase * cbWeightAdjust / pi_apr.noiseBeamletScale
        beamlet_noise_complex_phase_values = np.clip(beamlet_noise_complex_phase_values, -clip, clip)
        
        # Central
        clip = 2.0**pi_apr.W_channel_x_max - 1
        channel_x_wg_complex_phase_values = S_wg_complex_phase * cbWeightAdjust / pi_apr.wgChannelXcorScale
        channel_x_wg_complex_phase_values = np.clip(channel_x_wg_complex_phase_values, -clip, clip)
        channel_x_wg_complex_values = channel_x_wg_complex_phase_values * pi_apr.rcScale
        channel_x_noise_complex_phase_values = S_noise_complex_phase * cbWeightAdjust / pi_apr.noiseChannelXcorScale
        channel_x_noise_complex_phase_values = np.clip(channel_x_noise_complex_phase_values, -clip, clip)
        channel_x_noise_complex_values = channel_x_noise_complex_phase_values * pi_apr.rcScale
        
        if 'apertif' in apps:
            # Correlator
            XC_wg_auto_visibilities = channel_x_wg_complex_values**2 * N_int_chan_x
            XC_wg_auto_visibilities_dB = 10*np.log10(XC_wg_auto_visibilities)
            XC_noise_auto_visibilities = channel_x_noise_complex_values**2 * N_int_chan_x
            XC_noise_auto_visibilities_dB = 10*np.log10(XC_noise_auto_visibilities)
        
        if 'arts_sc1' in apps:
            # Arts SC1 beamlet voltage TAB
            if sc1_use_qua:
                beamlet_tab_wg_complex_phase_values_sc1 = beamlet_wg_complex_phase_values * tabWeightAdjustSc1
                beamlet_tab_noise_complex_phase_values_sc1 = beamlet_noise_complex_phase_values * tabWeightAdjustSc1
            else:
                beamlet_tab_wg_complex_phase_values_sc1 = beamlet_wg_complex_phase_values * nofDishes
                beamlet_tab_noise_complex_phase_values_sc1 = beamlet_noise_complex_phase_values * nofDishes
            clip = 2.0**pi_apr.W_volt_max - 1
            beamlet_tab_wg_complex_phase_values_sc1 = np.clip(beamlet_tab_wg_complex_phase_values_sc1, -clip, clip)
            beamlet_tab_noise_complex_phase_values_sc1 = np.clip(beamlet_tab_noise_complex_phase_values_sc1, -clip, clip)
        
        if 'arts_sc4' in apps:
            # Arts SC4 channel power TAB
            clip = 2.0**pi_apr.W_volt - 1
            power_tab_xx_wg_complex_values_sc4 = (channel_x_wg_complex_values * tabWeightAdjustSc4)**2    # Stokes I = XX* + 0
            power_tab_xx_wg_complex_values_sc4 = np.clip(power_tab_xx_wg_complex_values_sc4, 0, clip)  # Stokes I is unsigned >= 0
            power_tab_xx_noise_values_sc4 = (channel_x_noise_complex_values * tabWeightAdjustSc4)**2    # Stokes I = XX* + 0
            power_tab_xx_noise_values_sc4 = np.clip(power_tab_xx_noise_values_sc4, 0, clip)  # Stokes I is unsigned >= 0
            
            # Arts SC4 channel IAB
            clip = 2.0**pi_apr.W_pow - 1
            iab_xx_wg_complex_values_sc4 = channel_x_wg_complex_values**2 * iabGainAdjust   # Stokes I = XX* + 0
            iab_xx_wg_complex_values_sc4 = np.clip(iab_xx_wg_complex_values_sc4, 0, clip) # Stokes I is unsigned >= 0
            iab_xx_noise_values_sc4 = channel_x_noise_complex_values**2 * iabGainAdjust   # Stokes I = XX* + 0
            iab_xx_noise_values_sc4 = np.clip(iab_xx_noise_values_sc4, 0, clip) # Stokes I is unsigned >= 0
    
    
    if 'dynamic' in models:
        ################################################################################    
        # Dynamic model parameters
        ################################################################################    
        
        ################################################################################
        # Apertif / Arts settings
        # . BF
        fs            = 800e6
        fsub          = fs / N_fft
        Ts            = 1/fs
        Tsub          = 1/fsub
    
        # . XC, SC4
        fchan_x       = fsub / N_chan_x
        Tchanx        = 1/fchan_x
        
        # . SC4
        Tstokesx      = Tchanx
        M_int_ax      = pi_apr.M_int_ax         # = 16
        W_int_ax      = pi_apr.W_int_ax         # = 4 = log2(16)
        #M_int_ax      = 1
        #W_int_ax      = cm.ceil_log2(M_int_ax)
        
        ################################################################################
        # Dish CB
    
        # Number of bits for internal product
        W_cb_product = pi_apr.W_cb_in + pi_apr.W_cb_weight - 1  # = 31 = 16 + 16 -1, do -1 to skip double sign bit of product
        
        # Number of bits for fraction < 1
        W_beamlet_fraction = pi_apr.W_beamlet_fraction    # number of bits for fraction < 1
        
        if 'arts_sc1' in apps:
            ################################################################################    
            # Arts SC1 beamlet TAB
            
            # Number of bits for internal product
            W_tab_product_sc1 = pi_apr.W_beamlet + pi_apr.W_tab_weight_sc1 - 1  # = 23 = 8 + 16 -1, do -1 to skip double sign bit of product
            
            # Number of bits for fraction < 1
            W_tab_fraction_sc1 = W_beamlet_fraction - W_tab_unit_scale_sc1  # = 8  = 6 - -2  when sc1_use_qua = True
                                                                            # = 21 = 6 - -15 when sc1_use_qua = False
        
        if 'arts_sc4' in apps:
            ################################################################################    
            # Arts SC4 channel TAB
            
            # Fine channel TAB        
            # Number of bits for internal product
            W_tab_voltage_product_sc4 = pi_apr.W_channel_x + pi_apr.W_tab_weight_sc4 - 1  # = 17 = 9 + 9 -1, do -1 to skip double sign bit of product
            W_tab_power_product_sc4 = 2*W_tab_voltage_product_sc4 - 1  # = 33 = 2*17 -1, do -1 to skip double sign bit of product
            
            # For SC4 it is possible to use W_channel_x_fraction > 9, because the F_chan is internal
            # For SC3 it W_channel_x_fraction = 9, because the F_chan is in Apertif XC
            # Number of bits for fraction < 1
            W_tab_voltage_fraction_sc4 = pi_apr.W_channel_x_fraction - pi_apr.W_tab_unit_scale_sc4  # = 10 = 9 - -1
            
            # Arts channel TAB
            # Number of bits for fraction < 1
            W_tab_power_fraction_sc4 = 2*W_tab_voltage_fraction_sc4             # = 20 = 2*10
            W_tab_power_int_fraction_sc4 = W_tab_power_fraction_sc4 - W_int_ax  # = 16 = 20 - 4
                                                                                                                     
            ################################################################################    
            # Arts SC4 channel IAB
                    
            # Number of bits for internal product
            W_iab_power_product = 2*pi_apr.W_channel_x_fraction - 1  # = 17 = 2*9 -1, do -1 to skip double sign bit of product
            W_iab_gain_product = W_iab_power_product + pi_apr.W_iab_gain - 1  # = 32 = 17 + 16 - 1, do -1 to skip double sign bit of product
            
            # Number of bits for fraction < 1
            W_iab_power_fraction = 2*pi_apr.W_channel_x_fraction              # = 18 = 2*9
            W_iab_gain_fraction = W_iab_power_fraction - W_int_ax - pi_apr.W_iab_unit_scale    # = 16 = 18 - 4 - -2
    
        ################################################################################    
        # Dynamic model time series
        ################################################################################    
        
        # Time, frequency axis
        nofTsub = nofTchanx * N_chan_x
        nofTs = nofTsub * N_fft
        t = np.arange(nofTs) * Ts
        
        ################################################################################
        # Dish Compound Beamformer
        
        # ADC input
        if wgEnable:
            wgFreq = (selSub + (selChanx + selChanxFraction - N_chan_x/2.0) / N_chan_x) * fsub
            chirpFreq = 0
            if wgSweep:
                chirpFreq = wgSweep * fsub * np.arange(nofTs)/nofTs/2.0
            spSamples = wgAmpl * np.sin(2*np.pi*(wgFreq + chirpFreq) * t)
        else:
            spSamples = np.random.randn(nofTs)
            spSamples *= skySigma / rms(spSamples)
        
        if quantize:
            # The unit level is at 1 lsb of the ADC input, so the W_adc = 8 bit input
            # has no fraction and can thus directly be rounded to model quantization.
            spSamples = np.round(spSamples)
        rmsSp = rms(spSamples)  # real signal
        amplSp = rmsSp * pi_apr.wgScale
        
        # Fsub subbands full bandwidth
        spBlocks = spSamples.reshape((nofTsub, N_fft))  # [nofTsub][N_fft]
        fftSamples = np.fft.fft(spBlocks) / N_fft       # [nofTsub][N_fft]
        fSubSamples = fftSamples[:, 0:N_sub]            # [nofTsub][N_sub], real input, so only keep positive frequencies
        
        # SST for selSub
        sstSamples = fSubSamples[:, selSub].copy()
        if quantize:
            sstSamples *= 2.0**pi_apr.W_sst_fraction  # fixed point fraction
            sstSamples = np.round(sstSamples)
            sstSamples /= 2.0**pi_apr.W_sst_fraction  # back to normalized fraction
        rmsComplexSST = rms(sstSamples)       # complex signal
        rmsSST = rmsComplexSST / pi_apr.rcScale
        amplSST = rmsSST * pi_apr.wgScale
        SSTpower = N_int_sub * rmsComplexSST**2
        
        # BF samples full bandwidth
        bfSamples = fSubSamples * cbGain
        
        # BST for selSub
        bstSamples = bfSamples[:, selSub].copy()
        if quantize:
            bstSamples *= 2.0**pi_apr.W_bst_fraction  # fixed point fraction
            bstSamples = np.round(bstSamples)
            bstSamples /= 2.0**pi_apr.W_bst_fraction  # back to normalized fraction
        rmsComplexBST = rms(bstSamples)    # complex signal
        rmsBST = rmsComplexBST / pi_apr.rcScale
        amplBST = rmsBST * pi_apr.wgScale
        BSTpower = N_int_sub * rmsComplexBST**2
        
        # CB beamlet for selSub
        beamletSamples = bfSamples[:, selSub].copy()
        if quantize:
            beamletSamples *= 2.0**pi_apr.W_beamlet_fraction  # fixed point fraction
            beamletSamples = np.round(beamletSamples)
            beamletSamples /= 2.0**pi_apr.W_beamlet_fraction  # back to normalized fraction
        rmsComplexBeamlet = rms(beamletSamples)    # complex signal
        rmsBeamlet = rmsComplexBeamlet / pi_apr.rcScale
        amplBeamlet = rmsBeamlet * pi_apr.wgScale
        
        ################################################################################
        # Beamlet channelizer
        # . Fchan output
        transposeSamples = bfSamples.transpose()                                # [nofTsub][N_sub] --> [N_sub][nofTsub]
        transposeBlocks = transposeSamples.reshape(N_sub, nofTchanx, N_chan_x)  # [N_sub][nofTchanx][N_chan_x]
        fSubChanSamples = np.fft.fft(transposeBlocks) / N_chan_x                # [N_sub][nofTchanx][N_chan_x], complex input and output
        fSubChanSamples = np.fft.fftshift(fSubChanSamples, 2)
        fSubChanSamples = fSubChanSamples.transpose(1, 0, 2)                    # [nofTchanx][N_sub][N_chan_x]
        
        # . CB fine channel for selSub, selChanx
        fineChannelSamples = fSubChanSamples[:, selSub, selChanx].copy()
        if quantize:
            fineChannelSamples *= 2.0**pi_apr.W_channel_x_fraction  # fixed point fraction
            fineChannelSamples = np.round(fineChannelSamples)
            fineChannelSamples /= 2.0**pi_apr.W_channel_x_fraction  # back to normalized fraction
        rmsComplexFineChannel = rms(fineChannelSamples)    # complex signal
        rmsFineChannel = rmsComplexFineChannel / pi_apr.rcScale
        amplFineChannel = rmsFineChannel * pi_apr.wgScale
        
        if 'apertif' in apps:
            ################################################################################
            # Aperif XC channel auto-visibility
            XCpower = N_int_chan_x * rmsComplexFineChannel**2
        
        if 'arts_sc1' in apps:
            ################################################################################
            # Arts SC1 beamlet voltage TAB
            if wgEnable:
                # Coherent sum using beamlets from original ADC input
                beamletTabSamples = nofDishes * beamletSamples.copy()
            else:
                # Incoherent sum using copies with the same rms beamlet as for the original ADC input
                beamletTabSamples = np.zeros(nofTsub) + np.zeros(nofTsub) * 1j
                for tel in range(nofDishes):
                    complexSamples = np.random.randn(nofTsub) + np.random.randn(nofTsub) * 1j
                    complexSamples *= rmsComplexBeamlet / rms(complexSamples)
                    beamletTabSamples += complexSamples
            beamletTabSamples *= tabGainSc1
            if quantize:
                beamletTabSamples *= 2.0**W_tab_fraction_sc1  # fixed point fraction
                beamletTabSamples = np.round(beamletTabSamples)
                beamletTabSamples /= 2.0**W_tab_fraction_sc1  # back to normalized fraction
            rmsComplexBeamletTab = rms(beamletTabSamples)    # complex signal
            rmsBeamletTab = rmsComplexBeamletTab / pi_apr.rcScale
            amplBeamletTab = rmsBeamletTab * pi_apr.wgScale
        
        if 'arts_sc4' in apps:
            ################################################################################
            # Arts SC4 Power TAB for single polarization Stokes I = XX* + 0
            # . Fine channel voltage TAB
            if wgEnable:
                # Coherent sum using channel from original ADC input
                fineChannelTabSamples = nofDishes * fineChannelSamples.copy()
            else:
                # Incoherent sum using copies with the same rms channel as for the original ADC input
                fineChannelTabSamples = np.zeros(nofTchanx) + np.zeros(nofTchanx) * 1j
                for tel in range(nofDishes):
                    complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
                    complexSamples *= rmsComplexFineChannel / rms(complexSamples)
                    fineChannelTabSamples += complexSamples    # voltage sum
            fineChannelTabSamples *= tabGainSc4
            if quantize:
                fineChannelTabSamples *= 2.0**W_tab_voltage_fraction_sc4  # fixed point fraction
                fineChannelTabSamples = np.round(fineChannelTabSamples)
                fineChannelTabSamples /= 2.0**W_tab_voltage_fraction_sc4  # back to normalized fraction
            rmsComplexFineChannelTab = rms(fineChannelTabSamples)    # complex signal
            rmsFineChannelTab = rmsComplexFineChannelTab / pi_apr.rcScale
            amplFineChannelTab = rmsFineChannelTab * pi_apr.wgScale
            
            # . Fine channel power TAB
            fineChannelTabPowers = fineChannelTabSamples * np.conjugate(fineChannelTabSamples)  # powers
            fineChannelTabPowers = fineChannelTabPowers.real           # keep real, imag = 0
            meanFineChannelTabPower = np.mean(fineChannelTabPowers)    # = rmsComplexFineChannelTab**2, use mean for power data and rms for voltage data
            stdFineChannelTabPower = np.std(fineChannelTabPowers)
            
            # . Integrated channels power TAB
            if wgEnable:
                # The WG sinus only occurs in one fine channel if it is close to the center fine channel frequency, so
                # the other fine channels in N_chan_bf fine channels carry almost zero and do not contribute much
                integratedTabPowers = fineChannelTabPowers.copy()
            else:
                # Integrate using copies with the same rmsComplexFineChannelTab as for the orginal ADC input
                integratedTabPowers = np.zeros(nofTchanx)
                for ci in range(M_int_ax):
                    complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
                    complexSamples *= rmsComplexFineChannelTab / rms(complexSamples)
                    powerSamples = complexSamples * np.conjugate(complexSamples)
                    powerSamples = powerSamples.real           # keep real, imag = 0
                    integratedTabPowers += powerSamples
            if quantize:
                integratedTabPowers *= 2.0**W_tab_power_int_fraction_sc4  # fixed point fraction
                integratedTabPowers = np.round(integratedTabPowers)
                integratedTabPowers /= 2.0**W_tab_power_int_fraction_sc4  # back to normalized fraction
            meanIntegratedTabPower = np.mean(integratedTabPowers)    # use mean for power data and rms for voltage data
            stdIntegratedTabPower = np.std(integratedTabPowers)
            
            ################################################################################
            # Arts SC4 IAB for single polarization Stokes I = XX* + 0
            
            # . Fine channel powers
            fineChannelPowers = fineChannelSamples * np.conjugate(fineChannelSamples)
            fineChannelPowers = fineChannelPowers.real           # keep real, imag = 0
            meanFineChannelPower = np.mean(fineChannelPowers)    # = rmsComplexFineChannel**2, use mean for power data and rms for voltage data
            
            # . Integrated channel powers
            if wgEnable:
                # The WG sinus only occurs in one fine channel if it is close to the center fine channel frequency, so
                # the other fine channels in N_chan_bf fine channels carry almost zero and do not contribute much
                integratedChannelPowers = fineChannelPowers
            else:
                # Integrate using copies with the same rmsComplexFineChannel as for the orginal ADC input
                integratedChannelPowers = np.zeros(nofTchanx)
                for ci in range(M_int_ax):
                    complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
                    complexSamples *= rmsComplexFineChannel / rms(complexSamples)
                    powerSamples = complexSamples * np.conjugate(complexSamples)
                    powerSamples = powerSamples.real           # keep real, imag = 0
                    integratedChannelPowers += powerSamples
            meanIntegratedChannelPower = np.mean(integratedChannelPowers)    # use mean for power data and rms for voltage data
            
            # . Sum of powers is coherent sum
            if wgEnable:
                # Sum of powers using integrated channel from original ADC input
                integratedIabPowers = nofDishes * integratedChannelPowers
            else:
                # Sum of powers using copies with the same rmsComplexFineChannel as for the original ADC input
                integratedIabPowers = np.zeros(nofTchanx)
                for tel in range(nofDishes):
                    telPowerSamples = np.zeros(nofTchanx)
                    for ci in range(M_int_ax):
                        complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
                        complexSamples *= rmsComplexFineChannel / rms(complexSamples)
                        powerSamples = complexSamples * np.conjugate(complexSamples)
                        powerSamples = powerSamples.real      # keep real, imag = 0
                        telPowerSamples += powerSamples       # powers sum integration
                    integratedIabPowers += telPowerSamples    # powers sum IAB
            meanIntegratedIabPower = np.mean(integratedIabPowers)    # use mean for power data and rms for voltage data
            
            # . Output gain
            outputIabPowers = integratedIabPowers * iabGain
            if quantize:
                outputIabPowers *= 2.0**W_iab_gain_fraction    # fixed point fraction
                outputIabPowers = np.round(outputIabPowers)
                outputIabPowers /= 2.0**W_iab_gain_fraction    # back to normalized fraction
            meanOutputIabPower = np.mean(outputIabPowers)    # use mean for power data and rms for voltage data
            stdOutputIabPower = np.std(outputIabPowers)
        
        ################################################################################
        # Dynamic model: Normalized floating point fraction to fixed point fraction
        ################################################################################    
        
        # SST
        rmsComplexSST *= 2.0**pi_apr.W_sst_fraction
        rmsSST *= 2.0**pi_apr.W_sst_fraction
        amplSST *= 2.0**pi_apr.W_sst_fraction
        SSTpower *= 2.0**(2*pi_apr.W_sst_fraction)
        
        # BST
        rmsComplexBST *= 2.0**pi_apr.W_bst_fraction
        rmsBST *= 2.0**pi_apr.W_bst_fraction
        amplBST *= 2.0**pi_apr.W_bst_fraction
        BSTpower *= 2.0**(2*pi_apr.W_bst_fraction)
        
        # CB beamlet for selSub
        beamletSamples *= 2.0**pi_apr.W_beamlet_fraction
        rmsComplexBeamlet *= 2.0**pi_apr.W_beamlet_fraction
        rmsBeamlet *= 2.0**pi_apr.W_beamlet_fraction
        amplBeamlet *= 2.0**pi_apr.W_beamlet_fraction
        
        # CB fine channel for selSub, selChanx
        fineChannelSamples *= 2.0**pi_apr.W_channel_x_fraction
        rmsComplexFineChannel *= 2.0**pi_apr.W_channel_x_fraction
        rmsFineChannel *= 2.0**pi_apr.W_channel_x_fraction
        amplFineChannel *= 2.0**pi_apr.W_channel_x_fraction
        
        if 'apertif' in apps:
            # Aperif XC channel auto-visibility
            XCpower *= 2.0**(2*pi_apr.W_channel_x_fraction)
        
        if 'arts_sc1' in apps:
            # Arts SC1 beamlet voltage TAB
            beamletTabSamples *= 2.0**W_tab_fraction_sc1
            rmsComplexBeamletTab *= 2.0**W_tab_fraction_sc1
            rmsBeamletTab *= 2.0**W_tab_fraction_sc1
            amplBeamletTab *= 2.0**W_tab_fraction_sc1
        
        if 'arts_sc4' in apps:
            # Arts SC4 Power TAB for single polarization Stokes I = XX* + 0
            # . Fine channel voltage TAB
            fineChannelTabSamples *= 2.0**W_tab_voltage_fraction_sc4
            rmsComplexFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
            rmsFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
            amplFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
            
            # . Fine channel power TAB
            fineChannelTabPowers *= 2.0**W_tab_power_fraction_sc4
            meanFineChannelTabPower *= 2.0**W_tab_power_fraction_sc4
            stdFineChannelTabPower *= 2.0**W_tab_power_fraction_sc4
            
            # . Integrated channels power TAB
            integratedTabPowers *= 2.0**W_tab_power_int_fraction_sc4
            meanIntegratedTabPower *= 2.0**W_tab_power_int_fraction_sc4
            stdIntegratedTabPower *= 2.0**W_tab_power_int_fraction_sc4
            
            # Arts SC4 IAB for single polarization Stokes I = XX* + 0
            # . Fine channel powers
            fineChannelPowers *= 2.0**W_iab_power_fraction
            meanFineChannelPower *= 2.0**W_iab_power_fraction
            
            # . Integrated channel powers
            integratedChannelPowers *= 2.0**W_iab_power_fraction
            meanIntegratedChannelPower *= 2.0**W_iab_power_fraction
            
            # . Sum of powers is coherent sum
            integratedIabPowers *= 2.0**W_iab_power_fraction
            meanIntegratedIabPower *= 2.0**W_iab_power_fraction
            
            # . Output gain
            outputIabPowers *= 2.0**W_iab_gain_fraction
            meanOutputIabPower *= 2.0**W_iab_gain_fraction
            stdOutputIabPower *= 2.0**W_iab_gain_fraction
    
    ################################################################################
    # Logging
    ################################################################################    
    print '--------------------------------------------------------------------'
    print '-- User settings'
    print '--------------------------------------------------------------------'
    if wgEnable:
        print 'WG sinus input:'
        print '. wgSigma              = %f' % wgSigma
        print '. wgAmpl               = %f [fsub]' % wgAmpl
        print '. wgSweep              = %f [fsub]' % wgSweep
    else:
        print 'ADC noise input:'
        print '. skySigma             = %f' % skySigma
    print '. quantize             = %s' % quantize
    print '. nofDishes            = %d   (= %.1f bit)' % (nofDishes, math.log(nofDishes, 2))
    print ''
    
    if 'dynamic' in models:
        print 'Dynamic model:'
        print '. selSub               = %d' % selSub
        print '. selChanx             = %d' % selChanx
        print '. selChanxFraction     = %d' % selChanxFraction
        print '. nofTchanx            = %d' % nofTchanx
        print '. nofTsub              = %d' % nofTsub
        print '. nofTs                = %d' % nofTs
        print ''
    
    print '--------------------------------------------------------------------'
    print '-- Design settings'
    print '--------------------------------------------------------------------'
    print 'Apertif BF:'
    print '. N_fft                = %d' % N_fft
    print '. N_sub                = %d' % N_sub
    print '. N_int_sub            = %d' % N_int_sub
    print ''
    if 'apertif' in apps:
        print 'Apertif XC:'
        print '. N_chan_x             = %d' % N_chan_x
        print '. N_chan_bf            = %d' % N_chan_bf
        print '. N_int_chan_x         = %d' % N_int_chan_x
        print '. cbWeightAdjust       = %f' % cbWeightAdjust
        print '. cbMaxWeight          = %d' % pi_apr.cbMaxWeight
        print '. cbUnitWeight         = %d' % pi_apr.cbUnitWeight
        print '. cbWeight             = %d' % cbWeight
        print '. cbGain               = %f' % cbGain
        print ''
    if 'arts_sc1' in apps:
        print 'Arts SC1:'
        print '. sc1_use_qua          = %s' % sc1_use_qua
        print '. tabWeightAdjustSc1   = %f' % tabWeightAdjustSc1
        print '. tabUnitWeightSc1     = %d' % tabUnitWeightSc1
        print '. tabWeightSc1         = %d' % tabWeightSc1
        print '. tabGainSc1           = %.8f' % tabGainSc1
        print ''
    if 'arts_sc4' in apps:
        print 'Arts SC4:'
        print '. tabWeightAdjustSc4   = %f' % tabWeightAdjustSc4
        print '. iabGainAdjust        = %f' % iabGainAdjust
        print '. iabMaxGain           = %d' % pi_apr.iabMaxGain
        print '. iabUnitGain          = %d' % pi_apr.iabUnitGain
        print '. iabGainReg           = %d' % iabGainReg
        print '. iabGain              = %.8f' % iabGain
        print ''
    
    if 'dynamic' in models:
        print '--------------------------------------------------------------------'
        print '-- Dynamic model settings'
        print '--------------------------------------------------------------------'
        print 'Design:'
        print '. fs                           = %f [Hz]' % fs
        print '. Ts                           = %.2f [ns]' % (Ts * 1e9)
        print '. fsub                         = %f [Hz]' % fsub
        print '. Tsub                         = %.2f [us]' % (Tsub * 1e6)
        print '. W_fsub_gain                  = %d [bit]' % pi_apr.W_fsub_gain
        if 'apertif' in apps:
            print '. fchan_x                      = %f [Hz]' % fchan_x
            print '. Tchanx                       = %.2f [us]' % (Tchanx * 1e6)
            print '. W_fchan_x_gain               = %d [bit]' % pi_apr.W_fchan_x_gain
        if 'arts_sc4' in apps:
            print '. Tstokesx                     = %.2f [us]' % (Tstokesx * 1e6)
        print ''
    
        print 'Dish:'
        print '. W_adc                        = %d [bit]' % pi_apr.W_adc
        print '. W_adc_max                    = %d [bit], adc max = %d' % (pi_apr.W_adc_max, 2**pi_apr.W_adc_max-1)
        print '. W_sst_in                     = %d [bit]' % pi_apr.W_sst_in
        print '. W_sst_fraction               = %d [bit]' % pi_apr.W_sst_fraction
        print '. W_cb_weight                  = %d [bit]' % pi_apr.W_cb_weight
        print '. W_cb_unit_scale              = %d [bit]' % pi_apr.W_cb_unit_scale
        print '. W_cb_weight_fraction         = %d [bit]' % pi_apr.W_cb_weight_fraction
        print '. W_cb_in                      = %d [bit]' % pi_apr.W_cb_in
        print '. W_cb_product                 = %d [bit]' % W_cb_product
        print '. W_bst_in                     = %d [bit]' % pi_apr.W_bst_in
        print '. W_bst_fraction               = %d [bit]' % pi_apr.W_bst_fraction
        print '. W_beamlet_sc1                = %d [bit]' % pi_apr.W_beamlet_sc1
        print '. W_beamlet                    = %d [bit]' % pi_apr.W_beamlet
        print '. W_beamlet_fraction           = %d [bit]' % pi_apr.W_beamlet_fraction
        print ''
    
        if 'apertif' in apps:
            print 'Apertif XC:'
            print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
            print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
            print ''
            
        if 'arts_sc1' in apps:
            print 'Arts SC1 beamlet TAB:'
            print '. W_tab_weight_sc1             = %d [bit]' % pi_apr.W_tab_weight_sc1
            print '. W_tab_unit_scale_sc1         = %d [bit]' % W_tab_unit_scale_sc1
            print '. W_tab_weight_fraction_sc1    = %d [bit]' % W_tab_weight_fraction_sc1
            print '. W_beamlet_sc1                = %d [bit]' % pi_apr.W_beamlet_sc1
            print '. W_tab_product_sc1            = %d [bit]' % W_tab_product_sc1
            print '. W_tab_fraction_sc1           = %d [bit]' % W_tab_fraction_sc1
            print ''
        
        if 'arts_sc4' in apps:
            print 'Arts SC4 TAB:'
            print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
            print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
            print '. W_tab_weight_sc4             = %d [bit]' % pi_apr.W_tab_weight_sc4
            print '. W_tab_unit_scale_sc4         = %d [bit]' % pi_apr.W_tab_unit_scale_sc4
            print '. W_tab_weight_fraction_sc4    = %d [bit]' % pi_apr.W_tab_weight_fraction_sc4
            print '. tabUnitWeightSc4             = %d' % pi_apr.tabUnitWeightSc4
            print '. tabWeightSc4                 = %d' % tabWeightSc4
            print '. tabGainSc4                   = %.8f' % tabGainSc4
            print '. W_beamlet                    = %d [bit]' % pi_apr.W_beamlet
            print '. W_tab_voltage_product_sc4    = %d [bit]' % W_tab_voltage_product_sc4
            print '. W_tab_voltage_fraction_sc4   = %d [bit]' % W_tab_voltage_fraction_sc4
            print '. W_tab_power_product_sc4      = %d [bit]' % W_tab_power_product_sc4
            print '. W_tab_power_fraction_sc4     = %d [bit]' % W_tab_power_fraction_sc4
            print '. W_tab_power_int_fraction_sc4 = %d [bit]' % W_tab_power_int_fraction_sc4
            print '. W_int_ax                     = %d [bit]' % W_int_ax
            print '. M_int_ax                     = %d' % M_int_ax
            print ''
            
            print 'Arts SC4 IAB:'
            print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
            print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
            print '. W_int_ax                     = %d [bit]' % W_int_ax
            print '. M_int_ax                     = %d' % M_int_ax
            print '. W_iab_gain                   = %d [bit]' % pi_apr.W_iab_gain
            print '. W_iab_unit_scale             = %d [bit]' % pi_apr.W_iab_unit_scale
            print '. W_iab_gain_fraction          = %d [bit]' % pi_apr.W_iab_gain_fraction
            print '. W_iab_power_product          = %d [bit]' % W_iab_power_product
            print '. W_iab_power_fraction         = %d [bit]' % W_iab_power_fraction
            print '. W_iab_gain_product           = %d [bit]' % W_iab_gain_product
            print '. W_iab_gain_fraction          = %d [bit]' % W_iab_gain_fraction
            print ''
    
    if 'dynamic' in models and wgSweep == 0:
        print '--------------------------------------------------------------------'
        print '-- Dynamic model results'
        print '--------------------------------------------------------------------'
        if wgEnable:
            print 'Apertif BF internal levels (excluding sign bit) for coherent WG sinus input:'
            print '. ADC input rms                = %20.2f               (= %4.1f [bit])' % (rmsSp, math.log(rmsSp, 2))
            print '. ADC input ampl               = %20.2f               (= %4.1f [bit])' % (amplSp, math.log(amplSp, 2))
            print '. SST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexSST, math.log(rmsComplexSST, 2))
            print '. SST input ampl               = %20.2f               (= %4.1f [bit])' % (amplSST, math.log(amplSST, 2))
            print '. BST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexBST, math.log(rmsComplexBST, 2))
            print '. BST input ampl               = %20.2f               (= %4.1f [bit])' % (amplBST, math.log(amplBST, 2))
            print '. CB beamlet rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamlet, math.log(rmsComplexBeamlet, 2))
            print '. CB beamlet ampl              = %20.2f               (= %4.1f [bit])' % (amplBeamlet, math.log(amplBeamlet, 2))
        else:
            print 'Apertif BF internal levels (excluding sign bit) for incoherent ADC noise input:'
            print '. ADC input rms                = %20.2f               (= %4.1f [bit])' % (rmsSp, math.log(rmsSp, 2))
            print '. SST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexSST, math.log(rmsComplexSST, 2))
            print '. SST input rms                = %20.2f               (= %4.1f [bit])' % (rmsSST, math.log(rmsSST, 2))
            print '. BST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexBST, math.log(rmsComplexBST, 2))
            print '. BST input rms                = %20.2f               (= %4.1f [bit])' % (rmsBST, math.log(rmsBST, 2))
            print '. CB beamlet rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamlet, math.log(rmsComplexBeamlet, 2))
            print '. CB beamlet rms               = %20.2f               (= %4.1f [bit])' % (rmsBeamlet, math.log(rmsBeamlet, 2))
        print '. SST power                    = %20.2f  = %7.2f dB (= %4.1f [bit])' % (SSTpower, 10*math.log(SSTpower, 10), math.log(SSTpower, 2))
        print '. BST power                    = %20.2f  = %7.2f dB (= %4.1f [bit])' % (BSTpower, 10*math.log(BSTpower, 10), math.log(BSTpower, 2))
        print ''
        
        if 'apertif' in apps:
            if wgEnable:
                print 'Apertif XC output for coherent WG sinus input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
            else:
                print 'Apertif XC output for incoherent ADC noise input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
            print '. XC auto power                = %20.2f  = %7.2f dB (= %4.1f [bit])' % (XCpower, 10*math.log(XCpower, 10), math.log(XCpower, 2))
            print ''
        
        if 'arts_sc1' in apps:
            if wgEnable:
                print 'Arts SC1 beamlet TAB internal levels (excluding sign bit) for coherent WG sinus input:'
                print '. Output TAB rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamletTab, math.log(rmsComplexBeamletTab, 2))
                print '. Output TAB ampl              = %20.2f               (= %4.1f [bit])' % (amplBeamletTab, math.log(amplBeamletTab, 2))
            else:
                print 'Arts SC1 beamlet TAB internal levels (excluding sign bit) for incoherent ADC noise input:'
                print '. Output TAB rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamletTab, math.log(rmsComplexBeamletTab, 2))
                print '. Output TAB rms               = %20.2f               (= %4.1f [bit])' % (rmsBeamletTab, math.log(rmsBeamletTab, 2))
            print ''
        
        if 'arts_sc4' in apps:
            if wgEnable:
                print 'Arts SC4 TAB internal levels (excluding sign bit) for coherent WG sinus input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
                print '. Fine channel TAB rms complex = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannelTab, math.log(rmsComplexFineChannelTab, 2))
                print '. Fine channel TAB ampl        = %20.2f               (= %4.1f [bit])' % (amplFineChannelTab, math.log(amplFineChannelTab, 2))
                print '. Fine channel TAB power       = %20.2f               (= %4.1f [bit])' % (meanFineChannelTabPower, math.log(meanFineChannelTabPower, 2))
                print '. Output TAB power             = %20.2f               (= %4.1f [bit])' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2))
            else:
                print 'Arts SC4 TAB internal levels (excluding sign bit) for incoherent ADC noise input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
                print '. Fine channel TAB rms complex = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannelTab, math.log(rmsComplexFineChannelTab, 2))
                print '. Fine channel TAB rms         = %20.2f               (= %4.1f [bit])' % (rmsFineChannelTab, math.log(rmsFineChannelTab, 2))
                print '. Fine channel TAB power       = %20.2f               (= %4.1f [bit])' % (meanFineChannelTabPower, math.log(meanFineChannelTabPower, 2))
                print '. Fine channel TAB power std   = %20.2f               (= %4.1f [bit])' % (stdFineChannelTabPower, math.log(stdFineChannelTabPower, 2))
                print '. Output TAB power             = %20.2f               (= %4.1f [bit])' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2))
                print '. Output TAB power std         = %20.2f               (= %4.1f [bit])' % (stdIntegratedTabPower, math.log(stdIntegratedTabPower, 2))
            print ''
            
            if wgEnable:
                print 'Arts SC4 IAB internal levels (excluding sign bit) for coherent WG sinus input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
                print '. Fine channel power           = %20.2f               (= %4.1f [bit])' % (meanFineChannelPower, math.log(meanFineChannelPower, 2))
                print '. Integrated channel power     = %20.2f               (= %4.1f [bit])' % (meanIntegratedChannelPower, math.log(meanIntegratedChannelPower, 2))
                print '. Integrated IAB power         = %20.2f               (= %4.1f [bit])' % (meanIntegratedIabPower, math.log(meanIntegratedIabPower, 2))
                print '. Output IAB power             = %20.2f               (= %4.1f [bit])' % (meanOutputIabPower, math.log(meanOutputIabPower, 2))
            else:
                print 'Arts SC4 IAB internal levels (excluding sign bit) for incoherent ADC noise input:'
                print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
                print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
                print '. Fine channel power           = %20.2f               (= %4.1f [bit])' % (meanFineChannelPower, math.log(meanFineChannelPower, 2))
                print '. Integrated channel power     = %20.2f               (= %4.1f [bit])' % (meanIntegratedChannelPower, math.log(meanIntegratedChannelPower, 2))
                print '. Integrated IAB power         = %20.2f               (= %4.1f [bit])' % (meanIntegratedIabPower, math.log(meanIntegratedIabPower, 2))
                print '. Output IAB power             = %20.2f               (= %4.1f [bit])' % (meanOutputIabPower, math.log(meanOutputIabPower, 2))
                print '. Output IAB power std         = %20.2f               (= %4.1f [bit])' % (stdOutputIabPower, math.log(stdOutputIabPower, 2))
            print ''
    
    ################################################################################
    # Plotting
    ################################################################################    
    if useplot:
        tplot = test_plot.Testplot()
        tplot.close_plots()
        
        if 'static' in models:
            ################################################################################
            # Static model:
            # . WG, coherent input
            # . Noise, incoherent input
            ################################################################################
            
            if wgEnable:
                # WG, coherent input
                # Plot internal levels as function of WG amplitude
                Lx = A_wg
                if 'apertif' in apps or 'arts_sc4' in apps:
                    A = []
                    A.append(A_wg)
                    A.append(beamlet_wg_complex_phase_values * pi_apr.wgScale)
                    A.append(channel_x_wg_complex_phase_values * pi_apr.wgScale)
                    Alegend = ['adc', 'beamlet', 'channel_x']
                    lineFormats = ['-', '-', '-']
                    Title = 'Amplitudes for cbWeight=%d' % cbWeight
                    Ylim = None
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_amplitudes.png')   # Save current figure
    
                if 'arts_sc1' in apps:
                    A = []
                    A.append(A_wg)
                    A.append(beamlet_wg_complex_phase_values * pi_apr.wgScale)
                    A.append(beamlet_tab_wg_complex_phase_values_sc1 * pi_apr.wgScale)
                    Alegend = ['adc', 'beamlet', 'tab sc1']
                    lineFormats = ['-', '-', '*']
                    Title = 'Amplitudes for cbWeight=%d, tabWeightSc1=%d, nofDishes=%d' % (cbWeight, tabWeightSc1, nofDishes)
                    Ylim = None
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_amplitudes_sc1.png')   # Save current figure
            
                if 'arts_sc4' in apps:
                    A = []
                    A.append(power_tab_xx_wg_complex_values_sc4)
                    A.append(iab_xx_wg_complex_values_sc4)
                    Alegend = ['TAB XX* power', 'IAB XX* power']
                    lineFormats = ['-', '*']
                    Title = 'Powers for cbWeight=%d, tabWeightSc4=%d, iabGainReg=%d, nofDishes=%d' % (cbWeight, tabWeightSc4, iabGainReg, nofDishes)
                    Ylim = [0, 1.1 * 2.0**pi_apr.W_pow]
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_powers_sc4.png')   # Save current figure
                
                # Plot power statistics as function of WG amplitude
                A = []
                A.append(SST_wg_dBs)
                A.append(BST_wg_dBs)
                if 'apertif' in apps:
                    A.append(XC_wg_auto_visibilities_dB)
                    Alegend = ['SST', 'BST', 'XC']
                else:
                    Alegend = ['SST', 'BST']
                Title = 'Auto correlation power statistics for CB weight = %d' % cbWeight
                Ylim = None
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level [dB]', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_auto_powers.png')   # Save current figure
            
            else:
                # Noise, incoherent input
                # Plot internal levels as function of WG amplitude
                Lx = S_noise
                if 'apertif' in apps or 'arts_sc4' in apps:
                    A = []
                    A.append(S_noise)
                    A.append(beamlet_noise_complex_phase_values)
                    A.append(channel_x_noise_complex_phase_values)
                    Alegend = ['adc', 'beamlet', 'channel_x']
                    lineFormats = ['-', '*', '-']
                    Title = 'Sigmas for cbWeight=%d' % cbWeight
                    Ylim = None
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_sigmas.png')   # Save current figure
            
                if 'arts_sc1' in apps:
                    A = []
                    A.append(S_noise)
                    A.append(beamlet_noise_complex_phase_values)
                    A.append(beamlet_tab_noise_complex_phase_values_sc1)
                    Alegend = ['adc', 'beamlet', 'tab sc1']
                    lineFormats = ['-', '-', '*']
                    Title = 'Sigmas for cbWeight=%d, tabWeightSc1=%d, nofDishes=%d' % (cbWeight, tabWeightSc1, nofDishes)
                    Ylim = None
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_sigmas_sc1.png')   # Save current figure
    
                if 'arts_sc4' in apps:
                    A = []
                    A.append(power_tab_xx_noise_values_sc4)
                    A.append(iab_xx_noise_values_sc4)
                    Alegend = ['TAB XX* power', 'IAB XX* power']
                    lineFormats = ['-', '*']
                    Title = 'Powers for cbWeight=%d, tabWeightSc4=%d, iabGainReg=%d, nofDishes=%d' % (cbWeight, tabWeightSc4, iabGainReg, nofDishes)
                    Ylim = [0, 1.1 * 2.0**pi_apr.W_pow]
                    tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
                    tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_powers_sc4.png')   # Save current figure
                
                # Plot power statistics as function of WG amplitude
                A = []
                A.append(SST_noise_dBs)
                A.append(BST_noise_dBs)
                A.append(XC_noise_auto_visibilities_dB)
                Alegend = ['SST', 'BST', 'XC']
                Title = 'Auto correlation power statistics for CB weight = %d' % cbWeight
                Ylim = None
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level [dB]', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_auto_powers.png')   # Save current figure
    
            
        if 'dynamic' in models and wgSweep == 0:
            ###############################################################################
            # Dynamic model: Processing
            ###############################################################################
        
            # Plot dish ADC
            L = spSamples[0:nof_periods]  # do not plot all to save time
            if wgEnable:
                Title = 'Dish ADC input (sigma = %.2f, ampl = %.2f)' % (rmsSp, wgAmpl)
                ylim = 5 * math.ceil(0.2 * 1.1 * amplSp)
            else:
                Title = 'ADC input (sigma = %.2f)' % rmsSp
                ylim = 5 * math.ceil(0.2 * 1.1 * rmsSp * nofSigma)
            Ylim = [-ylim, ylim]
            tplot.plot_one_dimensional_list(3, L, Lx=None, representation='', offset=0, Title=Title, Xlabel='Time [Tsub]', Ylabel='Voltage', Xlim=None, Ylim=Ylim)
            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_adc.png')   # Save current figure
            
            # Plot dish output beamlets for selSub
            A = []
            A.append(beamletSamples.real[0:nof_periods])
            A.append(beamletSamples.imag[0:nof_periods])
            Alegend = ['real', 'imag']
            if wgEnable:
                Title = 'Dish CB output samples (ampl = %.2f)' % amplBeamlet
                ylim = 5 * math.ceil(0.2 * 1.1 * amplBeamlet)
            else:
                Title = 'Dish CB output samples (sigma = %.2f)' % rmsBeamlet
                ylim = 5 * math.ceil(0.2 * 1.1 * rmsBeamlet * nofSigma)
            Ylim = [-ylim, ylim]
            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_beamlet.png')   # Save current figure
            
            if 'apertif' in apps or 'arts_sc4' in apps:
                # Plot correlator fine channels for selSub, selChanx
                A = []
                A.append(fineChannelSamples.real[0:nof_periods])
                A.append(fineChannelSamples.imag[0:nof_periods])
                Alegend = ['real', 'imag']
                if wgEnable:
                    Title = 'Central fine channel samples (ampl = %.2f)' % amplFineChannel
                    ylim = 5 * math.ceil(0.2 * 1.1 * amplFineChannel)
                else:
                    Title = 'Central fine channel samples (sigma = %.2f)' % rmsFineChannel
                    ylim = 5 * math.ceil(0.2 * 1.1 * rmsFineChannel * nofSigma)
                Ylim = [-ylim, ylim]
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_fine_channel.png')   # Save current figure
            
            if 'arts_sc1' in apps:
                # Plot Arts SC1 output TAB samples for selSub and nofDishes
                # . Independent of nofDishes when tabWeightSc1 is scaled by nofDishes or sqrt(nofDishes)
                A = []
                A.append(beamletTabSamples.real[0:nof_periods])
                A.append(beamletTabSamples.imag[0:nof_periods])
                Alegend = ['real', 'imag']
                if wgEnable:
                    Title = 'Arts SC1 output TAB samples (ampl = %.2f, nofDishes = %d)' % (amplBeamletTab, nofDishes)
                    ylim = 5 * math.ceil(0.2 * 1.1 * amplBeamletTab)
                else:
                    Title = 'Arts SC1 output TAB samples (sigma = %.2f, nofDishes = %d)' % (rmsBeamletTab, nofDishes)
                    ylim = 5 * math.ceil(0.2 * 1.1 * rmsBeamletTab * nofSigma)
                Ylim = [-ylim, ylim]
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_beamlet_tab.png')   # Save current figure
        
            if 'arts_sc4' in apps:
                # Plot Arts SC4 output TAB powers for selSub, selChanx and nofDishes
                # . Independent of nofDishes, because tabWeightSc4 is scaled by nofDishes or sqrt(nofDishes)
                A = []
                A.append(integratedTabPowers[0:nof_periods])
                Alegend = None
                Title = 'Arts SC4 output TAB powers (mean = %.2f = %.1f [bit], nofDishes = %d)' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2), nofDishes)
                if wgEnable:
                    ylim = 5 * math.ceil(0.2 * 1.1 * meanIntegratedTabPower)
                else:
                    ylim = 5 * math.ceil(0.2 * 1.1 * meanIntegratedTabPower * nofSigma)
                Ylim = [0, ylim]
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tchanx]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_power_tab.png')   # Save current figure
                
                # Plot Arts SC4 output IAB powers for selSub, selChanx and nofDishes
                # . Independent of nofDishes, because tabWeightSc4 is scaled by nofDishes or sqrt(nofDishes)
                A = []
                A.append(outputIabPowers)
                Alegend = None
                Title = 'Arts SC4 output IAB powers (mean = %.2f = %.1f [bit], nofDishes = %d)' % (meanOutputIabPower, math.log(meanOutputIabPower, 2), nofDishes)
                if wgEnable:
                    ylim = 5 * math.ceil(0.2 * 1.1 * meanOutputIabPower)
                else:
                    ylim = 5 * math.ceil(0.2 * 1.1 * meanOutputIabPower * nofSigma)
                Ylim = [0, ylim]
                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tchanx]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
                tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_power_iab.png')   # Save current figure
            
            
        if 'dynamic' in models and wgSweep != 0:
            ###############################################################################
            # Dynamic model: WG frequency sweep
            ###############################################################################
    
            # Plot subband spectrum abs(fSubSamples[nofTsub][N_sub])
            A = np.abs(fSubSamples)
            tplot.image_two_dimensional_list(3, A, transpose=True, grid=True, Title='Subband magnitudes', Xlabel='time [Tsub]', Ylabel='freq [fsub]')
            tplot.save_figure('plots/apertif_arts_firmware_model_wg_sweep_subband_spectrum.png')   # Save current figure
            
            # Plot fine channels spectrum abs(fSubChanSamples[N_sub][nofTchanx][N_chan_x])
            A = np.abs(fSubChanSamples.reshape(nofTchanx, N_sub*N_chan_x))   # [nofTchanx][N_sub*N_chan_x]
            chanLo = 0
            chanHi = N_sub*N_chan_x - 1
            nofSubBorder = 0
            nofSubBorder = 1 + int(abs(wgSweep))
            if nofSubBorder>0:
                if wgSweep>0:
                    chanLo = N_chan_x * (selSub - nofSubBorder)
                    chanHi = N_chan_x * (selSub + nofSubBorder + wgSweep + 1) - 1
                else:
                    chanLo = N_chan_x * (selSub - nofSubBorder + wgSweep)
                    chanHi = N_chan_x * (selSub + nofSubBorder + 1) - 1
                A = A[:, chanLo:chanHi+1]
            extent = [0, nofTchanx-1, 1.0 * chanLo / N_chan_x - 0.5, 1.0 * chanHi / N_chan_x - 0.5]
            #print chanLo, chanHi, extent
            #extent = None
            tplot.image_two_dimensional_list(3, A, transpose=True, grid=True, Title='Channel magnitudes', Xlabel='time [Tchanx]', Ylabel='freq [fsub]', extent=extent)
            tplot.save_figure('plots/apertif_arts_firmware_model_wg_sweep_fine_channel_spectrum.png')   # Save current figure
        
        tplot.show_plots()