Skip to content
Snippets Groups Projects
Commit d227ec85 authored by Bas van der Tol's avatar Bas van der Tol
Browse files

Clarify code: remove commented out code, update comments, rename files

parent 2a87c980
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 147 additions and 77 deletions
.gitlab-ci.yml
+ 1
− 1
View file @ d227ec85
......@@ -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
......
cpp/msv3readutils.cc
+ 0
− 54
View file @ d227ec85
......@@ -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);
......
cpp/oskar/oskarelementresponse.cc
+ 9
− 1
View file @ d227ec85
......@@ -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);
}
......
cpp/telescope/oskar.cc
+ 13
− 7
View file @ d227ec85
......@@ -19,6 +19,12 @@ OSKAR::OSKAR(MeasurementSet &ms, const Options &options)
std::unique_ptr<griddedresponse::GriddedResponse> OSKAR::GetGriddedResponse(
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));
......
demo/comparison-oskar/CMakeLists.txt
+ 8
− 8
View file @ d227ec85
#------------------------------------------------------------------------------
# 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
)
demo/comparison-oskar/compare_stationresponse.py
+ 1
− 1
View file @ d227ec85
......@@ -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')
......
demo/comparison-oskar/generate_basefunction_plots.py
+ 1
− 1
View file @ d227ec85
......@@ -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")
......
demo/comparison-oskar/stationresponse.cpp → demo/comparison-oskar/make_station_response_image.cpp
+ 105
− 0
View file @ d227ec85
......@@ -20,7 +20,6 @@ int main(int argc, char** argv){
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(
......@@ -51,6 +50,8 @@ int main(int argc, char** argv){
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();
......@@ -95,14 +96,6 @@ int main(int argc, char** argv){
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];
}
}
......@@ -110,127 +103,3 @@ int main(int argc, char** argv){
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);
}
*/
demo/comparison-oskar/read_oskar_beams.py
+ 5
− 4
View file @ d227ec85
......@@ -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)]])
......
demo/comparison-oskar/skalowmini-coef.tm/README.md 0 → 100644
+ 4
− 0
View file @ d227ec85
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment