Skip to content
Snippets Groups Projects
Select Git revision
  • 2f8a54373b7068f9955284c742e9f8bb85f48cbd
  • master default protected
  • gec-113-static-ska-sdp-func
  • rtd-ubuntu24-plucky
  • ast-1644-temp-fix
  • ast-1600-fix-beam-for-meerkat-ska-mid
  • readthedocs-c++17
  • mwa_python_wrapper
  • ast-1384-remove-element-index-argument
  • ast-1509-fix-polarization-orientation-in-gridded-response
  • add-test-for-squared-mueller-matrices
  • 1493-extend-python-bindings
  • ast-1493-implement-response-dishpoint-1
  • ast-1325-prototype-ska-beam-model-interface
  • ast-1416-oskar-ska-sdp-func-1
  • ast-1386-create-default-element
  • ast-1384-fix-sorted-frequencies-check-sphericalharmonicsresponse-1
  • ast-1111-add-vector-bindings
  • ast-973-add-test-for-lobes-coefficients
  • ast-645-add-beam-normalisation-mode-preapplied
  • disable-element-beam-1
  • v0.7.1
  • v0.7.0
  • v0.6.2
  • v0.6.1
  • v0.6.0
  • v0.5.8
  • v0.5.7
  • v0.5.6
  • v0.5.5
  • v0.5.4
  • v0.5.3
  • v0.5.2
  • v0.5.1
  • v0.4.0
  • v0.3.1
  • v0.3.0
  • v0.2.0
  • v0.1.3
  • v0.1.2
  • v0.1.1
41 results

tgmrt.cc

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    tgmrt.cc 7.32 KiB
    // Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
    // SPDX-License-Identifier: GPL-3.0-or-later
    
    #include <boost/test/unit_test.hpp>
    #include <boost/test/tools/floating_point_comparison.hpp>
    
    #include "config.h"
    
    #include "../load.h"
    #include "../griddedresponse/dishgrid.h"
    #include "../telescope/telescope.h"
    #include "../telescope/dish.h"
    
    namespace {
    
    const double kTime = 5.08859e+09;
    // There's a "reference" fits file at 383.325 MHz that I've used to
    // determine some properties such as the (approx) FWHM that are later
    // used in these tests
    const double kFrequency = 383.325e6;
    const size_t kWidth = 64;
    const size_t kHeight = 64;
    
    aocommon::CoordinateSystem GetCoordinateSystem() {
      aocommon::CoordinateSystem coord_system;
      coord_system.width = kWidth;
      coord_system.height = kHeight;
      coord_system.ra = -156.837226063 * (M_PI / 180.0);
      coord_system.dec = 50.4691390331 * (M_PI / 180.0);
      coord_system.dl = 1.0 * (M_PI / 180.0 / 60.0);
      coord_system.dm = 1.0 * (M_PI / 180.0 / 60.0);
      coord_system.l_shift = 0.0;
      coord_system.m_shift = 0.0;
      return coord_system;
    }
    
    struct GmrtFixture {
      GmrtFixture() : ms(GMRT_MOCK_MS), telescope(everybeam::Load(ms, options)) {
        grid_response = telescope->GetGriddedResponse(GetCoordinateSystem());
      }
      ~GmrtFixture(){};
      casacore::MeasurementSet ms;
      everybeam::Options options;
      std::unique_ptr<everybeam::telescope::Telescope> telescope;
      std::unique_ptr<everybeam::griddedresponse::GriddedResponse> grid_response;
    
      size_t Index(size_t x, size_t y) const { return 4 * (kWidth * y + x); }
    
      void CheckStation(const std::complex<float>* single_antenna_buffer) const {
        const size_t size = kWidth * kHeight * 2 * 2;
        for (size_t i = 0; i != size; i += 4) {
          // Check if off diagonals are zero:
          BOOST_CHECK_LT(std::fabs(single_antenna_buffer[i + 1].real()), 1e-4);
          BOOST_CHECK_LT(std::fabs(single_antenna_buffer[i + 2].real()), 1e-4);
          // Because X and Y are calculated in the same way, they should always be
          // identical
          BOOST_CHECK_CLOSE_FRACTION(single_antenna_buffer[i].real(),
                                     single_antenna_buffer[i + 3].real(), 1e-4);
        }
        // Check if all phases (/imaginary values) are zero:
        for (size_t i = 0; i != size; ++i) {
          BOOST_CHECK_LT(std::fabs(single_antenna_buffer[i].imag()), 1e-4);
        }
    
        // Check if the centre pixel is approximately the unit matrix
        // (Due to interpolation the diagonal gains are not exactly 1; it was
        // 0.999789 when I checked).
        const size_t centre_index = Index(kWidth / 2, kHeight / 2);
        BOOST_CHECK_CLOSE_FRACTION(single_antenna_buffer[centre_index].real(), 1.0,
                                   1e-2);
    
        // The FWHM is approx at the edge. These are (un-squared) voltages, hence
        // its value is about sqrt(0.5) (=~0.7).
        const float value_at_edge = 0.697981417;
        const size_t edge_index = Index(kWidth / 2, 0);
        BOOST_CHECK_CLOSE_FRACTION(single_antenna_buffer[edge_index].real(),
                                   value_at_edge, 1e-2);
        BOOST_CHECK_CLOSE_FRACTION(single_antenna_buffer[edge_index + 3].real(),
                                   value_at_edge, 1e-2);
      }
    };
    
    void CheckUndersampledIntegratedResponse(size_t undersampling_factor) {
      GmrtFixture gmrt;
      const std::vector<double> time_array{kTime};
      const size_t n_antennas = gmrt.telescope->GetNrStations();
      const size_t n_baselines = n_antennas * (n_antennas + 1) / 2;
      const std::vector<double> baseline_weights(n_baselines * time_array.size(),
                                                 1.0);
    
      const std::vector<aocommon::HMC4x4> result =
          gmrt.grid_response->UndersampledIntegratedResponse(
              everybeam::CorrectionMode::kFull, time_array, kFrequency, 0,
              undersampling_factor, baseline_weights);
      const size_t undersampled_width = kWidth / undersampling_factor;
      const size_t undersampled_height = kHeight / undersampling_factor;
      BOOST_REQUIRE_EQUAL(result.size(), undersampled_width * undersampled_height);
    
      for (size_t element_index = 0; element_index != 4; ++element_index) {
        const float lowres_centre_pixel =
            result[undersampled_width * (undersampled_height / 2) +
                   undersampled_width / 2]
                .Data(element_index);
        const size_t off_centre_y =
            undersampled_height / 2 - 8 / undersampling_factor;
        const float lowres_off_centre_pixel =
            result[undersampled_width * off_centre_y + undersampled_width / 2].Data(
                element_index);
        const float lowres_edge_pixel =
            result[undersampled_width / 2].Data(element_index);
    
        std::vector<float> buffer(kWidth * kHeight);
        gmrt.grid_response->UpsampleResponse(buffer.data(), element_index, kWidth,
                                             kHeight, result, undersampling_factor);
    
        const float upsampled_centre_pixel =
            buffer[kWidth * (kHeight / 2) + kWidth / 2];
        const float upsampled_edge_pixel = buffer[kWidth / 2];
    
        if (element_index == 0 || element_index == 3) {
          BOOST_CHECK_CLOSE_FRACTION(lowres_centre_pixel, 1.0, 1e-2);
          BOOST_CHECK_CLOSE_FRACTION(lowres_off_centre_pixel, 0.959, 1e-2);
          BOOST_CHECK_CLOSE_FRACTION(lowres_edge_pixel, 0.487, 1e-2);
          BOOST_CHECK_CLOSE_FRACTION(upsampled_centre_pixel, 1.0, 1e-2);
          BOOST_CHECK_CLOSE_FRACTION(upsampled_edge_pixel, 0.487, 1e-2);
        } else {
          BOOST_CHECK_LT(std::fabs(lowres_centre_pixel), 1e-2);
          BOOST_CHECK_LT(std::fabs(lowres_off_centre_pixel), 1e-2);
          BOOST_CHECK_LT(std::fabs(lowres_edge_pixel), 1e-2);
          BOOST_CHECK_LT(std::fabs(upsampled_centre_pixel), 1e-2);
          BOOST_CHECK_LT(std::fabs(upsampled_edge_pixel), 1e-2);
        }
      }
    }
    
    }  // namespace
    
    BOOST_AUTO_TEST_SUITE(tgmrt)
    
    BOOST_FIXTURE_TEST_CASE(single_antenna_gridded_response, GmrtFixture) {
      BOOST_CHECK(nullptr != dynamic_cast<everybeam::griddedresponse::DishGrid*>(
                                 grid_response.get()));
    
      // Define buffer and get gridded responses
      std::vector<std::complex<float>> single_antenna_buffer(
          grid_response->GetStationBufferSize(1));
      grid_response->Response(everybeam::CorrectionMode::kFull,
                              single_antenna_buffer.data(), kTime, kFrequency, 0,
                              0);
      BOOST_REQUIRE_EQUAL(single_antenna_buffer.size(),
                          std::size_t(kWidth * kHeight * 2 * 2));
    
      CheckStation(single_antenna_buffer.data());
    }
    
    BOOST_FIXTURE_TEST_CASE(all_antennas_gridded_response, GmrtFixture) {
      std::vector<std::complex<float>> buffer(
          grid_response->GetStationBufferSize(telescope->GetNrStations()));
      grid_response->ResponseAllStations(everybeam::CorrectionMode::kFull,
                                         buffer.data(), kTime, kFrequency, 0);
      const size_t single_antenna_size = kWidth * kHeight * 4;
      const size_t expected_size = telescope->GetNrStations() * single_antenna_size;
      BOOST_REQUIRE_EQUAL(buffer.size(), expected_size);
    
      for (size_t antenna_index = 0; antenna_index != telescope->GetNrStations();
           ++antenna_index) {
        const std::complex<float>* single_antenna_buffer =
            &buffer[single_antenna_size * antenna_index];
        CheckStation(single_antenna_buffer);
      }
    }
    
    BOOST_AUTO_TEST_CASE(undersampled_integrated_response_factor_1) {
      CheckUndersampledIntegratedResponse(1);
    }
    
    BOOST_AUTO_TEST_CASE(undersampled_integrated_response_factor_2) {
      CheckUndersampledIntegratedResponse(2);
    }
    
    BOOST_AUTO_TEST_SUITE_END()