source_components

Tests for the functions in WODEN/src/source_components.cu. These functions calculate visibilities for the different COMPONENT types, as well as calling the various beam models to be applied to the visibilities.

test_apply_beam_gains.c

This calls source_components::test_kern_apply_beam_gains, which tests source_components::kern_apply_beam_gains. This kernel applies beam gain and leakage terms to Stokes visibilities to create linear Stokes polarisation visibilities via:

\[\begin{split}\begin{eqnarray} \mathrm{V}^{XX}_{12} = (g_{1x}g_{2x}^{\ast} + D_{1x}D_{2x}^{\ast})\mathrm{V}^{I}_{12} + (g_{1x}g_{2x}^{\ast} - D_{1x}D_{2x}^{\ast})\mathrm{V}^{Q}_{12} \\ + (g_{1x}D_{2x}^{\ast} + D_{1x}g_{2x}^{\ast})\mathrm{V}^{U}_{12} + i(g_{1x}D_{2x}^{\ast} - D_{1x}g_{2x}^{\ast})\mathrm{V}^{V}_{12} \end{eqnarray}\end{split}\]

\[\begin{split}\begin{eqnarray} \mathrm{V}^{XY}_{12} = (g_{1x}D_{2y}^{\ast} + D_{1x}g_{2y}^{\ast})\mathrm{V}^{I}_{12} + (g_{1x}D_{2y}^{\ast} - D_{1x}g_{2y}^{\ast})\mathrm{V}^{Q}_{12} \\ + (g_{1x}g_{2y}^{\ast} + D_{1x}D_{2y}^{\ast})\mathrm{V}^{U}_{12} + i(g_{1x}g_{2y}^{\ast} - D_{1x}D_{2y}^{\ast})\mathrm{V}^{V}_{12} \end{eqnarray}\end{split}\]

\[\begin{split}\begin{eqnarray} \mathrm{V}^{XY}_{12} = (g_{1x}D_{2y}^{\ast} + D_{1x}g_{2y}^{\ast})\mathrm{V}^{I}_{12} + (g_{1x}D_{2y}^{\ast} - D_{1x}g_{2y}^{\ast})\mathrm{V}^{Q}_{12} \\ + (g_{1x}g_{2y}^{\ast} + D_{1x}D_{2y}^{\ast})\mathrm{V}^{U}_{12} + i(g_{1x}g_{2y}^{\ast} - D_{1x}D_{2y}^{\ast})\mathrm{V}^{V}_{12} \end{eqnarray}\end{split}\]

\[\begin{split}\begin{eqnarray} \mathrm{V}^{YY}_{12} = (D_{1y}D_{2y}^{\ast} + g_{1y}g_{2y}^{\ast})\mathrm{V}^{I}_{12} + (D_{1y}D_{2y}^{\ast} - g_{1y}g_{2y}^{\ast})\mathrm{V}^{Q}_{12} \\ + (D_{1y}g_{2y}^{\ast} + g_{1y}D_{2y}^{\ast})\mathrm{V}^{U}_{12} + i(D_{1y}g_{2y}^{\ast} - g_{1y}D_{2y}^{\ast})\mathrm{V}^{V}_{12} \end{eqnarray}\end{split}\]

where \(\ast\) means complex conjugate, \(g_x, D_x, D_y, g_y\) are beam gain and leakage terms, with subscript 1 and 2 meaning antenna 1 and 2, and \(\mathrm{V}^I, \mathrm{V}^Q, \mathrm{V}^U, \mathrm{V}^V\) the Stokes visibilities. These tests try a number of combinations of values that have a simple outcome, and tests that the function returns the expected values. The combinations are shown in the table below. For all combinations, the beam gain and leakage is used for both antennas in the above equations. Each entry is a complex values and should be read as real,imag. Both FLOAT and DOUBLE code are just tested to a 32 bit accuracy here as these a simple numbers that require little precision.

\(g_x\)	\(D_x\)	\(D_y\)	\(g_y\)	\(\mathrm{V}^I\)	\(\mathrm{V}^Q\)	\(\mathrm{V}^U\)	\(\mathrm{V}^V\)	\(\mathrm{V}^{XX}\)	\(\mathrm{V}^{XY}\)	\(\mathrm{V}^{YX}\)	\(\mathrm{V}^{YY}\)
1,0	0,0	0,0	1,0	1,0	0,0	0,0	0,0	1,0	0,0	0,0	1,0
1,0	0,0	0,0	1,0	0,0	1,0	0,0	0,0	1,0	0,0	0,0	-1,0
1,0	0,0	0,0	1,0	0,0	0,0	1,0	0,0	0,0	1,0	1,0	0,0
1,0	0,0	0,0	1,0	0,0	0,0	0,0	1,0	0,0	0,1	0,-1	0,0
0,0	1,0	1,0	0,0	1,0	0,0	0,0	0,0	1,0	0,0	0,0	1,0
0,0	1,0	1,0	0,0	0,0	1,0	0,0	0,0	-1,0	0,0	0,0	1,0
0,0	1,0	1,0	0,0	0,0	0,0	1,0	0,0	0,0	1,0	1,0	0,0
0,0	1,0	1,0	0,0	0,0	0,0	0,0	1,0	0,0	0,-1	0,1	0,0
2,0	2,0	2,0	2,0	1,0	0,0	0,0	0,0	8,0	8,0	8,0	8,0
2,0	2,0	2,0	2,0	0,0	1,0	0,0	0,0	0,0	0,0	0,0	0,0
2,0	2,0	2,0	2,0	0,0	0,0	1,0	0,0	8,0	8,0	8,0	8,0
2,0	2,0	2,0	2,0	0,0	0,0	0,0	1,0	0,0	0,0	0,0	0,0
1,2	3,4	5,6	7,8	1,0	0,0	0,0	0,0	30,0	70,8	70,-8	174,0
1,2	3,4	5,6	7,8	0,0	1,0	0,0	0,0	-20,0	-36,0	-36,0	-52,0
1,2	3,4	5,6	7,8	0,0	0,0	1,0	0,0	22,0	62,8	62,-8	166,0
1,2	3,4	5,6	7,8	0,0	0,0	0,0	1,0	-4,0	-4,-16	-4,16	-4,0

test_calc_measurement_equation.c

This calls source_components::test_kern_calc_measurement_equation, which calls source_components::kern_calc_measurement_equation, which in turn is testing the device code source_components::calc_measurement_equation (which is used internally in WODEN in another kernel that calls multiple device functions, so I had to write a new kernel to test it alone.)

The following methodology is also written up the JOSS paper (TODO link JOSS paper once it’s accepted), where it’s used in a sky model + array layout through to visibility end-to-end version. Here we directly test the fuction above.

calc_measurement_equation calculates the phase-tracking measurement equation:

\[V(u,v,w) = \exp \left[ 2\pi i\left( ul + vm + w(n-1) \right) \right]\]

We can use Euler’s formula to split this into real and imaginary components. If I label the phase for a particular source and baseline as

\[\phi = 2\pi \left( ul + vm + w(n-1)\right)\]

then the real and imaginary parts of the visibility \(V_{re}\), \(V_{im}\) are

\[\begin{split}V_{re} = \cos(\phi) \\ V_{im} = \sin(\phi)\end{split}\]

A definitive test of the calc_measurement_equation function then is to then set \(\phi\) to a number of values which produce known sine and cosine outputs, by selecting specific combinations of u,v,w and l,m,n. First of all, consider the case when u,v,w = 1,1,1. In that case,

\[\frac{\phi_{\mathrm{simple}}}{2\pi} = l + m + (n-1).\]

if we further set l == m, we end up with

\[\begin{split}\frac{\phi_{\mathrm{simple}}}{2\pi} = 2l + (n-1), \\ l = \sqrt{\left( \frac{1 - n^2}{2} \right)}\end{split}\]

I shoved those two equations into Wolfram Alpha who assured me that a solution for n here is

\[n = \frac{\sqrt{2}\sqrt{-\phi_{\mathrm{simple}}^2 - 4\pi\phi_{\mathrm{simple}} + 8\pi^2} + \phi_{\mathrm{simple}} + 2\pi}{6\pi}\]

which we can then use to calculate values for l,m through

\[l = m = \sqrt{\frac{1 - n^2}{2}}.\]

By selecting the following values for \(\phi\), we can create the following set of l,m,n coords, which have the a known set of outcomes:

\(\phi_{\mathrm{simple}}\)	l,m	n	\(\cos(\phi)\)	\(\sin(\phi)\)
\(0\)	0.0	1.0	\(1.0\)	\(0\)
\(\pi/6\)	0.0425737516338956	0.9981858300655398	\(\sqrt{3}/2\)	\(0.5\)
\(\pi/4\)	0.0645903244635131	0.9958193510729726	\(\sqrt{2}/2\)	\(\sqrt{2}/2\)
\(\pi/3\)	0.0871449863555500	0.9923766939555675	\(0.5\)	\(\sqrt{3}/2\)
\(\pi/2\)	0.1340695840364469	0.9818608319271057	\(0.0\)	\(1.0\)
\(2\pi/3\)	0.1838657911209207	0.9656017510914922	\(-0.5\)	\(\sqrt{3}/2\)
\(3\pi/4\)	0.2100755148372292	0.9548489703255412	\(-\sqrt{2}/2\)	\(\sqrt{2}/2\)
\(5\pi/6\)	0.2373397982598921	0.9419870701468823	\(-\sqrt{3}/2\)	\(0.5\)
\(\pi\)	0.2958758547680685	0.9082482904638630	\(-1.0\)	\(0.0\)
\(7\pi/6\)	0.3622725654470420	0.8587882024392495	\(-\sqrt{3}/2\)	\(-0.5\)
\(5\pi/4\)	0.4003681253515569	0.8242637492968862	\(-\sqrt{2}/2\)	\(-\sqrt{2}/2\)

Note

If you try and go higher in \(\phi\) then because I set \(l == m\) you no longer honour \(\sqrt{l^2 + m^2 + n^2} <= 1.0\) I think this range of angles is good enough coverage though.

This is a great test for when \(u,v,w = 1\), but we want to test a range of baseline lengths to check our function is consistent for short and long baselines. We can play another trick, and set all baseline coords to be equal, i.e. \(u = v = w = b\) where \(b\) is baseline length. In this form, the phase including the baseline length \(\phi_{b}\) is

\[\phi_{b} = 2\pi b\left( l + m + n - 1 \right) = b\phi_{\mathrm{simple}}.\]

As sine/cosine are periodic functions, the following is true:

\[\phi_{\mathrm{simple}} = \phi_{\mathrm{simple}} + 2\pi \mathrm{n}\]

where \(\mathrm{n}\) is some integer. This means for a given \(\phi_{\mathrm{simple}}\) from the table above, we can find an appropriate \(b\) that should still result in the expected sine and cosine outputs by setting

\[\begin{split}b\phi_{\mathrm{simple}} = \phi_{\mathrm{simple}} + 2\pi \mathrm{n} \\ b = \frac{\phi_{\mathrm{simple}} + 2\pi \mathrm{n}}{\phi_{\mathrm{simple}}}\end{split}\]

for a range of \(\mathrm{n}\) values. The values of \(\mathrm{n}\) and the resultant size of \(b\) that I use in testing are shown in the table below (note for \(\phi_{\mathrm{simple}} = 0\) I just set \(b = 2\pi \mathrm{n}\) as the effects of \(l,m,n\) should set everything to zero regardless of baseline coords).

\(\phi_{\mathrm{simple}}\)	\(b(\mathrm{n=0})\)	\(b(\mathrm{n=1})\)	\(b(\mathrm{n=10})\)	\(b(\mathrm{n=100})\)	\(b(\mathrm{n=1000})\)	\(b(\mathrm{n=10000})\)
\(0\)	0.0	6.3	62.8	628.3	6283.2	62831.9
\(\pi/6\)	1.0	13.0	121.0	1201.0	12001.0	120001.0
\(\pi/4\)	1.0	9.0	81.0	801.0	8001.0	80001.0
\(\pi/3\)	1.0	7.0	61.0	601.0	6001.0	60001.0
\(\pi/2\)	1.0	5.0	41.0	401.0	4001.0	40001.0
\(2\pi/3\)	1.0	4.0	31.0	301.0	3001.0	30001.0
\(3\pi/4\)	1.0	3.7	27.7	267.7	2667.7	26667.7
\(5\pi/6\)	1.0	3.4	25.0	241.0	2401.0	24001.0
\(\pi\)	1.0	3.0	21.0	201.0	2001.0	20001.0
\(7\pi/6\)	1.0	2.7	18.1	172.4	1715.3	17143.9
\(5\pi/4\)	1.0	2.6	17.0	161.0	1601.0	16001.0

This gives a range of baseline lengths from 1 to \(> 10^4\) wavelengths.

In this test, I run every combination of \(l,m,n\) and \(u,v,w = b\) for each \(\phi_{\mathrm{simple}}\) from the tables above, and assert that the real and imaginary of every output visibility match the expected values of \(\sin(\phi_{\mathrm{simple}})\) and \(\cos(\phi_{\mathrm{simple}})\). When compiling with FLOAT precision, I assert the outputs must be within an absolute tolerance of 2e-3, and for DOUBLE a tolerance of 2e-9.

The error scales with the length of baseline, as shown in this plot below. Here, I have plotted the fractional offset of the recovered value of \(\sin(\phi_{\mathrm{simple}})\) (imaginary part of the visibility) and \(\cos(\phi_{\mathrm{simple}})\) (real part of the visibility), compared to their analytically expected outcome. I’ve plotted each \(\phi_{\mathrm{simple}}\) as a different symbol, with the FLOAT in blue and DOUBLE in orange.

You can see as you increase baseline length, a general trend of increasing error is seen. Note these are the absolute differences plotted here, to work on a log10 scale. For the DOUBLE results, it’s close to 50% a negative or positive offset from expected. The FLOAT results that perform the worst (where \(\phi_{\mathrm{simple}} = 3\pi/4,\, 7\pi/6\)) correspond to values of \(b\) with large fractional values (see the table above), showing how the 32 bit precision fails to truthfully report large fractional numbers.

As a second test, I setup 10,000 u,v,w ranging from -1000 to 1000 wavelengths, and 3600 l,m,n coordinates that span the entire sky, and run them through source_components::test_kern_calc_measurement_equation and check they equal the equivalent measurement equation as calculated by C in 64 bit precisions. This checks the kernel works for a range of input coordinates. I assert the CUDA outputs must be within an absolute tolerance of 1e-7 for the FLOAT code, and 1e-15 for the DOUBLE code.

test_extrap_stokes.c

This calls source_components::test_extrap_stokes_all_models, which calls source_components::extrapolate_Stokes, which handles extrapolating a the flux density of component to a given set of frequencies. This covers all types of flux behaviours (currently POWER_LAW, CURVED_POWER_LAW, and LIST types).

Six components of each type of flux behaviour are tested, each with randomly assigned values (these values are generated using the script WODEN/cmake_testing/source_components/write_test_extrap_stokes_header.py). For the LIST flux types, each list has a random number of entries, as well as each list entry begin given a random flux. The values are copied into a source_t struct, passed through the CUDA code, and extrapolated to 25 frequencies between 50 and 300 MHz. The outputs are tested against equivalent C functions in double precision. The woden_double code is tested to an absolute precision of 1e-11 Jy, with the woden_float a 9e-4 Jy precision (note some of the extrapolated fluxes are of order 1e3).

To visualise the results, run WODEN/cmake_testing/source_components/test_extrap_stokes.py, which will produce the following plots. The plots only show the Stokes I outputs, but all Stokes values are tested against the C test code.

For the POWER_LAW type components:

For the CURVED_POWER_LAW type components (where I have set the curve of the peak to be within 100 to 200 MHz so they should all curve in these plots):

Finally, for the LIST type components. Note here, the black line is the information contained in the sky model, the little orange crosses are a python implementation of the linear interpolation between points to double check everything, and the cyan squares are what is output by the CUDA code.

Note

Most SEDs are assumed to be close to a power-law, and so are linear in log space. The linear interpolation is therefore done in log space. However, for some complicated modelling, negative flux values are required. For these values, the interpolation is done in linear space. Hence in the plot below, where the fluxes dip into negative territory, the extrapolated fluxes not longer lie on the straight lines, as these are log-log plots.

test_get_beam_gains.c

This calls source_components::test_kern_get_beam_gains, which calls source_components::kern_get_beam_gains, which in turn is testing the device code source_components::get_beam_gains. This function handles grabbing the pre-calculated beam gains for a specific beam model, time, and frequency (assuming the beam gains have already been calculated). kern_get_beam_gains is setup to call get_beam_gains for multiple inputs and recover them into a set of output arrays.

Beam gain calculations are stored in primay_beam_J* arrays, including all frequency and time steps, as well as all directions on the sky. This test sets all real entries in the four primay_beam_J* beam gain arrays to the value of their index. In this way, we can easily predict the expected value in the outputs as being the index of the beam gain we wanted to select. I’ve set the imaginary to zero.

This test runs with two time steps, two frequency steps, three baselines, and four beam models. Three different outcomes are expected given the beam model:

• ANALY_DIPOLE, GAUSS_BEAM: The values of the gains are tested to match the expected index. The leakage terms are tested to be zero as the models have no leakage terms

• FEE_BEAM, FEE_BEAM_INTERP, MWA_ANALY: Both the beam gain and leakage terms are tested as these models include leakage terms

• NO_BEAM: The gain terms are tested to be 1.0, and leakage to be 0.0

