primary_beam_cuda

API documentation for primary_beam_cuda.cu.

Device methods to calculate perfect Gaussian and Analytic Dipole primary beam responses. Currently, the analytic dipole is fixed to being an MWA dipole. Both models assume there is no leakage and beams are purely real.

Functions

__device__ void twoD_Gaussian(user_precision_t x, user_precision_t y, user_precision_t xo, user_precision_t yo, user_precision_t sigma_x, user_precision_t sigma_y, user_precision_t cos_theta, user_precision_t sin_theta, user_precision_t sin_2theta, user_precision_t *d_beam_real, user_precision_t *d_beam_imag)

Calculate a two dimensional Gaussian.

Returns the Gaussian as defined on Wikipedia here via the equation:

\[ G(x,y) = \exp \left( -\left( a(x-x_o)^2 + 2b(x-x_o)(y-y_o) + c(y-y_o)^2 \right) \right) \]

where

\[\begin{split}\begin{eqnarray*} a &=& \frac{\cos(\theta)^2}{2\sigma_x^2} + \frac{\sin(\theta)^2}{2\sigma_y^2} \\ b &=& -\frac{\sin(2\theta)}{4\sigma_x^2} + \frac{\sin(2\theta)}{4\sigma_y^2} \\ c &=& \frac{\sin(\theta)^2}{2\sigma_x^2} + \frac{\cos(\theta)^2}{2\sigma_y^2} \end{eqnarray*}\end{split}\]

with \((x_0,y_0)\) the central coordinates, \((\sigma_x,\sigma_y)\) the standard deviations.

Parameters:
  • x[in] x coordinate

  • y[in] y coordinate

  • xo[in] Central x coordinate

  • yo[in] Central y coordinate

  • sigma_x[in] Square root of variance in x

  • sigma_y[in] Square root of variance in x

  • cos_theta[in] Cosine of the rotation angle

  • sin_theta[in] Sine of the rotation angle

  • sin_2theta[in] Sine of two times the rotation angle

  • *d_beam_real[inout] Real part of the Gaussian repsonse

  • *d_beam_imag[inout] Imaginary part of the Gaussian reponse

__global__ void kern_gaussian_beam(double *d_beam_ls, double *d_beam_ms, double beam_ref_freq, double *d_freqs, user_precision_t fwhm_lm, user_precision_t cos_theta, user_precision_t sin_theta, user_precision_t sin_2theta, int num_freqs, int num_times, int num_components, cuUserComplex *d_primay_beam_J00, cuUserComplex *d_primay_beam_J11)

Kernel to calculate a Gaussian primary beam response at the given interferometric sky \((l,m)\) coords and frequency.

The primary beam is to be calculated for each sky direction, each time step, and each frequency. The size of the beam repsonse on the sky changes with frequency, so need a reference frequency beam_ref_freq and full-width half-maximum (fwhm_lm, in \(l,m\) units) to scale the beam width with frequency. The Gaussian beam is held to point at a given az/za, so the sky positions to calculate in d_beam_ls, d_beam_ms should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The outputs are stored in d_primay_beam_J00, d_primay_beam_J11, where 00 refers to the north-south polarisation, 11 the east-west polarisation, in order of time, frequency, COMPONENT.

When called with dim3 grid, threads, kernel should be called with both grid.x and grid.y defined, where:

  • grid.x * threads.x >= num_components * num_time_steps

  • grid.y * threads.y >= num_freqs

Parameters:
  • *d_beam_ls[in] Array of \(l\) coords to calculate beam at

  • *d_beam_ms[in] Array pf \(m\) coords to calculate beam at

  • beam_ref_freq[in] Reference frequency at which the FWHM is applicable (Hz)

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • fwhm_lm[in] FWHM of the beam in \(l,m\) coords

  • cos_theta[in] Cosine of the rotation angle

  • sin_theta[in] Sine of the rotation angle

  • sin_2theta[in] Sine of two times the rotation angle

  • num_freqs[in] Number of frequencies being calculated

  • num_times[in] Number of time steps being calculated

  • num_components[in] Number of COMPONENTS the beam is calculated for

  • *d_primay_beam_J00[inout] Device array to store the beam Jones complex J[0,0] response in

  • *d_primay_beam_J11[inout] Device array to store the beam Jones complex J[1,1] response in

void calculate_gaussian_beam(int num_components, int num_time_steps, int num_freqs, user_precision_t ha0, user_precision_t sdec0, user_precision_t cdec0, user_precision_t fwhm_lm, user_precision_t cos_theta, user_precision_t sin_theta, user_precision_t sin_2theta, double beam_ref_freq, double *d_freqs, double *beam_has, double *beam_decs, cuUserComplex *d_primay_beam_J00, cuUserComplex *d_primay_beam_J11)

Calculate the Gaussian primary beam response at the given hour angle and declinations beam_point_has, beam_point_decs. Note the XX and YY repsonses are equal in this toy example.

The primary beam is to be calculated for each sky direction, each time step, and each frequency. The size of the beam repsonse on the sky changes with frequency, so need a reference frequency beam_ref_freq and full-width half-maximum (fwhm_lm, in \(l,m\) units) to scale the beam width with frequency. The Gaussian beam is held to point at a given az/za, so the sky positions to calculate in d_beam_ls, d_beam_ms should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The outputs are stored in d_primay_beam_J00, d_primay_beam_J11, where 00 refers to the north-south polarisation, 11 the east-west polarisation, in order of time, frequency, COMPONENT. This function uses the beam centre pointing ha0, dec0 to calculate an \(l,m\) coord system in which to calculate the Gaussian beam, using kern_gaussian_beam.

Parameters:
  • num_components[in] Number of COMPONENTS the beam is calculated for

  • num_time_steps[in] Number of time steps being calculated

  • num_freqs[in] Number of frequencies being calculated

  • ha0[in] Hour angle of beam pointing centre (radians)

  • sdec0[in] Sine of Declination of the beam pointing centre (radians)

  • cdec0[in] Cosine of Declination of the beam pointing centre (radians)

  • fwhm_lm[in] FWHM of the beam in \(l,m\) coords

  • cos_theta[in] Cosine of the rotation angle

  • sin_theta[in] Sine of the rotation angle

  • sin_2theta[in] Sine of two times the rotation angle

  • beam_ref_freq[in] Reference frequency at which the FWHM is applicable (Hz)

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • *beam_has[in] Array of Hour Angles to calculate the beam toward

  • *beam_decs[in] Array of Declinations to calculate the beam toward

  • *d_primay_beam_J00[inout] Device array to store the beam Jones complex J[0,0] response in

  • *d_primay_beam_J11[inout] Device array to store the beam Jones complex J[1,1] response in

__device__ void analytic_dipole(user_precision_t az, user_precision_t za, user_precision_t wavelength, cuUserComplex *d_beam_X, cuUserComplex *d_beam_Y)

Calculate the beam response of a north-south (X) and east-west (Y) analytic dipole on an infinite ground screen, for the given sky direciton az,za and wavelength.

Dipoles are assumed to be MWA, and given a length of 0.3 metres. Beam size on the sky scales with frequency hence the need for wavelength

Parameters:
  • az[in] Azimuth (radians)

  • za[in] Zenith Angle (radians)

  • wavelength[in] Wavelength (metres)

  • d_beam_X[inout] Complex beam value for north-south dipole

  • d_beam_Y[inout] Complex beam value for east-west dipole

__global__ void kern_analytic_dipole_beam(user_precision_t *d_azs, user_precision_t *d_zas, double *d_freqs, int num_freqs, int num_times, int num_components, cuUserComplex *d_primay_beam_J00, cuUserComplex *d_primay_beam_J11)

Kernel to calculate an Analytic MWA Dipole over an infinite ground screen at the given Azimuth and Zenith Angles d_azs, d_zas and frequencies d_freqs.

The primary beam is to be calculated for each sky direction, each time step, and each frequency. The Analytic dipole beam is stationary on the sky, so the Azimuth and Zenith Angles in azs,zas should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The outputs are stored in d_primay_beam_J00, d_primay_beam_J11, where 00 refers to the north-south polarisation, 11 the east-west polarisation, in order of time, frequency, COMPONENT. Beam outputs are normalised to zenith

When called with dim3 grid, threads, kernel should be called with both grid.x and grid.y defined, where:

  • grid.x * threads.x >= num_components * num_time_steps

  • grid.y * threads.y >= num_freqs

Todo:

Make the zenith normalisation an option

Parameters:
  • *d_azs[in] Array of Azimuth angles to calculate the beam towards (radians)

  • *d_zas[in] Array of Zenith Angles to calculate the beam towards (radians)

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • num_freqs[in] Number of frequencies being calculated

  • num_times[in] Number of time steps being calculated

  • num_components[in] Number of COMPONENTS the beam is calculated for

  • *d_primay_beam_J00[inout] Device array to store the beam Jones complex J[0,0] response in

  • *d_primay_beam_J11[inout] Device array to store the beam Jones complex J[1,1] response in

void calculate_analytic_dipole_beam(int num_components, int num_time_steps, int num_freqs, user_precision_t *azs, user_precision_t *zas, double *d_freqs, cuUserComplex *d_primay_beam_J00, cuUserComplex *d_primay_beam_J11)

Calculate the Analytic Dipole over an infinite ground screen primary beam response at the given Azimuth and Zenith Angles azs,zas, and frequencies d_freqs.

The primary beam is to be calculated for each sky direction, each time step, and each frequency. The Analytic dipole beam is stationary on the sky, so the Azimuth and Zenith Angles in azs,zas should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The outputs are stored in d_primay_beam_J00, d_primay_beam_J11, where 00 refers to the north-south polarisation, 11 the east-west polarisation, in order of time, frequency, COMPONENT. Beam outputs are normalised to zenith. Note eveything starting with d_ should be in device memory.

When called with dim3 grid, threads, kernel should be called with both grid.x and grid.y defined, where:

  • grid.x * threads.x >= num_components * num_time_steps

  • grid.y * threads.y >= num_freqs

Parameters:
  • num_components[in] Number of COMPONENTS the beam is calculated for

  • num_time_steps[in] Number of time steps being calculated

  • num_freqs[in] Number of frequencies being calculated

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • *azs[in] Array of Azimuth angles to calculate the beam towards (radians)

  • *zas[in] Array of Zenith Angles to calculate the beam towards (radians)

  • *d_primay_beam_J00[inout] Device array to store the beam Jones complex J[0,0] response in

  • *d_primay_beam_J11[inout] Device array to store the beam Jones complex J[1,1] response in

__device__ void RTS_MWA_beam(user_precision_t az, user_precision_t za, double ha, double dec, double wavelength, double *d_metre_delays, double latitude, int norm, cuUserComplex *gx, cuUserComplex *Dx, cuUserComplex *Dy, cuUserComplex *gy)

Calculate the analytic MWA beam response to a single az/za direction for a give wavelength and set of delays (d_metre_delays)

Based on the RTS code found in mwa_tile::local_FillMatrices. I think this code inherently does some kind of parallatic rotation, hence needs the ha/dec along with the az/za. The delays added to the paths of individual dipoles allow the beam to be ‘pointed’. The physical length of these paths should be given in d_metre_delays (this conversion from the delays given in the MWA metafits is handled by primary_beam_cuda::calculate_RTS_MWA_analytic_beam)

Parameters:
  • az[in] Azimuth to calculate the beam toward (radians)

  • za[in] Zenith angle to calculate the beam toward (radians)

  • ha[in] Corresponding hour angle to the given az, za

  • dec[in] Corresponding Declination to the given az, za

  • wavelength[in] Wavelength to calculate the beam at (metres)

  • d_metre_delays[in] The physical delays lengths added to each dipole to steer the beam (metres)

  • latitude[in] Latitude of the array (radians)

  • norm[in] Whether to normalise to zenith or not (1 to normalise, 0 to not)

  • *gx[inout] The gain for the north-south beam

  • *Dx[inout] The leakage for the north-south beam

  • *Dy[inout] The leakage for the east-west beam

  • *gy[inout] The gain for the east-west beam

__global__ void kern_RTS_analytic_MWA_beam(user_precision_t *d_azs, user_precision_t *d_zas, user_precision_t *d_beam_has, user_precision_t *d_beam_decs, double *d_metre_delays, double *d_freqs, double latitude, int norm, int num_freqs, int num_times, int num_components, cuUserComplex *d_gxs, cuUserComplex *d_Dxs, cuUserComplex *d_Dys, cuUserComplex *d_gys)

Kernel to calculate the analytic MWA primary beam to a set of sky directions d_azs and d_zas for a given set of delays d_metre_delays and frequencies d_freqs.

Kernel calls primary_beam_cuda::RTS_MWA_beam. The MWA primary beam is stationary on the sky for a given set of delays, so the Azimuth and Zenith Angles in azs,zas should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The delays added to the paths of individual dipoles allow the beam to be ‘pointed’. The physical length of these paths should be given in d_metre_delays (this conversion from the delays given in the MWA metafits is handled by primary_beam_cuda::calculate_RTS_MWA_analytic_beam)

When called with dim3 grid, threads, kernel should be called with both grid.x and grid.y defined, where:

  • grid.x * threads.x >= num_components * num_time_steps

  • grid.y * threads.y >= num_freqs

Parameters:
  • *d_azs[in] Array of Azimuth angles to calculate the beam towards (radians)

  • *d_zas[in] Array of Zenith Angles to calculate the beam towards (radiany)

  • *d_beam_has[in] Hour angles corresponding to azs,zas (radians)

  • *d_beam_decs[in] Declinations corresponding to azs,zas (radians)

  • *d_metre_delays[in] Delays that specify the pointing (path length, metres)

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • latitude[in] Latitude of the array (radians)

  • norm[in] Whether to normalise to zenith or not (1 to normalise, 0 to not)

  • num_freqs[in] Number of frequencies being calculated

  • num_times[in] Number of time steps being calculated

  • num_components[in] Number of COMPONENTS the beam is calculated for

  • *d_gxs[inout] The gains for the north-south beam

  • *d_Dxs[inout] The leakages for the north-south beam

  • *d_Dys[inout] The leakages for the east-west beam

  • *d_gys[inout] The gains for the east-west beam

void calculate_RTS_MWA_analytic_beam(int num_components, int num_time_steps, int num_freqs, user_precision_t *azs, user_precision_t *zas, user_precision_t *delays, double latitude, int norm, double *beam_has, double *beam_decs, double *d_freqs, cuUserComplex *d_gxs, cuUserComplex *d_Dxs, cuUserComplex *d_Dys, cuUserComplex *d_gys)

Calculate the analytic MWA primary beam to a set of sky directions azs and zas for a given set of delays delays as listed in an MWA metafits file, and frequencies d_freqs.

Uses the kernel primary_beam_cuda::kern_RTS_analytic_MWA_beam. The MWA primary beam is stationary on the sky for a given set of delays, so the Azimuth and Zenith Angles in azs,zas should contain num_components*num_times values, as the COMPONENTs move through the beam with time. The delays added to the paths of individual dipoles allow the beam to be ‘pointed’. The delays as listed in the metafits (given as delays) are listed in units of the delay time added internally to the tile (in seconds). This function coverts them into a path length (metres), as needed by primary_beam_cuda::RTS_MWA_beam`.

Parameters:
  • num_components[in] Number of COMPONENTS the beam is calculated for

  • num_time_steps[in] Number of time steps being calculated

  • num_freqs[in] Number of frequencies being calculated

  • *azs[in] Array of Azimuth angles to calculate the beam towards (radians)

  • *zas[in] Array of Zenith Angles to calculate the beam towards (radiany)

  • *delays[in] Delays that specificy the pointing (as reported in the MWA metafits)

  • latitude[in] Latitude of the array (radians)

  • norm[in] Whether to normalise to zenith or not (1 to normalise, 0 to not)

  • *beam_has[in] Hour angles corresponding to azs,zas (radians)

  • *beam_decs[in] Declinations corresponding to azs,zas (radians)

  • *d_freqs[in] Array of frequencies to calculate beam at (Hz)

  • *d_gxs[inout] The gains for the north-south beam

  • *d_Dxs[inout] The leakages for the north-south beam

  • *d_Dys[inout] The leakages for the east-west beam

  • *d_gys[inout] The gains for the east-west beam

void run_hyperbeam_cuda(int num_components, int num_time_steps, int num_freqs, uint8_t parallactic, struct FEEBeamGpu *cuda_fee_beam, double *azs, double *zas, double *latitudes, cuUserComplex *d_primay_beam_J00, cuUserComplex *d_primay_beam_J01, cuUserComplex *d_primay_beam_J10, cuUserComplex *d_primay_beam_J11)

Calculate the FEE MWA primary beam model to a set of sky directions azs and zas for a given initialised mwa_hyperbeam device beam object *cuda_fee_beam. NOTE that the azs, zas need to increment by component (fastest changing), then time (slowest changing). This is the OPPOSITE of what happens for all other functions, but is needed for pointer arithmatic that must be done to feed things into mwa_hyperbeam efficiently. Soz boz.

Calls mwa_hyperbeam::calc_jones_cuda_device to calculate the beam responses on the device. This function requires an initialised struct FEEBeamGpu *cuda_fee_beam object (initialised using mwa_hyperbeam::new_cuda_fee_beam), which in turn needs a struct FEEBeam *fee_beam (initialised using mwa_hyperbeam::new_fee_beam). Running these functions gathers the spherical harmnoic coefficients for the requested frequencies, as well as the delays to point the beam. If these aren’t setup correctly, this will fall flat on it’s face.

Once the beam repsonses have been calculated, split them up into the WODEN d_primay_beam_J* arrays using the kernel primary_beam_cuda::kern_map_hyperbeam_gains.

Parameters:
  • num_components[in] Number of COMPONENTS the beam is calculated for

  • num_time_steps[in] Number of time steps being calculated

  • num_freqs[in] Number of frequencies being calculated

  • parallactic[in] Whether to rotate by parallactic angle or not

  • *cuda_fee_beam[in] An initialised mwa_hyperbeam struct FEEBeamGpu

  • *azs[in] Array of Azimuth angles to calculate the beam towards (radians)

  • *zas[in] Array of Zenith Angles to calculate the beam towards (radians)

  • *latitudes[in] The latitude of the array for each time step (radians); this can be NULL is parallactic = 0

  • *d_primay_beam_J00[inout] The gains for the north-south beam

  • *d_primay_beam_J01[inout] The leakages for the north-south beam

  • *d_primay_beam_J10[inout] The leakages for the east-west beam

  • *d_primay_beam_J11[inout] The gains for the east-west beam