Skip to content

GitLab

Commit bf1f35cb authored by Daniel van der Schuur's avatar Daniel van der Schuur
Browse files
Options

Copied Matlab dir from SVN. Added readme.md

parents
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 3664 additions and 0 deletions
matlab/BFweights.m 0 → 100755
View file @ bf1f35cb
%% BF weights for TABs ARTS
%
% This script implements the equation on the third line on slide 16 of the
% presentation "Beamforming for ARTS" given by Stefan Wijnholds during an
% ARTS project meeting on April 4, 2018. The script produces a matrix
% containing the phases of the beamformer weights of size Ndish x NTAB,
% where Ndish is the number of dishes and NTAB is the number of TABs. It is
% assumed that the dishes for a Uniform Linear Array and that the TABs are
% equidistantly spaced between center of the compound beam and the first
% grating response. It is also assumed that the delay is already
% compensated for the center of the compound beam.
%
% SJW, 4 April 2018
%% start with a clean workspace
clear
close all
%% calculate phases
Ndish = 10;
NTAB = 12;
phase = zeros(Ndish, NTAB);
% determine beta, taking into account the TAB indexing conventions for ARTS
TABidx = -floor(NTAB/2):floor((NTAB-1)/2);
beta = TABidx / NTAB;
% define dish locations measured in integer multiples of the common
% quotient baseline
n = 0:(Ndish-1);
% calculate phases
phase = 2 * pi * n.' * beta;
phase_n = n.' * beta;
% show result
imagesc(phase_n)
colorbar
title(['BF weight phases * 2*pi [rad]'])
xlabel(['TABs']);
ylabel(['Dishes']);
matlab/FIR_LPF_ApertifCF.m 0 → 100644
View file @ bf1f35cb
%-----------------------------------------------------------------------------
%
% Copyright (C) 2016
% 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: R. de Wild, 2015 (Original)
% E. Kooistra, 2016
%
% Purpose : Calculate FIR filter coefficients for Apertif channel
% filterbank using fircls1.
% Description :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% workspace initiation
clear all
close all
fig=0;
% preset FIR-filter parameters
K = 64; % Number of channels
Taps = 8; % Number of taps, even 4 could be sufficient (erko)
M = Taps*K; % Number of filter coefficients (taps)
h_ix = 1:M; % coefficient index
nbit = 9; % bits/coefficient = 1 sign + (nbit-1) mantissa
fs = 781250; % sampling clock = 800.0 MHz / 1024
f_nyq = fs/2; % Nyquist frequency
f_pb = fs/K/2; % end of passband
f_chan = fs/K;
% normalized frequency axis & normalized amplitude axis
w_pb = f_pb/f_nyq; % = f_chan
% The fircls1 pass band deviation and stop band deviation are positive
% fractions of 1. The ripple in dB needs to be divided by 10 to get from
% dB = 10*log10() the log10() and divided by 2 because it is a relative
% voltage value (instead of a relative power value).
% When the pass band ripple in dB is positive then the linear factor is
% sligthly larger than 1, so subtract 1 to get the pass band deviation.
% The pass band deviation defines the +- fraction around 1. The stop band
% 'ripple' in dB is negative and defines the attenuation. The stop band
% deviation is the fraction of 1 that remains in the stop band. Therefore
% both the pass band ripple and the stop band attenuation result both in
% linear deviations that are close to 0.
% The fircls1 description defines DEVP is the maximum passband deviation
% or ripple (in linear units) from 1, and DEVS is the maximum stopband
% deviation or ripple (in linear units) from 0.
r_pb = 10^(0.025)-1; % passband ripple = 0.5 dB
r_sb = 10^(-3.0); % stopband ripple = -60.0 dB
% compute M filter-coefficients
h = fircls1(M-1, w_pb, r_pb, r_sb);
hfs_abs = abs(fftshift(fft(h ,M)));
% compute quantized filter coefficients
hmax = max(h);
hq = double( uencode(h, nbit, hmax, 'signed') );
hqfs_abs = abs(fftshift (fft( hq ,M))) * sum(h)/sum(hq); % normalize to response for DC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Print FIR-filter parameters
sprintf([ ...
' fs = sample rate = %6.2f [Hz] \n' , ...
' f_chan = channel bandwidth = %6.2f [Hz]\n' , ...
' r_pb = pass band deviation = %9.7f = %6.2f dB_W ripple\n' , ...
' r_sb = stop band deviation = %9.7f = %6.2f dB_W attenuation\n' , ...
' M = number of coefficients = %d'], ...
fs, f_chan, r_pb, 20*log10(1+r_pb), r_sb, 20*log10(r_sb), M)
% Save FIR-filter-coefficients
% N.B. - read '.txt'-file with 'wordpad' or use \r\n !!!
file_name = ['FIR_LPF_ApertifCF.txt'];
FIRtable = [ h(h_ix) ] ;
fid = fopen(file_name , 'w');
fprintf(fid ,'%14.10f \n', FIRtable);
fclose(fid);
% Plot FIR-filter coefficients
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(hq);
title(['FIR filter coefficients (', num2str(nbit), ' bits)'] );
grid on;
% Plot FIR-filter amplitude characteristic
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
f = (h_ix - M/2-1) * fs/M / f_chan;
plot ( f, 20.0 * log10 ( hfs_abs ), 'r-', ...
f, 20.0 * log10 ( hqfs_abs ), 'k-') ;
xlo = f(1);
xhi = f(M);
xlo = -3;
xhi = 3;
xlim( [xlo xhi] );
grid on;
zoom on;
title('FIR filter transfer function');
text( xlo, -5, ['sampling rate = ', num2str(fs), ' [Hz]']);
text( xlo, -10, ['channel bandwidth = ', num2str(f_chan), ' [Hz]']);
text( xlo, -15, ['number of taps = ', num2str(Taps)]);
text( xlo, -20, ['number of channels = ', num2str(K)]);
text( xlo, -25, ['number of coefficients = ', num2str(M)]);
text( xlo, -30, ['number of bits/coefficient = ', num2str(nbit)]);
xlabel('Frequency [channels]');
ylabel('Power [dB]');
print -dpng 'arts_fir';
matlab/apertif_arts_firmware_model.py 0 → 100644
View file @ bf1f35cb
This diff is collapsed.
matlab/apertif_matlab_readme.txt 0 → 100644
View file @ bf1f35cb
This directory contains Matlab code for modelling DSP aspects of Apertif:
- Functions are used in scripts.
- Scripts can run in Matlab
- Figure files are saved figure plots of scripts
- Data files are input read by a scripts, or output saved by a the script
1) Fringe stopping:
fringe.m Script Plot the fringe due to the varying difference in geometrical delay.
fringe.jpg Figure Fringe plot saved by fringe.m (figure in ASTRON-MEM-199)
phase_tracking.m Script Plot the phase tracking (PT) phase and phase error as function of the
varying geometrical delay (see ASTRON-MEM-199).
phase_tracking_phases.txt Data Input reference data (the PT piecewise linear coefficients) for HDL
simulation saved by phase_tracking.m
phase_tracking_coefficients.txt Data Output reference data (the PT expected phases) for HDL simulation
saved by phase_tracking.m
phase_lookup.m Script Plot and save lookup table for conversion from phi to pair of re, im
2) Quantization
try_uencode.m Script Try Matlab uencode to determine optimum quantisation, better use
quantize.m
quantize.m Function Quantize signed input data
try_quantize.m Script Run quantize.m
compensation_and_selection_gain.m Script Determine gain_w and gain value for optimum requantisation
try_power_requantize.py Python Investigate requantization of data and auto power for e.g. Stokes I
and IQUV in Arts
apertif_arts_firmware_model.py Python Model fixed point signal levels in Apertif and Arts, see
../doc/ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf
3) Low pass FIR filter coefficients
See also:
- $RADIOHDL/libraries/dsp/doc/filterbank.txt
- $RADIOHDL/libraries/dsp/filter/src/python/ for diff* FIR filter files and recreate* FIR files from dat
to mif/hex format
Fsub original from LOFAR:
- data/Coeffs16384Kaiser-quant.dat Data Original PFIR filter coefficients from LOFAR 1.0 subband filterbank
(16 tap, 1024 point FFT), the script that created these PFIR
coefficients is not available anymore (Erko looked in old
repositories back to 2005 of Lofar software and Station firmware
but could not find it.)
- data/Coeffs16384Kaiser-quant-withdc.dat = data/Coeffs16384Kaiser-quant.dat
Fsub adjusted to have no DC in PFIR coefficients per polyphase
- data/Coeffs16384Kaiser-quant-nodc.dat = Used in Apertif Fsub obtained from data/Coeffs16384Kaiser-quant.dat
using run_pfir_coeff.m -r 18333 as commited, so with application =
'lofar_subband', which selects:
. config.hp_adjust = false
. config.dc_adjust = true
Fsub bypass (equivalent to not having a PFIR section, useful to isolate issues between PFB and only FFT)
- data/run_pfir_coeff_m_bypass_16taps_1024points_16b.dat
Fsub bypass per polyphase, used to diagnose quantization and spurious issues per polyphase, but probably not
so useful
- data/run_pfir_coeff_m_bypass_16taps_1024points_16b_polyphase0.dat -- only polyphase 0 has coeff != 0
- data/run_pfir_coeff_m_bypass_16taps_1024points_16b_polyphase1.dat -- only polyphase 1 has coeff != 0
- data/run_pfir_coeff_m_bypass_16taps_1024points_16b_polyphase2.dat -- only polyphase 2 has coeff != 0
- data/run_pfir_coeff_m_bypass_16taps_1024points_16b_polyphase3.dat -- only polyphase 3 has coeff != 0
Fchan bypass (equivalent to not having a PFIR section, useful to isolate issues between PFB and only FFT)
- data/run_pfir_coeff_m_bypass_8taps_32points_9b.dat = for Arts SC4: N = 32, L = 8, coeff_w = 9.
- data/run_pfir_coeff_m_bypass_8taps_64points_9b.dat = for Apertif
bypass PFIR for Fchan using ones at tap 0 and zeros at other taps.
Obtained using run_pfir_coeff.m -r 18333, but with application =
'test_bypass' which selects N = 64, L = 8, coeff_w = 9.
Fchan with half power at +-fchan/2
- data/run_pfir_coeff_m_flat_hp_fircls1_8taps_32points_9b.dat = for Arts SC4
- data/run_pfir_coeff_m_flat_hp_fircls1_8taps_64points_9b.dat = for Apertif
PFIR for Fchan in Apertif obtained using run_pfir_coeff.m -r 18333
but with application = 'apertif_channel' which selects:
. config.hp_adjust = true
. config.dc_adjust = false
Fchan with half power at +-fchan/2 and no DC in PFIR coefficients per polyphase
- data/run_pfir_coeff_m_flat_hp_no_dc_fircls1_8taps_32points_9b.dat = for Arts SC4
- data/run_pfir_coeff_m_flat_hp_no_dc_fircls1_8taps_64points_9b.dat = for Apertif
PFIR for Fchan in Apertif obtained using run_pfir_coeff.m -r 18746
but with application = 'apertif_channel' which for -r 18746
selects:
. config.hp_adjust = true
. config.dc_adjust = true (was false in -r 18333)
FIR_LPF_ApertifCF.m Script Calculate FIR filter coefficients for Apertif channel filterbank
using fircls1 (original by Ronald de Wild).
FIR_LPF_ApertifCF.txt Data Floating point FIR coefficients saved by FIR_LPF_ApertifCF.m
arts_fir.png Figure FIR filter transfer function plot saved by FIR_LPF_ApertifCF.m
pfir_coeff.m Function Compute polyphase filterbank coefficients (same purpose as
FIR_LPF_ApertifCF.m but with more options)
pfir_coeff_adjust.m Function Compute polyphase filterbank coefficients with optionbal half power
bandwidth adjustment (calls pfir_coeff.m)
pfir_coeff_dc_adjust.m Function Apply DC adjustment to polyphase filterbank coefficients
run_pfir_coeff.m Script Run pfir_coeff_adjust.m and pfir_coeff.m to show filter response in
time and frequency domain, and optionally calls pfir_coeff_dc_adjust.m
4) FFT
sinusoids_and_fft_frequency_bins.m Script Code on DFT and DTFT from Matlab blog on internet
fft_sinus_noise_scaling.m Script to investigate sinus and noise levels of FFT input and output
try_fft.m Script Try DFT and DTFT
clk_spectrum Script Investigate EMI of clock and random data toggling
5) Data path functions
wg.m Function Waveform generator per block of data
dt.m Function Delay tracking per block of data
bsn_source.m Function Block Sequence Number source per block of data
pfir.m Function Polyphase FIR filter per block of data
pfft.m Function FFT per block of data
ds.m Function Doppler frequency shift per block of data
reorder_serial.m Function Select data from block of data
corner_turn.m Function Collect M blocks of in_data(1:K) and then transpose [M] and [K] for
corner turn
6) Polyphase filterbank
lofar/polyphase_demo_unquantized.m Script original by Jan Stemerdink
one_pfb.m Script Run data path model with subband polyphase filterbank (real input)
delay_tracking_pfb.m Script Based on one_pfb.m to simulate polyphase filterbank response on a delay step
psk_pfb.m Script Based on one_pfb.m to simulate polyphase filterbank response on PSK input
two_pfb.m Script Run data path model with subband polyphase filterbank followed by a channel
polyphase filterbank
7) HDL reference data
run_pfir.m Script Run the PFIR with real input to create HDL reference data (based on one_pfb.m)
run_pfft.m Script Run the PFFT with real input to create HDL reference data (based on one_pfb.m)
run_pfb.m Script Run the PFB with real input to create HDL reference data (based on one_pfb.m)
run_pfft_complex.m Script Run the PFFT with complex input to create HDL reference data (based on one_pfb.m)
run_pfb_complex.m Script Run the PFB with complex input to create HDL reference data (based on one_pfb.m)
8) Correlator
correlator_null_clip_or_wrap.m Script Investigate impact of input nulling, clipping or wrapping on correlator output
9) Beamformer
BFweights.m Script BF weights for TABs ARTS (original by Stefan Wijnholds).
beamformer_wrap_clip.m Script Compare wrapping or clipping of intermediate beamformer sum
\ No newline at end of file
matlab/beamformer_wrap_clip.m 0 → 100644
View file @ bf1f35cb
%-----------------------------------------------------------------------------
%
% Copyright (C) 2019
% 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: E. Kooistra, 2019
%
% Purpose : Model wrapping in the beamformer
%
% The internal beamformer sum can have significant bit growth. With 14 bit
% input at the ADC, log2(sqrt(1024)) = 5 bit extra for the PFB processing
% gain and log2(96) = 6.6 bit extra for the BF gain, the maximum internal
% beamformer sum is about 25 bit.
% For the final beamformer sum only w = 8 bits will be output. Therefor
% it seems a waste to still have to internally use 25 bits that require
% logic resources and IO resources on the ring. If the wrapped sum is used
% then less bits need to be used. Some extra bits are still needed to be
% able to distinguish wrapped final sums from correct final sums. Instead
% of 25 bits, e.g. use and transport 16 bit beamlet sums within the
% beamformer.
%
% The suppression of the beamformer outside the main beam makes that RFI
% outside the main beam may still fit in the w = 8 bit final sum.
% However intermediate beamformer sums can exceed this w = 8 bit output
% range. If the intermediate sum is clipped, then the output will be
% corrupted, because the clipping also flattens the weak astronomical
% signal that is on top of the RFI. If the intermediate sums are wrapped,
% then the information of the astronomical signal is preserved. In the
% final sum the RFI is then still attenuated and the astronomical signal
% remains undisturbed.
%
% The astronomical signal is weak in the sense that without RFI it fits
% in the w = 8 bit intermediate sum and final beamlet sum. The RFI signal
% is strong in the sense that it causes intermediate beamlet sums to
% overflow, but it is weak enough to still fits in the w = 8 bit beamlet
% output sum, provided that the RFI is not in the direction of the main
% beam.
clear all;
close all;
fig=0;
%% Settings
c = 3e8; % speed of light
N_theta= 1000; % Number of beam angles to evaluate over 90 degrees
N_ant = 96; % Number of antennas in array
f = 50e6; % RF frequency in LBA band
L = c/f; % RF wave length
b = 0.8*L; % distance between antennas, choose < L to have one main lobe
w = 8; % number of bits in beamlet
A = 3; % amplitude of weak astronomical signal, choose A << R
R = 30; % amplitude of strong RFI signal, choose R < 2**(w-1)
disp(sprintf(['. RF frequency = %5.1f MHz\n', ...
'. RF wave length = %5.1f m\n', ...
'. Antenna distance = %5.1f m\n', ...
'. Number of antennas = %d\n'], ...
f/1e6, ...
L, ...
b, ...
N_ant));
%% Random logic signal toggling at clock frequency fclk
s_period = 2^w;
s_mod = 2^(w-1);
s_max = s_mod-1;
s_min = -s_mod;
thetas = 90 * (0:N_theta-1)/N_theta;
DTs = b / c * sin(thetas*2*pi/360);
phis = 2*pi*f*DTs;
ants = 0:N_ant-1;
bf_sums_maxs = zeros(1, N_theta);
bf_sums = zeros(1, N_theta);
bf_sum_wraps = zeros(1, N_theta);
bf_sum_clips = zeros(1, N_theta);
d_wraps = zeros(1, N_theta);
d_clips = zeros(1, N_theta);
for I = 1:N_theta
% Calculate beamformer sum, using intermediate wrapping or clipping
phi = phis(I);
kphis = phi*ants;
% Sort kphis to have worst case order where the BF inputs first add
% all up constructively and then all down constructively, to have
% the largest intermediate BF sum results.
kphis = sort(mod(phi*ants, 2*pi));
s = 0;
s_wrap = 0;
s_clip = 0;
for k = 1:N_ant
a = round(A + R*cos(kphis(k)));
%a = round(A + R*sin(kphis(k)));
s = s + a;
if abs(s) > bf_sums_maxs(I)
bf_sums_maxs(I) = abs(s);
end
s_wrap = wrap(s_wrap + a, s_period);
s_clip = clip(s_clip + a, s_min, s_max);
end
bf_sums(I) = s;
bf_sum_wraps(I) = s_wrap;
bf_sum_clips(I) = s_clip;
% Derive difference with full range final beamformwer sum
d_final_wraps(I) = s - s_wrap;
d_final_clips(I) = s - s_clip;
% Derive difference with w bit output beamformwer sum
d_output_wraps(I) = wrap(s, s_period) - s_wrap;
d_output_clips(I) = clip(s, s_min, s_max) - s_clip;
end
%% Plots
xfig = 300;
yfig = 200;
xfigw = 1000;
yfigw = 800;
dfig = 20;
%% Plot final beamformer sums
fig=fig+1;
figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
figure(fig);
plot(thetas, bf_sums, 'k', thetas, bf_sum_wraps, 'r--', thetas, bf_sum_clips, 'b--');
title('Final beamformer sums');
xlabel('Beam angle [degrees]');
ylabel('Final sum');
grid on;
%% Plot abs maximum of intermediate beamformer sums
fig=fig+1;
figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
figure(fig);
plot(thetas, bf_sums_maxs, 'r');
title('Maximum of intermediate beamformer sums');
xlabel('Beam angle [degrees]');
ylabel('Intermediate maximum sum');
grid on;
%% Plot differences between final sum and wrapped sum, clipped sum
fig=fig+1;
figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
figure(fig);
plot(thetas, d_final_wraps, 'r', thetas, d_final_clips, 'b--');
title('Difference final beamformer sums');
xlabel('Beam angle [degrees]');
ylabel('Difference of final sum');
grid on;
%% Plot differences between output sum and wrapped sum, clipped sum
fig=fig+1;
figure('position', [xfig+fig*dfig yfig-fig*dfig xfigw yfigw]);
figure(fig);
plot(thetas, d_output_wraps, 'r', thetas, d_output_clips, 'b--');
title('Difference output beamformer sums');
xlabel('Beam angle [degrees]');
ylabel('Difference of final sum');
grid on;
matlab/bsn_source.m 0 → 100644
View file @ bf1f35cb
%-----------------------------------------------------------------------------
%
% Copyright (C) 2016
% 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: E. Kooistra, 2016
%
% Purpose : BSN source per block of data
% Description : Increment ctrl.bsn for every block.
function [state] = bsn_source(ctrl)
% Increment BSN for next call
state = ctrl;
state.bsn = ctrl.bsn + 1;
\ No newline at end of file
matlab/clip.m 0 → 100644
View file @ bf1f35cb
%-----------------------------------------------------------------------------
%
% Copyright (C) 2019
% 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: E. Kooistra, 2019
%
% Purpose : Clip input s into range [clip_min, clip_max]
function [s_clip] = clip(s, clip_min, clip_max)
s_clip = s;
if s > clip_max
s_clip = clip_max;
elseif s < clip_min
s_clip = clip_min;
end
end
matlab/clk_spectrum.m 0 → 100644
View file @ bf1f35cb
%-----------------------------------------------------------------------------
%
% Copyright (C) 2019
% 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