Both FLOAT and DOUBLE code are tested to a 32 bit accuracy here as these are simple numbers that require little precision.

test_source_component_common.c

This calls source_components::test_source_component_common, which calls source_components::source_component_common. source_component_common is run by all visibility calculation functions (the functions kern_calc_visi_point, kern_calc_visi_gauss, kern_calc_visi_shape). It handles calculating the l,m,n coordinates and beam response for all COMPONENTs in a sky model, regardless of the type of COMPONENT.

Similarly to the tests in test_lmn_coords.c, I setup a slice of 9 RA coordinates, and hold the Dec constant, set the phase centre to RA\(_{\textrm{phase}}\), Dec\(_{\textrm{phase}}\) = \(0^\circ, 0^\circ\). This way I can analytically predict what the l,m,n calculated coordinates should be (which are tested to be within 1e-15 of expected values).

In these tests I run with three time steps, two frequency steps (100 and 200 MHz), and five baselines (the coordinates of which don’t matter, but change the size of the outputs, so good to have a non-one value). I input a set of az,za coords that match the RA,Dec coords for an LST = 0.0. As I have other tests that check the sky moves with time, I just set the sky to be stationary with time here, to keep the test clean.

For each primary beam type, I run the 9 COMPONENTs through the test, and check the calcualted l,m,n are correct, and check that the calculated beam values match a set of expected values, which are stored in test_source_component_common.h. As with previous tests varying the primary beam, I check that leakage terms should be zero when the model doesn’t include them.

The absolute tolerance values used for the different beam models, for the two different precisions are shown in the table below. Note I’ve only stored the expected values for some of the beams to 1e-7 accuracy, as the absolute accuracy of these functions beam functions is tested elsewhere. The GAUSS_BEAM values are calculated analytically in the same test as described in test_gaussian_beam.c.

Beam type	FLOAT tolerance	DOUBLE tolerance
GAUSS_BEAM	1e-7	1e-12
ANALY_DIPOLE	1e-6	1e-7
FEE_BEAM	1e-7	1e-7
FEE_BEAM_INTERP	1e-7	1e-7
MWA_ANALY	1e-7	1e-12

Of the nine components tested, 3 are given POWER_LAW, 3 are CURVED_POWER_LAW, and 3 are LIST flux styles. Similarly to test_extrap_stokes.c, each component is given a random set of parameters, and compared to a C version of the functions to extrapolate the Stokes parameters. The woden_double code is tested to an absolute precision of 1e-12 Jy, with the woden_float a 1e-4 Jy precision (note some of the extrapolated fluxes are of order 1e3).

test_kern_calc_visi_point.c

This calls source_components::test_kern_calc_visi_all, which in turn calls source_components::kern_calc_visi_point_or_gauss (in this case being used for POINT components). This kernel calculates the visibility response for POINT COMPONENTs for a number of sky directions, for all time and frequency steps, and all baselines.

I set up a grid of 25 l,m coords with l,m coords ranging over -0.5, -0.25, 0.0, 0.25, 0.5. I run a simulation with 10 baselines, where I set u,v,w to:

\[\begin{split}u,v = 100(b_{\mathrm{ind}} + 1) \\ w = 10(b_{\mathrm{ind}} + 1)\end{split}\]

where \(b_{\mathrm{ind}}\) is the baseline index, meaning the test covers the baseline range \(100 < u,v <= 1000\) and \(10 < w <= 100\). The test also runs three frequncies, 150, 175, 200 MHz, and two time steps. As I am providing predefined u,v,w, I don’t need to worry about LST effects, but I simulate with two time steps to make sure the resultant visibilities end up in the right order.

Overall, I run three groups of tests here:

• Keeping the beam gains and flux densities constant at 1.0

• Varying the flux densities with frequency and keeping the beam gains constant at 1.0. I set all components to just be Stokes I, as the other polarisations are tested elsewhere. There are 10 POWER_LAW, 10 CURVED_POWER_LAW, and 5 LIST type components, each with randomly generated values. Equivalent C code is used to calculate the expected fluxes

• Varying the beam gains with frequency and keeping the flux densities constant at 1.0. As the beam can vary with time, frequency, and direction on sky, I assign each beam gain a different value. As num_freqs*num_times*num_components = 375, I set the real of all beam gains to \(\frac{1}{375}(B_{\mathrm{ind}} + 1)\), where \(B_{\mathrm{ind}}\) is the beam value index. This way we get a unique value between 0 and 1 for all beam gains, allowing us to test time/frequency is handled correctly by the function under test

Each set of tests is run for all six primary beam types, so a total of 18 tests are called. Each test calls kern_calc_visi_point, which should calculate the measurement equation for all baselines, time steps, frequency steps, and COMPONENTs. It should also sum over COMPONENTs to get the resultant visibility for each baseline, time, and freq. To test the outputs, I have created equivalent C functions at 64 bit precision in test_kern_calc_visi_common.c to calculate the measurement equation for the given inputs. For all visibilities, for the FLOAT version I assert the CUDA code output must match the C code output to within an fractional tolerance of 7e-5 to the C value, for both the real and imaginary parts. For the DOUBLE code, the fractional tolerance is 1e-13. I’ve switched to fractional tolerance here as the range of magnitudes covered by these visibilities means a small absolute tolernace will fail a large magnitude visibility when it reports a value that is correct to 1e-11%.

test_kern_calc_visi_gauss.c

This calls source_components::test_kern_calc_visi_all, which in turn calls source_components::kern_calc_visi_point_or_gauss (in this case being used for GAUSSIAN components). This kernel calculates the visibility response for GAUSSIAN COMPONENTs for a number of sky directions, for all time and frequency steps, and all baselines.

This runs all tests as described by test_kern_calc_visi_point.c, plus a fourth set of tests that varies the position angle, major, and minor axis of the input GAUSSIAN components, for a total of 24 tests. Again, I have C code to test the CUDA code against. I assert the CUDA code output must match the C code output to within an fractional tolerance of 7e-5 to the C value, for both the real and imaginary parts. For the DOUBLE code, the fractional tolerance is 1e-12.

test_kern_calc_visi_shape.c

This calls source_components::test_kern_calc_visi_all, which calls source_components::kern_calc_visi_shape. This kernel calculates the visibility response for SHAPELET COMPONENTs for a number of sky directions, for all time and frequency steps, and all baselines.

This runs all tests as described by test_kern_calc_visi_gauss.c, plus a fifth set of tests that gives multiple shapelet basis function parameters to the input SHAPELET components, for a total of 30 tests. Again, I have C code to test the CUDA code against.

The final 5th test really pushes the FLOAT code hard, as the range of magnitudes of the visibilities is large. As a result, FLOAT code is tested to to within a fractional tolerance of 5e-3 to the C values (which happens mostly when the expected value is around 1e-5 Jy, so a fractional offset of 5e-3 is an absolute offset of 5e-8 Jy), for both the real and imaginary parts. For the DOUBLE code, the fractional tolerance is 1e-12.

test_update_sum_visis.c

This calls source_components::test_kern_update_sum_visis, which in turn calls source_components::kern_update_sum_visis. This kernel gathers pre-calculated primary beam values, unpolarised visibilities, Stokes flux densities, and combines them all into linear polarisation Stokes visibilities. It then sums all COMPONENTs together onto the final set of linear polarisation visibilities.

This code runs three sets of tests, each with three baselines, four time steps, three frequencies, and ten COMPONENTs. For each set of tests, the input arrays that are being tested have their values set to the index of the value, making the summation able to test whether the correct time, frequency, and beam indexes are being summed into the resultant visibilities. The three sets of tests consist of:

• Varying the beam gains, while keeping the flux densities and unpolarised measurement equations constant.

• Varying the flux densities, while keeping the beam gains and unpolarised measurement equations constant.

• Varying the unpolarised measurement equations, while keeping the beam gains and flux densities constant.

Each set of tests is run for all primary beam types, for a total of 18 tests. The different beam models have different expected values depending on whether they include leakage terms or not.

test_kern_calc_autos.c

This calls source_components::test_kern_calc_autos, which in turn calls source_components::kern_calc_autos. This kernel gathers pre-calculated primary beam values, and Stokes flux densities, and performs a dot-product to calculate the auto-correlations of all antennas, for all time steps.

This test runs by creating a Stokes I sky model of 4 components, extrapolated to three frequencies and two time steps. Complex beam gains are generated for all directions, times, and frequencies. The auto-correlation calculation is then performed by C code, which the outputs of the GPU code are tested against.

As some primary beam models are real only, and only some have leakage, the function is tested for all beam model types separately, with the expected outcomes calculated accordingly.

The woden_double code is tested to an absolute precision of 1e-12 Jy, with the woden_float a 1e-2 Jy precision as compared to the C code (note some of the resultant fluxes are >10,000 Jy, hence the large absolute error. You should probably just always use the double precision version).