Skip to content
Snippets Groups Projects
Commit 30854f3b authored by Andre Offringa's avatar Andre Offringa
Browse files

Implement Airy disk function

parent 7e747661
No related branches found
Tags v0.5.3
1 merge request!292Implement Airy disk function
Pipeline #60474 passed with warnings
Showing
with 370 additions and 65 deletions
......@@ -7,7 +7,7 @@ cmake_minimum_required(VERSION 3.8)
#------------------------------------------------------------------------------
# Set version name and project number
set(EVERYBEAM_VERSION 0.5.2)
set(EVERYBEAM_VERSION 0.5.3)
if(EVERYBEAM_VERSION MATCHES "^([0-9]+)\\.([0-9]+)\\.([0-9]+)")
set(EVERYBEAM_VERSION_MAJOR "${CMAKE_MATCH_1}")
set(EVERYBEAM_VERSION_MINOR "${CMAKE_MATCH_2}")
......
......@@ -56,12 +56,35 @@ add_library(
sphericalharmonicsresponse.cc
sphericalharmonicsresponsefixeddirection.cc
station.cc
telescope/alma.cc
telescope/lofar.cc
telescope/dish.cc
telescope/lofar.cc
telescope/mwa.cc
telescope/phasedarray.cc
telescope/oskar.cc
telescope/skamid.cc)
telescope/skamid.cc
griddedresponse/aartfaacgrid.cc
griddedresponse/airygrid.cc
griddedresponse/dishgrid.cc
griddedresponse/griddedresponse.cc
griddedresponse/mwagrid.cc
griddedresponse/skamidgrid.cc
pointresponse/airypoint.cc
pointresponse/dishpoint.cc
pointresponse/mwapoint.cc
pointresponse/phasedarraypoint.cc
pointresponse/skamidpoint.cc
circularsymmetric/atcacoefficients.cc
circularsymmetric/gmrtcoefficients.cc
circularsymmetric/meerkatcoefficients.cc
circularsymmetric/vlacoefficients.cc
circularsymmetric/voltagepattern.cc
# Phased array telescopes (SKA, LOFAR)
griddedresponse/phasedarraygrid.cc
# MWA(beam) related
mwabeam/tilebeam2016.cc
mwabeam/beam2016implementation.cc)
# Make sure that when other targets within this project link against the everybeam target,
# they can find the include files.
......
......@@ -11,17 +11,27 @@ using aocommon::ImageCoordinates;
using aocommon::UVector;
using everybeam::circularsymmetric::VoltagePattern;
namespace {
double LmMaxSquared(double frequency_hz, double maximum_radius_arc_min) {
const double factor =
(180.0 / M_PI) * 60.0 * frequency_hz * 1e-9; // arcminutes * GHz
const double rmax = maximum_radius_arc_min / factor;
return rmax * rmax;
}
} // namespace
void VoltagePattern::EvaluatePolynomial(const UVector<double>& coefficients,
double reference_frequency,
bool invert) {
// This comes from casa's: void PBMath1DIPoly::fillPBArray(), wideband case
const size_t nsamples = 10000;
constexpr size_t nsamples = 10000;
const size_t nfreq = frequencies_.size();
const size_t ncoef = coefficients.size() / nfreq;
values_.resize(nsamples * nfreq);
inverse_increment_radius_ = double(nsamples - 1) / maximum_radius_arc_min_;
double* output = values_.data();
const double referenced_increment =
std::sqrt(reference_frequency_ / 1e9) / inverse_increment_radius_;
std::sqrt(reference_frequency / 1e9) / inverse_increment_radius_;
for (size_t n = 0; n != nfreq; n++) {
const double* freq_coefficients = &coefficients[n * ncoef];
for (size_t i = 0; i < nsamples; i++) {
......@@ -49,6 +59,43 @@ void VoltagePattern::EvaluatePolynomial(const UVector<double>& coefficients,
}
}
void VoltagePattern::EvaluateAiryDisk(double dish_diameter_in_m,
double blocked_diameter_in_m) {
// Number of evaluated grid points in the 1D radial function
constexpr size_t n_samples = 10000;
values_.resize(n_samples);
inverse_increment_radius_ = double(n_samples - 1) / maximum_radius_arc_min_;
// This scales the maximum radius from arcminutes on the sky at
// 1 GHz for a 24.5 m unblocked aperture to the J1 Bessel function
// coordinates (7.016 at the 2nd null).
const double dimensionless_max_in_rad = maximum_radius_arc_min_ * 7.016 /
(1.566 * 60.) * dish_diameter_in_m /
24.5;
const double dimensionless_inverse_in_rad =
static_cast<double>(n_samples - 1) / dimensionless_max_in_rad;
if (blocked_diameter_in_m == 0.0) {
values_[0] = 1.0;
for (size_t i = 1; i < n_samples; i++) {
const double x = static_cast<double>(i) / dimensionless_inverse_in_rad;
values_[i] = 2.0 * std::cyl_bessel_j(1, x) / x;
}
} else {
const double length_ratio = dish_diameter_in_m / blocked_diameter_in_m;
const double area_ratio = length_ratio * length_ratio;
const double area_norm = area_ratio - 1.0;
values_[0] = 1.0;
for (size_t i = 1; i < n_samples; i++) {
const double x = static_cast<double>(i) / dimensionless_inverse_in_rad;
values_[i] =
(area_ratio * 2.0 * std::cyl_bessel_j(1, x) / x -
2.0 * std::cyl_bessel_j(1, x * length_ratio) / (x * length_ratio)) /
area_norm;
}
}
}
UVector<double> VoltagePattern::InterpolateValues(double freq) const {
UVector<double> result;
size_t ifit = 0;
......@@ -84,20 +131,13 @@ const double* VoltagePattern::InterpolateValues(
}
}
double VoltagePattern::LmMaxSquared(double frequency_hz) const {
const double factor =
(180.0 / M_PI) * 60.0 * frequency_hz * 1e-9; // arcminutes * GHz
const double rmax = maximum_radius_arc_min_ / factor;
return rmax * rmax;
}
void VoltagePattern::Render(std::complex<float>* aterm, size_t width,
size_t height, double pixel_scale_x,
double pixel_scale_y, double phase_centre_ra,
double phase_centre_dec, double pointing_ra,
double pointing_dec, double phase_centre_dl,
double phase_centre_dm, double frequency_hz) const {
const double lm_max_sq = LmMaxSquared(frequency_hz);
const double lm_max_sq = LmMaxSquared(frequency_hz, maximum_radius_arc_min_);
UVector<double> interpolated_values;
const double* vp = InterpolateValues(frequency_hz, interpolated_values);
......@@ -144,7 +184,7 @@ void VoltagePattern::Render(std::complex<float>* aterm, size_t width,
void VoltagePattern::Render(std::complex<float>* aterm, double phase_centre_ra,
double phase_centre_dec, double pointing_ra,
double pointing_dec, double frequency_hz) const {
const double lm_max_sq = LmMaxSquared(frequency_hz);
const double lm_max_sq = LmMaxSquared(frequency_hz, maximum_radius_arc_min_);
UVector<double> interpolated_values;
const double* vp = InterpolateValues(frequency_hz, interpolated_values);
......
......@@ -10,13 +10,12 @@
namespace everybeam {
namespace circularsymmetric {
//! Holds the information for a symmetric voltage pattern
//! Holds the information for a circularly symmetric voltage pattern
class VoltagePattern {
public:
VoltagePattern(aocommon::UVector<double> frequencies,
double maximum_radius_arc_min, double reference_frequency)
double maximum_radius_arc_min)
: maximum_radius_arc_min_(maximum_radius_arc_min),
reference_frequency_(reference_frequency),
frequencies_(std::move(frequencies)){};
size_t NSamples() const { return values_.size() / frequencies_.size(); }
......@@ -26,15 +25,23 @@ class VoltagePattern {
}
void EvaluatePolynomial(const aocommon::UVector<double>& coefficients,
bool invert);
double reference_frequency, bool invert);
void EvaluateAiryDisk(double dish_diameter_in_m,
double blocked_diameter_in_m);
/**
* Interpolate response onto a 2D grid.
*/
void Render(std::complex<float>* aterm, size_t width, size_t height,
double pixel_scale_x, double pixel_scale_y,
double phase_centre_ra, double phase_centre_dec,
double pointing_ra, double pointing_dec, double phase_centre_dl,
double phase_centre_dm, double frequency_hz) const;
// Specialization for single point
/**
* Interpolate response for a single point.
*/
void Render(std::complex<float>* aterm, double phase_centre_ra,
double phase_centre_dec, double pointing_ra, double pointing_dec,
double frequency_hz) const;
......@@ -47,11 +54,8 @@ class VoltagePattern {
double frequency_hz,
aocommon::UVector<double>& interpolated_values) const;
double LmMaxSquared(double frequency_hz) const;
double maximum_radius_arc_min_;
double inverse_increment_radius_;
const double maximum_radius_arc_min_;
const double reference_frequency_;
// These are the radial (one-dimensional) values of the beam
// It is an array of size nsamples x nfrequencies, where the sample index is
......
......@@ -4,3 +4,12 @@
install(FILES griddedresponse.h aartfaacgrid.h dishgrid.h mwagrid.h
phasedarraygrid.h skamidgrid.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/griddedresponse")
install(
FILES aartfaacgrid.h
airygrid.h
dishgrid.h
griddedresponse.h
mwagrid.h
phasedarraygrid.h
skamidgrid.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/griddedresponse")
// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
#include "airygrid.h"
#include "../telescope/alma.h"
#include "../circularsymmetric/voltagepattern.h"
using aocommon::HMC4x4;
using aocommon::UVector;
namespace everybeam {
namespace griddedresponse {
void AiryGrid::Response([[maybe_unused]] BeamMode beam_mode,
std::complex<float>* buffer,
[[maybe_unused]] double time, double frequency,
size_t station_idx, size_t field_id) {
// While the AiryGrid class could be made generic, since the Alma
// telescope is currently the only telescope making use of this Airy disc
// code, we use the specific Alma object here.
// TODO merge this implementation with the SkaMidTelescope class, which
// also implements an Airy disc (see also AST-1376).
const telescope::Alma& alma_telescope =
static_cast<const telescope::Alma&>(*telescope_);
double pdir_ra;
double pdir_dec;
std::tie(pdir_ra, pdir_dec) = alma_telescope.GetFieldPointing()[field_id];
const telescope::AiryParameters parameters =
alma_telescope.GetAiryParameters(station_idx);
circularsymmetric::VoltagePattern vp({frequency},
parameters.maximum_radius_arc_min);
vp.EvaluateAiryDisk(parameters.dish_diameter_in_m,
parameters.blocked_diameter_in_m);
vp.Render(buffer, width_, height_, dl_, dm_, ra_, dec_, pdir_ra, pdir_dec,
l_shift_, m_shift_, frequency);
}
void AiryGrid::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer, double time,
double frequency, size_t field_id) {
const telescope::Alma& alma_telescope =
static_cast<const telescope::Alma&>(*telescope_);
if (alma_telescope.IsHomogeneous())
HomogeneousAllStationResponse(beam_mode, buffer, time, frequency, field_id);
else
InhomogeneousAllStationResponse(beam_mode, buffer, time, frequency,
field_id);
}
} // namespace griddedresponse
} // namespace everybeam
// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
#ifndef EVERYBEAM_GRIDDEDRESPONSE_AIRYGRID_H_
#define EVERYBEAM_GRIDDEDRESPONSE_AIRYGRID_H_
#include "griddedresponse.h"
namespace everybeam {
namespace griddedresponse {
/**
* @brief Class for computing the gridded response of (possibly heterogenous)
* telescopes with an Airy disk response, e.g. ALMA.
*/
class [[gnu::visibility("default")]] AiryGrid final : public GriddedResponse {
public:
AiryGrid(const telescope::Telescope* telescope_ptr,
const aocommon::CoordinateSystem coordinate_system)
: GriddedResponse(telescope_ptr, coordinate_system){};
void Response(BeamMode beam_mode, std::complex<float> * buffer, double time,
double frequency, size_t station_idx, size_t field_id) override;
void ResponseAllStations(BeamMode beam_mode, std::complex<float> * buffer,
double time, double frequency, size_t field_id)
override;
void IntegratedResponse(
BeamMode beam_mode, float* buffer,
[[maybe_unused]] const std::vector<double>& time_array, double frequency,
size_t field_id, size_t undersampling_factor,
[[maybe_unused]] const std::vector<double>& baseline_weights) override {
// The integrated response of an Airy disk telescope is independent of time,
// so call IntegratedResponse as if it were one time step
GriddedResponse::IntegratedResponse(beam_mode, buffer, 0.0, frequency,
field_id, undersampling_factor, {0.0});
}
bool PerformUndersampling() const override { return false; }
private:
/**
* @brief Make integrated snapshot, specialized/simplified for dish
* telescopes.
*
* @param matrices Vector of Mueller matrices
* @param frequency Frequency (Hz)
* @param field_id Field id
*/
void MakeIntegratedDishSnapshot(std::vector<aocommon::HMC4x4> & matrices,
double frequency, size_t field_id);
};
} // namespace griddedresponse
} // namespace everybeam
#endif // EVERYBEAM_GRIDDEDRESPONSE_DISHGRID_H_
......@@ -4,7 +4,6 @@
#include "dishgrid.h"
#include "../telescope/dish.h"
#include "../circularsymmetric/voltagepattern.h"
#include "../circularsymmetric/vlacoefficients.h"
#include <aocommon/uvector.h>
#include <aocommon/matrix2x2.h>
......@@ -34,26 +33,14 @@ void DishGrid::Response([[maybe_unused]] BeamMode beam_mode,
dish_telescope.GetDishCoefficients()->ReferenceFrequency();
circularsymmetric::VoltagePattern vp(
dish_telescope.GetDishCoefficients()->GetFrequencies(frequency),
max_radius_arc_min, reference_frequency);
max_radius_arc_min);
const aocommon::UVector<double> coefs_vec =
dish_telescope.GetDishCoefficients()->GetCoefficients(frequency);
vp.EvaluatePolynomial(coefs_vec, false);
vp.EvaluatePolynomial(coefs_vec, reference_frequency, false);
vp.Render(buffer, width_, height_, dl_, dm_, ra_, dec_, pdir_ra, pdir_dec,
l_shift_, m_shift_, frequency);
}
void DishGrid::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer, double,
double frequency, size_t field_id) {
Response(beam_mode, buffer, 0.0, frequency, 0, field_id);
const size_t station_buffer = width_ * height_ * 4;
// Just repeat nstations times
for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
std::copy_n(buffer, station_buffer, buffer + i * station_buffer);
}
}
void DishGrid::IntegratedResponse(
BeamMode /* beam_mode */, float* buffer, double, double frequency,
size_t field_id, size_t undersampling_factor,
......
......@@ -12,9 +12,8 @@ namespace everybeam {
namespace griddedresponse {
/**
* @brief Class for computing the gridded response of dish telescopes,
* e.g. VLA, ATCA, GMRT, SKA-MID
*
* @brief Class for computing the gridded response of homogeneous dish
* telescopes, e.g. VLA, ATCA, GMRT, MeerKAT, SKA-MID.
*/
class [[gnu::visibility("default")]] DishGrid : public GriddedResponse {
public:
......@@ -27,7 +26,9 @@ class [[gnu::visibility("default")]] DishGrid : public GriddedResponse {
void ResponseAllStations(BeamMode beam_mode, std::complex<float> * buffer,
double time, double frequency, size_t field_id)
final override;
final override {
HomogeneousAllStationResponse(beam_mode, buffer, time, frequency, field_id);
}
void IntegratedResponse(
BeamMode beam_mode, float* buffer, double time, double frequency,
......
......@@ -200,6 +200,30 @@ void GriddedResponse::MakeIntegratedSnapshot(
}
}
void GriddedResponse::HomogeneousAllStationResponse(BeamMode beam_mode,
std::complex<float>* buffer,
double time,
double frequency,
size_t field_id) {
Response(beam_mode, buffer, time, frequency, 0, field_id);
const size_t station_buffer = width_ * height_ * 4;
// Repeat nstations times
for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
std::copy_n(buffer, station_buffer, buffer + i * station_buffer);
}
}
void GriddedResponse::InhomogeneousAllStationResponse(
BeamMode beam_mode, std::complex<float>* buffer, double time,
double frequency, size_t field_id) {
const size_t station_buffer_size = width_ * height_ * 4;
for (size_t i = 0; i != telescope_->GetNrStations(); ++i) {
Response(beam_mode, buffer, time, frequency, i, field_id);
buffer += station_buffer_size;
}
}
void GriddedResponse::DoFFTResampling(
float* destination, int width_in, int height_in, int width_out,
int height_out, const std::vector<aocommon::HMC4x4>& matrices) {
......
......@@ -69,8 +69,8 @@ class [[gnu::visibility("default")]] GriddedResponse {
double frequency, size_t field_id) = 0;
/**
* @brief Calculate integrated/undersampled beam for a single time step.
* This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
* @brief Calculate baseline-integrated (undersampled) beam for a single time
* step. This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
* may override MakeIntegratedSnapshot() to implement a more efficient
* version.
*
......@@ -90,7 +90,7 @@ class [[gnu::visibility("default")]] GriddedResponse {
const std::vector<double>& baseline_weights);
/**
* @brief Calculate integrated beam over multiple time steps.
* @brief Calculate baseline-integrated beam over multiple time steps.
* This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
* may override MakeIntegratedSnapshot() to implement a more efficient
* version. This function stores all Mueller matrices for the full size
......@@ -179,6 +179,29 @@ class [[gnu::visibility("default")]] GriddedResponse {
int width_out, int height_out,
const std::vector<aocommon::HMC4x4>& matrices);
/**
* Fills the buffer by querying the response of the first antenna and copies
* it n_stations times. Telescopes with homogeneous antennas can use this
* function to implement @ref ResponseAllStations().
*
* @param buffer should be of size width * height * 4 * n_stations.
*/
void HomogeneousAllStationResponse(BeamMode beam_mode,
std::complex<float> * buffer, double time,
double frequency, size_t field_id);
/**
* Fills the buffer by querying the response of the antennas one by one.
* Telescopes with heterogeneous antennas can use this function to implement
* @ref ResponseAllStations(). Homogeneous arrays can use the faster
* function @ref HomogeneousAllStationResponse() instead.
*
* @param buffer should be of size width * height * 4 * n_stations.
*/
void InhomogeneousAllStationResponse(
BeamMode beam_mode, std::complex<float> * buffer, double time,
double frequency, size_t field_id);
const telescope::Telescope* telescope_;
size_t width_, height_;
double ra_, dec_, dl_, dm_, l_shift_, m_shift_;
......@@ -187,7 +210,7 @@ class [[gnu::visibility("default")]] GriddedResponse {
/**
* @brief Calculate a baseline-integrated snapshot.
* By default, this function will request the response for all antennas and
* perform a weighted average. (Partly) homogenous arrays may implement a
* perform a weighted average. (Partly) homogeneous arrays may implement a
* faster implementation by overriding this method.
*/
virtual void MakeIntegratedSnapshot(
......
......@@ -60,18 +60,6 @@ void MWAGrid::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
}
}
void MWAGrid::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer, double time,
double frequency,
[[maybe_unused]] size_t field_id) {
Response(beam_mode, buffer, time, frequency, 0, 0);
const size_t station_buffer = width_ * height_ * 4;
// Repeated copy for nstations
for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
std::copy_n(buffer, station_buffer, buffer + i * station_buffer);
}
}
void MWAGrid::MakeIntegratedSnapshot(BeamMode beam_mode,
std::vector<aocommon::HMC4x4>& matrices,
double time, double frequency,
......@@ -105,4 +93,4 @@ void MWAGrid::MakeIntegratedSnapshot(BeamMode beam_mode,
}
}
} // namespace griddedresponse
} // namespace everybeam
\ No newline at end of file
} // namespace everybeam
......@@ -24,7 +24,9 @@ class [[gnu::visibility("default")]] MWAGrid final : public GriddedResponse {
void ResponseAllStations(BeamMode beam_mode, std::complex<float> * buffer,
double time, double frequency, size_t field_id)
override;
override {
HomogeneousAllStationResponse(beam_mode, buffer, time, frequency, field_id);
}
private:
std::unique_ptr<everybeam::mwabeam::TileBeam2016> tile_beam_;
......
......@@ -3,8 +3,9 @@
#include "load.h"
#include "telescope/lofar.h"
#include "telescope/alma.h"
#include "telescope/dish.h"
#include "telescope/lofar.h"
#include "telescope/lwa.h"
#include "telescope/mwa.h"
#include "telescope/oskar.h"
......@@ -34,6 +35,8 @@ TelescopeType GetTelescopeType(const casacore::MeasurementSet& ms) {
return kAARTFAAC;
} else if (telescope_name.compare(0, 4, "ATCA") == 0) {
return kATCATelescope;
} else if (telescope_name == "ALMA") {
return kALMATelescope;
} else if (telescope_name.compare(0, 4, "EVLA") == 0) {
return kVLATelescope;
} else if (telescope_name == "GMRT") {
......@@ -65,6 +68,9 @@ std::unique_ptr<telescope::Telescope> Load(const casacore::MeasurementSet& ms,
case kLofarTelescope:
telescope = std::make_unique<telescope::LOFAR>(ms, options);
break;
case kALMATelescope:
telescope = std::make_unique<telescope::Alma>(ms, options);
break;
case kATCATelescope: {
auto coefs = std::make_unique<circularsymmetric::ATCACoefficients>();
telescope =
......
......@@ -19,6 +19,7 @@ enum TelescopeType {
kUnknownTelescope,
kAARTFAAC,
kATCATelescope,
kALMATelescope,
kGMRTTelescope,
kLofarTelescope,
kMeerKATTelescope,
......
# Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
# SPDX-License-Identifier: GPL-3.0-or-later
install(FILES pointresponse.h phasedarraypoint.h aartfaacpoint.h dishpoint.h
mwapoint.h skamidpoint.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/pointresponse")
install(
FILES aartfaacpoint.h
airypoint.h
dishpoint.h
mwapoint.h
pointresponse.h
phasedarraypoint.h
skamidpoint.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/pointresponse")
// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
#include "airypoint.h"
#include "../telescope/alma.h"
#include "../circularsymmetric/voltagepattern.h"
#include "../circularsymmetric/vlacoefficients.h"
#include <aocommon/uvector.h>
namespace everybeam::pointresponse {
void AiryPoint::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
double ra, double dec, double frequency,
size_t station_idx, size_t field_id) {
const telescope::Alma& alma_telescope =
static_cast<const telescope::Alma&>(GetTelescope());
double pdir_ra;
double pdir_dec;
std::tie(pdir_ra, pdir_dec) = alma_telescope.GetFieldPointing()[field_id];
const telescope::AiryParameters parameters =
alma_telescope.GetAiryParameters(station_idx);
circularsymmetric::VoltagePattern vp({frequency},
parameters.maximum_radius_arc_min);
vp.EvaluateAiryDisk(parameters.dish_diameter_in_m,
parameters.blocked_diameter_in_m);
vp.Render(buffer, ra, dec, pdir_ra, pdir_dec, frequency);
}
void AiryPoint::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer, double ra,
double dec, double freq, size_t field_id) {
const telescope::Alma& alma_telescope =
static_cast<const telescope::Alma&>(GetTelescope());
if (alma_telescope.IsHomogeneous())
HomogeneousAllStationsResponse(beam_mode, buffer, ra, dec, freq, field_id);
else
ResponseAllStations(beam_mode, buffer, ra, dec, freq, field_id);
}
} // namespace everybeam::pointresponse
#ifndef EVERYBEAM_POINTRESPONSE_AIRY_POINT_H_
#define EVERYBEAM_POINTRESPONSE_AIRY_POINT_H_
#include "pointresponse.h"
namespace everybeam {
namespace pointresponse {
/**
* @brief Class for computing the directional response of telescopes with
* an Airy Disc response, e.g. ALMA.
*/
class [[gnu::visibility("default")]] AiryPoint final : public PointResponse {
public:
AiryPoint(const telescope::Telescope* telescope_ptr, double time)
: PointResponse(telescope_ptr, time){};
void Response(BeamMode beam_mode, std::complex<float> * buffer, double ra,
double dec, double freq, size_t station_idx, size_t field_id)
override;
void ResponseAllStations(BeamMode beam_mode, std::complex<float> * buffer,
double ra, double dec, double freq, size_t field_id)
override;
};
} // namespace pointresponse
} // namespace everybeam
#endif // EVERYBEAM_POINTRESPONSE_DISHPOINT_H_
......@@ -25,10 +25,10 @@ void DishPoint::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
dish_telescope.GetDishCoefficients()->ReferenceFrequency();
circularsymmetric::VoltagePattern vp(
dish_telescope.GetDishCoefficients()->GetFrequencies(freq),
max_radius_arc_min, reference_frequency);
max_radius_arc_min);
const aocommon::UVector<double> coefs_vec =
dish_telescope.GetDishCoefficients()->GetCoefficients(freq);
vp.EvaluatePolynomial(coefs_vec,
vp.EvaluatePolynomial(coefs_vec, reference_frequency,
dish_telescope.GetDishCoefficients()->AreInverted());
vp.Render(buffer, ra, dec, pdir_ra, pdir_dec, freq);
}
......
......@@ -7,6 +7,7 @@
#ifndef EVERYBEAM_POINTRESPONSE_POINTRESPONSE_H_
#define EVERYBEAM_POINTRESPONSE_POINTRESPONSE_H_
#include <algorithm>
#include <complex>
#include <mutex>
......@@ -147,6 +148,18 @@ class PointResponse {
void ClearTimeUpdate() { has_time_update_ = false; }
void HomogeneousAllStationsResponse(BeamMode beam_mode,
std::complex<float>* buffer, double ra,
double dec, double freq,
size_t field_id) {
Response(beam_mode, buffer, ra, dec, freq, 0, field_id);
// Repeat the same response for all other stations
for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
std::copy_n(buffer, 4, buffer + i * 4);
}
}
private:
const telescope::Telescope* telescope_;
double time_;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment