diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d5425a4d15292c78995ff5fb66aabc5c1ef0c0f3..e451d797ee0ba1fd4ded534b5448d0f935580d1e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -94,7 +94,7 @@ build-compare-oskar:
# Run OSKAR comparison, set some env variables for this session
- export NPIXELS=8 && export APPLY_TRANSPOSE=OFF && MAX_ORDER=3 && TOLERANCE=1e-12
- cd /opt/everybeam/build
- - make VERBOSE=1 comparison-oskar
+ - make VERBOSE=1 comparison-oskar-basefunctions
# Run OSKAR stationresponse comparison
- export NPIXELS=32 && TOLERANCE=1e-5
- make VERBOSE=1 comparison-oskar-stationresponse
diff --git a/cpp/msv3readutils.cc b/cpp/msv3readutils.cc
index 3e80ae69e61ed9342fa2cc5ad138af124bb2b797..ddcfe8519bdeae8677b11a02a93a5e39cd536265 100644
--- a/cpp/msv3readutils.cc
+++ b/cpp/msv3readutils.cc
@@ -73,44 +73,12 @@ using namespace casacore;
vector3r_t TransformToFieldCoordinates(
const vector3r_t &position, const Antenna::CoordinateSystem::Axes &axes);
-// inline Antenna::CoordinateSystem ReadCoordinateSystem(
-// const casacore::Table &table, unsigned int id) {
-// casacore::ArrayQuantColumn<casacore::Double> c_position(table, "POSITION",
-// "m");
-// casacore::ArrayQuantColumn<casacore::Double> c_axes(table,
-// "COORDINATE_SYSTEM",
-// "m");
-//
-// // Read antenna field center (ITRF).
-// casacore::Vector<casacore::Quantity> aips_position = c_position(id);
-// assert(aips_position.size() == 3);
-//
-// vector3r_t position = {{aips_position(0).getValue(),
-// aips_position(1).getValue(),
-// aips_position(2).getValue()}};
-//
-// // Read antenna field coordinate axes (ITRF).
-// casacore::Matrix<casacore::Quantity> aips_axes = c_axes(id);
-// assert(aips_axes.shape().isEqual(casacore::IPosition(2, 3, 3)));
-//
-// vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(),
-// aips_axes(2, 0).getValue()}};
-// vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(),
-// aips_axes(2, 1).getValue()}};
-// vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(),
-// aips_axes(2, 2).getValue()}};
-//
-// Antenna::CoordinateSystem coordinate_system = {position, {p, q, r}};
-// return coordinate_system;
-// }
-
BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int id,
ElementResponse::Ptr element_response) {
Antenna::CoordinateSystem coordinate_system =
common::ReadCoordinateSystem(table, id);
BeamFormer::Ptr beam_former(
new BeamFormerIdenticalAntennas(coordinate_system));
- // BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
@@ -142,20 +110,6 @@ BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int id,
return beam_former;
}
-vector3r_t ReadStationPhaseReference(const Table &table, unsigned int id);
-// {
-// vector3r_t phase_reference = {0.0, 0.0, 0.0};
-// const string columnName("LOFAR_PHASE_REFERENCE");
-// if (common::HasColumn(table, columnName)) {
-// ROScalarMeasColumn<MPosition> c_reference(table, columnName);
-// MPosition mReference =
-// MPosition::Convert(c_reference(id), MPosition::ITRF)();
-// MVPosition mvReference = mReference.getValue();
-// phase_reference = {mvReference(0), mvReference(1), mvReference(2)};
-// }
-// return phase_reference;
-// }
-
Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
const ElementResponseModel model) {
ROMSAntennaColumns antenna(ms.antenna());
@@ -173,9 +127,6 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
// Create station.
Station::Ptr station(new Station(name, position, model));
- // Read phase reference position (if available).
- // station->SetPhaseReference(ReadStationPhaseReference(ms.antenna(), id));
-
Table tab_phased_array = common::GetSubTable(ms, "PHASED_ARRAY");
// The Station will consist of a BeamFormer that combines the fields
@@ -184,10 +135,6 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
auto beam_former =
ReadMSv3AntennaField(tab_phased_array, id, station->GetElementResponse());
- // TODO
- // If There is only one field, the top level beamformer is not needed
- // and the station antenna can be set the the beamformer of the field
-
station->SetAntenna(beam_former);
size_t field_id = 0;
@@ -195,7 +142,6 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
Antenna::CoordinateSystem coordinate_system =
common::ReadCoordinateSystem(tab_phased_array, field_id);
auto element_response = station->GetElementResponse();
- // TODO: rotate coordinate system for antenna
auto element = Element::Ptr(
new Element(coordinate_system, element_response, element_id));
station->SetElement(element);
diff --git a/cpp/oskar/oskarelementresponse.cc b/cpp/oskar/oskarelementresponse.cc
index 2081171df960823d6ce5caa1c700c45817641061..f1669f35db033da1254d7a9369d5165f2c8b516b 100644
--- a/cpp/oskar/oskarelementresponse.cc
+++ b/cpp/oskar/oskarelementresponse.cc
@@ -51,7 +51,15 @@ void OSKARElementResponseSphericalWave::Response(
std::complex<double>* alpha_ptr = dataset->GetAlphaPtr(element_id);
double phi_x = phi;
- double phi_y = phi; // + M_PI_2;
+ double phi_y = phi;
+
+ // TODO: phi_x and phi_y can have different values if there is only one set
+ // of coefficients that is is used for both dipoles.
+ // In that case it is assumed the Y dipole rotated 90deg with respect
+ // to the X dipole, so then phi_y = phi+ M_PI_2.
+ // That case needs to be detected when the coefficients are read,
+ // and here phi_y needs to be set accordingly.
+
oskar_evaluate_spherical_wave_sum_double(1, &theta, &phi_x, &phi_y, l_max,
alpha_ptr, response_ptr);
}
diff --git a/cpp/telescope/oskar.cc b/cpp/telescope/oskar.cc
index 0e62d3f9acfcb4c3af0005db62ff95e997c4fdd4..f3b9c151d11c7a9e7a36f26f413a23117f10924b 100644
--- a/cpp/telescope/oskar.cc
+++ b/cpp/telescope/oskar.cc
@@ -18,13 +18,19 @@ OSKAR::OSKAR(MeasurementSet &ms, const Options &options)
}
std::unique_ptr<griddedresponse::GriddedResponse> OSKAR::GetGriddedResponse(
- const coords::CoordinateSystem &coordinate_system){
- // Get and return GriddedResponse ptr
- // std::unique_ptr<griddedresponse::GriddedResponse> grid(
- // new griddedresponse::LOFARGrid(this, coordinate_system));
- // // griddedresponse::GriddedResponse grid(LOFARGrid(this,
- // coordinate_system));
- // return grid;
+ const coords::CoordinateSystem &coordinate_system) {
+ throw std::runtime_error(
+ "GetGriddedResponse() is not implemented for OSKAR Telescope");
+
+ // TODO: return an OSKARGrid here, in a similar way to the commented out code
+ // below
+
+ // Get and return GriddedResponse ptr
+ // std::unique_ptr<griddedresponse::GriddedResponse> grid(
+ // new griddedresponse::LOFARGrid(this, coordinate_system));
+ // // griddedresponse::GriddedResponse grid(LOFARGrid(this,
+ // coordinate_system));
+ // return grid;
};
Station::Ptr OSKAR::ReadStation(const MeasurementSet &ms, std::size_t id,
diff --git a/demo/comparison-oskar/CMakeLists.txt b/demo/comparison-oskar/CMakeLists.txt
index 5cb5f03e127bde12fdf3bdfdf78e814d7d62b68a..c749680334aca52d4cbaf71f9ab3f84efb2e7ac7 100644
--- a/demo/comparison-oskar/CMakeLists.txt
+++ b/demo/comparison-oskar/CMakeLists.txt
@@ -1,13 +1,13 @@
#------------------------------------------------------------------------------
# CMake file for compiling a comparison between OSKAR and EveryBeam
-add_executable(comparison-oskar-generate-beampattern main.cpp)
-target_link_libraries(comparison-oskar-generate-beampattern oskar)
+add_executable(make_element_response_image make_element_response_image.cpp)
+target_link_libraries(make_element_response_image oskar)
-add_executable(comparison-oskar-generate-station-response stationresponse.cpp)
-target_link_libraries(comparison-oskar-generate-station-response everybeam)
+add_executable(make_station_response_image make_station_response_image.cpp)
+target_link_libraries(make_station_response_image everybeam)
# Required to get the config.h header
-target_include_directories(comparison-oskar-generate-station-response PRIVATE "${CMAKE_BINARY_DIR}")
+target_include_directories(make_station_response_image PRIVATE "${CMAKE_BINARY_DIR}")
file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/layout.txt"
@@ -25,17 +25,17 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy_directory
#------------------------------------------------------------------------------
# comparison-oskar knits together the cpp code and the python scripts
-add_custom_target(comparison-oskar
+add_custom_target(comparison-oskar-basefunctions
COMMAND ${CMAKE_COMMAND} -E env
EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts"
"${CMAKE_CURRENT_SOURCE_DIR}/generate_basefunction_plots.sh"
- DEPENDS comparison-oskar-generate-beampattern
+ DEPENDS make_element_response_image
)
add_custom_target(comparison-oskar-station-response
COMMAND ${CMAKE_COMMAND} -E env
EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts:${CMAKE_SOURCE_DIR}/scripts/misc"
"${CMAKE_CURRENT_SOURCE_DIR}/compare_stationresponse.sh"
- DEPENDS comparison-oskar-generate-station-response
+ DEPENDS make_station_response_image
)
diff --git a/demo/comparison-oskar/compare_stationresponse.py b/demo/comparison-oskar/compare_stationresponse.py
index 3628300c3d084ebecf2298f5f1ca284ebb982b91..f73319133efef7329f0de868262856d2d8cea134 100644
--- a/demo/comparison-oskar/compare_stationresponse.py
+++ b/demo/comparison-oskar/compare_stationresponse.py
@@ -13,7 +13,7 @@ npixels = int(os.environ["NPIXELS"]) if "NPIXELS" in os.environ else 256
run_oskar_simulation.main(npixels)
subprocess.check_call(["add_beaminfo.py", "skalowmini-coef.MS", "skalowmini-coef.tm"])
subprocess.check_call(["oskar_csv_to_hdf5.py", "skalowmini-coef.tm", "oskar.h5"])
-subprocess.check_call(["./comparison-oskar-generate-station-response", str(npixels)])
+subprocess.check_call(["./make_station_response_image", str(npixels)])
A = read_oskar_beams()
B = np.load('station-response.npy')
diff --git a/demo/comparison-oskar/generate_basefunction_plots.py b/demo/comparison-oskar/generate_basefunction_plots.py
index b9d048b11a96e2ae8948f408a6f3c6fe79eeeca8..3c11c12f436a9d7f94dc4c8910b23632c3a159d0 100644
--- a/demo/comparison-oskar/generate_basefunction_plots.py
+++ b/demo/comparison-oskar/generate_basefunction_plots.py
@@ -68,7 +68,7 @@ for em_idx in range(2):
generate_oskar_csv(l * l - 1 + l - m, em_idx)
subprocess.check_call(["oskar_csv_to_hdf5.py", "telescope.tm", "oskar.h5"])
- subprocess.check_call(["comparison-oskar-generate-beampattern", str(npixels)])
+ subprocess.check_call(["make_element_response_image", str(npixels)])
B = np.load("response.npy")
diff --git a/demo/comparison-oskar/main.cpp b/demo/comparison-oskar/make_element_response_image.cpp
similarity index 100%
rename from demo/comparison-oskar/main.cpp
rename to demo/comparison-oskar/make_element_response_image.cpp
diff --git a/demo/comparison-oskar/make_station_response_image.cpp b/demo/comparison-oskar/make_station_response_image.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ebbaacdc2d3d19178584fcff2b8052225946b6b5
--- /dev/null
+++ b/demo/comparison-oskar/make_station_response_image.cpp
@@ -0,0 +1,105 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+
+#include <oskarelementresponse.h>
+
+#include "../../external/npy.hpp" // to save arrays in numpy format
+
+#include "load.h"
+#include "options.h"
+#include "config.h"
+#include "coords/coordutils.h"
+#include "coords/itrfconverter.h"
+
+using namespace everybeam;
+
+int main(int argc, char** argv){
+
+ ElementResponseModel response_model = ElementResponseModel::kOSKARSphericalWave;
+ Options options;
+ options.element_response_model = response_model;
+
+ casacore::MeasurementSet ms("skalowmini-coef.MS");
+
+ casacore::ScalarMeasColumn<casacore::MDirection> referenceDirColumn(
+ ms.field(),
+ casacore::MSField::columnName(casacore::MSFieldEnums::REFERENCE_DIR));
+
+ auto reference_dir = referenceDirColumn(0);
+
+ std::cout << ms.nrow() << std::endl;;
+ auto unique_times_table = ms.sort("TIME", casacore::Sort::Ascending, casacore::Sort::NoDuplicates);
+ std::cout << unique_times_table.nrow() << std::endl;
+
+ casacore::ScalarColumn<double> time_column(ms, "TIME");
+
+ real_t time = time_column(0);
+ std::cout << "time: " << time << std::endl;
+
+ std::cout << reference_dir << std::endl;
+
+ vector3r_t station0 ;
+ vector3r_t tile0;
+ coords::ITRFConverter itrf_converter(time);
+ coords::SetITRFVector(itrf_converter.ToDirection(reference_dir), station0);
+ coords::SetITRFVector(itrf_converter.ToDirection(reference_dir), tile0);
+
+ // Load OSKAR Telescope
+ auto telescope = Load(ms, options);
+ auto &oskar_telescope = *dynamic_cast<telescope::OSKAR*>(telescope.get());
+
+ Station::Ptr station = oskar_telescope.GetStation(0);
+ auto element_response = std::make_shared<everybeam::OSKARElementResponseSphericalWave>("oskar.h5");
+ station->SetResponse(element_response);
+
+ Antenna::Ptr antenna = station->GetAntenna();
+
+ auto p = antenna->coordinate_system_.axes.p;
+ auto q = antenna->coordinate_system_.axes.q;
+ auto r = antenna->coordinate_system_.axes.r;
+
+
+ double freq = 50e6;
+
+ int N;
+ if (argc == 1){
+ N = 256;
+ }
+ else{
+ N = atoi(argv[1]);
+ }
+
+ std::vector<std::complex<double>> result(N*N*2*2);
+ typedef std::complex<double>result_arr_t[N][N][2][2];
+ result_arr_t &result_arr = * (result_arr_t*) result.data();
+
+
+ for(int i=0; i<N; ++i) {
+ std::cout << i << std::endl;
+ double x = (2.0*i)/(N-1) - 1.0;
+ for(int j=0; j<N; ++j) {
+ double y = (2.0*j)/(N-1) - 1.0;
+
+ double z = sqrt(1.0 - x*x - y*y);
+
+ const vector3r_t direction = {
+ x*p[0] + y*q[0] + z*r[0],
+ x*p[1] + y*q[1] + z*r[1],
+ x*p[2] + y*q[2] + z*r[2]
+ };
+ real_t freq0 = 50e6;
+
+ auto result = station->Response(time, freq, direction, freq0, station0, tile0, false);
+ result_arr[i][j][0][0] = result[0][0];
+ result_arr[i][j][0][1] = result[0][1];
+ result_arr[i][j][1][0] = result[1][0];
+ result_arr[i][j][1][1] = result[1][1];
+
+ }
+ }
+
+ const long unsigned leshape [] = {(long unsigned int) N, (long unsigned int) N, 2, 2};
+ npy::SaveArrayAsNumpy("station-response.npy", false, 4, leshape, result);
+}
+
diff --git a/demo/comparison-oskar/read_oskar_beams.py b/demo/comparison-oskar/read_oskar_beams.py
index 53dc898736fd3e8ab5f18722dffd736b37c9b12c..424c3d4ea2dd54e33331918b157b9f0000493622 100644
--- a/demo/comparison-oskar/read_oskar_beams.py
+++ b/demo/comparison-oskar/read_oskar_beams.py
@@ -27,7 +27,12 @@ def read_oskar_beams():
if A is None:
N = d.shape[-1]
A = np.zeros((N, N, 2, 2), dtype=np.complex128)
+
+ # One axis in the OSKAR fits files runs in opposite direction compared to the
+ # numpy arrays made by make_element_response_image.cpp
+ # we flip that axis here to make them the same
A[:, :, i, j] = d[::-1, :]
+
N = A.shape[0]
for i in range(N):
@@ -37,10 +42,6 @@ def read_oskar_beams():
theta = np.arcsin(np.sqrt(x * x + y * y))
phi = np.arctan2(y, x)
- # Need to swap the sign of phi to get the transformation right
- # Still need to figure out why
- # phi = -phi
-
e_theta = np.array([[np.cos(phi)], [np.sin(phi)]])
e_phi = np.array([[-np.sin(phi)], [np.cos(phi)]])
diff --git a/demo/comparison-oskar/skalowmini-coef.tm/README.md b/demo/comparison-oskar/skalowmini-coef.tm/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..5c4d065d41369e692a60a977f68725ca107adfb8
--- /dev/null
+++ b/demo/comparison-oskar/skalowmini-coef.tm/README.md
@@ -0,0 +1,4 @@
+The element_pattern_spherical_wave*.txt files in this directory are copied
+from the [OSKAR example data](https://github.com/OxfordSKA/OSKAR/files/1984210/OSKAR-2.7-Example-Data.zip)
+
+
diff --git a/demo/comparison-oskar/stationresponse.cpp b/demo/comparison-oskar/stationresponse.cpp
deleted file mode 100644
index 0f34e492b51bedde6f9f7cf5c4902f41fd4d1d7e..0000000000000000000000000000000000000000
--- a/demo/comparison-oskar/stationresponse.cpp
+++ /dev/null
@@ -1,236 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include <cstdlib>
-
-#include <oskarelementresponse.h>
-
-#include "../../external/npy.hpp" // to save arrays in numpy format
-
-#include "load.h"
-#include "options.h"
-#include "config.h"
-#include "coords/coordutils.h"
-#include "coords/itrfconverter.h"
-
-using namespace everybeam;
-
-int main(int argc, char** argv){
-
- ElementResponseModel response_model = ElementResponseModel::kOSKARSphericalWave;
- Options options;
- options.element_response_model = response_model;
-
-// casacore::MeasurementSet ms("/home/vdtol/skalowmini/skalowmini-coef1.MS");
- casacore::MeasurementSet ms("skalowmini-coef.MS");
-
- casacore::ScalarMeasColumn<casacore::MDirection> referenceDirColumn(
- ms.field(),
- casacore::MSField::columnName(casacore::MSFieldEnums::REFERENCE_DIR));
-
- auto reference_dir = referenceDirColumn(0);
-
- std::cout << ms.nrow() << std::endl;;
- auto unique_times_table = ms.sort("TIME", casacore::Sort::Ascending, casacore::Sort::NoDuplicates);
- std::cout << unique_times_table.nrow() << std::endl;
-
- casacore::ScalarColumn<double> time_column(ms, "TIME");
-
- real_t time = time_column(0);
- std::cout << "time: " << time << std::endl;
-
- std::cout << reference_dir << std::endl;
-
- vector3r_t station0 ;
- vector3r_t tile0;
- coords::ITRFConverter itrf_converter(time);
- coords::SetITRFVector(itrf_converter.ToDirection(reference_dir), station0);
- coords::SetITRFVector(itrf_converter.ToDirection(reference_dir), tile0);
-
- // Load OSKAR Telescope
- auto telescope = Load(ms, options);
- auto &oskar_telescope = *dynamic_cast<telescope::OSKAR*>(telescope.get());
-
- Station::Ptr station = oskar_telescope.GetStation(0);
-
- Antenna::Ptr antenna = station->GetAntenna();
-
- auto p = antenna->coordinate_system_.axes.p;
- auto q = antenna->coordinate_system_.axes.q;
- auto r = antenna->coordinate_system_.axes.r;
-
-
- double freq = 50e6;
-
- int N;
- if (argc == 1){
- N = 256;
- }
- else{
- N = atoi(argv[1]);
- }
-
- std::vector<std::complex<double>> result(N*N*2*2);
- typedef std::complex<double>result_arr_t[N][N][2][2];
- result_arr_t &result_arr = * (result_arr_t*) result.data();
-
-
- for(int i=0; i<N; ++i) {
- std::cout << i << std::endl;
- double x = (2.0*i)/(N-1) - 1.0;
- for(int j=0; j<N; ++j) {
- double y = (2.0*j)/(N-1) - 1.0;
-
- double z = sqrt(1.0 - x*x - y*y);
-
- const vector3r_t direction = {
- x*p[0] + y*q[0] + z*r[0],
- x*p[1] + y*q[1] + z*r[1],
- x*p[2] + y*q[2] + z*r[2]
- };
- real_t freq0 = 50e6;
-
- auto result = station->Response(time, freq, direction, freq0, station0, tile0, false);
- result_arr[i][j][0][0] = result[0][0];
- result_arr[i][j][0][1] = result[0][1];
- result_arr[i][j][1][0] = result[1][0];
- result_arr[i][j][1][1] = result[1][1];
-
-// std::cout << result_arr[i][j][0][0] << ", " << result_arr[i][j][1][0] << std::endl;
-
-// auto result = station->ArrayFactor(time, freq, direction, freq0, station0, tile0);
-// result_arr[i][j][0][0] = result[0];
-// result_arr[i][j][0][1] = 0.0;
-// result_arr[i][j][1][0] = 0.0;
-// result_arr[i][j][1][1] = result[1];
-
- }
- }
-
- const long unsigned leshape [] = {(long unsigned int) N, (long unsigned int) N, 2, 2};
- npy::SaveArrayAsNumpy("station-response.npy", false, 4, leshape, result);
-}
-
-
-/*
-
-#include "../options.h"
-#include "../griddedresponse/lofargrid.h"
-#include "../elementresponse.h"
-#include "../../external/npy.hpp"
-
-#include <complex>
-#include <cmath>
-
-using namespace everybeam;
-
-BOOST_AUTO_TEST_CASE(load_lofar) {
- ElementResponseModel response_model = ElementResponseModel::kHamaker;
- Options options;
-
- casacore::MeasurementSet ms(LOFAR_MOCK_MS);
-
- // Load LOFAR Telescope
- std::unique_ptr<telescope::Telescope> telescope =
- Load(ms, response_model, options);
-
- // Assert if we indeed have a LOFAR pointer
- BOOST_CHECK(nullptr != dynamic_cast<telescope::LOFAR*>(telescope.get()));
-
- // Assert if correct number of stations
- std::size_t nstations = 70;
- BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
-
- // Assert if GetStation(stationd_id) behaves properly
- BOOST_CHECK_EQUAL(telescope->GetStation(0)->GetName(), "CS001HBA0");
-
- // Properties extracted from MS
- double time = 4929192878.008341;
- double frequency = 138476562.5;
- std::size_t width(4), height(4);
- double ra(2.15374123), dec(0.8415521), dl(0.5 * M_PI / 180.),
- dm(0.5 * M_PI / 180.), shift_l(0.), shift_m(0.);
-
- coords::CoordinateSystem coord_system = {.width = width,
- .height = height,
- .ra = ra,
- .dec = dec,
- .dl = dl,
- .dm = dm,
- .phase_centre_dl = shift_l,
- .phase_centre_dm = shift_m};
- std::unique_ptr<griddedresponse::GriddedResponse> grid_response =
- telescope->GetGriddedResponse(coord_system);
- BOOST_CHECK(nullptr !=
- dynamic_cast<griddedresponse::LOFARGrid*>(grid_response.get()));
-
- // Define buffer and get gridded responses
- std::vector<std::complex<float>> antenna_buffer_single(
- grid_response->GetBufferSize(1));
- grid_response->CalculateStation(antenna_buffer_single.data(), time, frequency,
- 23);
- BOOST_CHECK_EQUAL(antenna_buffer_single.size(),
- std::size_t(width * height * 2 * 2));
-
- // LOFARBeam output at pixel (2,2):
- std::vector<std::complex<float>> lofar_p22 = {{-0.175908, -0.000478397},
- {-0.845988, -0.00121503},
- {-0.89047, -0.00125383},
- {0.108123, -5.36076e-05}};
-
- // Compare with everybeam
- std::size_t offset_22 = (2 + 2 * height) * 4;
- for (std::size_t i = 0; i < 4; ++i) {
- BOOST_CHECK(std::abs(antenna_buffer_single[offset_22 + i] - lofar_p22[i]) <
- 1e-4);
- }
-
- // LOFARBeam output at pixel (1,3):
- std::vector<std::complex<float>> lofar_p13 = {{-0.158755, -0.000749433},
- {-0.816165, -0.00272568},
- {-0.863389, -0.00283979},
- {0.0936919, 0.000110673}};
-
- // Compare with everybeam
- std::size_t offset_13 = (1 + 3 * height) * 4;
- for (std::size_t i = 0; i < 4; ++i) {
- BOOST_CHECK(std::abs(antenna_buffer_single[offset_13 + i] - lofar_p13[i]) <
- 1e-4);
- }
-
- // std::vector<std::complex<float>> antenna_buffer_all(
- // grid_response->GetBufferSize(telescope->GetNrStations()));
- // grid_response->CalculateAllStations(antenna_buffer_all.data(), time,
- // frequency);
- // BOOST_CHECK_EQUAL(
- // antenna_buffer_all.size(),
- // std::size_t(telescope->GetNrStations() * width * height * 2 * 2));
-
- // Test with differential beam, single
- Options options_diff_beam;
- options_diff_beam.use_differential_beam = true;
-
- // Load LOFAR Telescope
- std::unique_ptr<telescope::Telescope> telescope_diff_beam =
- Load(ms, response_model, options_diff_beam);
-
- std::unique_ptr<griddedresponse::GriddedResponse> grid_response_diff_beam =
- telescope_diff_beam->GetGriddedResponse(coord_system);
-
- std::vector<std::complex<float>> antenna_buffer_diff_beam(
- grid_response_diff_beam->GetBufferSize(1));
- grid_response_diff_beam->CalculateStation(antenna_buffer_diff_beam.data(),
- time, frequency, 15);
-
- double norm_jones_mat = 0.;
- for (std::size_t i = 0; i < 4; ++i) {
- norm_jones_mat += std::norm(antenna_buffer_diff_beam[offset_22 + i]);
- }
- BOOST_CHECK(std::abs(norm_jones_mat - 2.) < 1e-6);
-
- // Print to np array
- // const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
- // npy::SaveArrayAsNumpy("station_responses.npy", false, 4, leshape,
- // antenna_buffer_single);
-}
-
-*/