// Station.cc: Representation of the station beam former.
//
// Copyright (C) 2013
// ASTRON (Netherlands Institute for Radio Astronomy)
// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
//
// This file is part of the LOFAR software suite.
// The LOFAR software suite is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// The LOFAR software suite is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
//
// $Id$
#include "station.h"
#include "common/mathutils.h"
#include "hamaker/hamakerelementresponse.h"
#include "oskar/oskarelementresponse.h"
#include "lobes/lobeselementresponse.h"
using namespace everybeam;
Station::Station(const std::string &name, const vector3r_t &position,
const ElementResponseModel model)
: name_(name),
position_(position),
phase_reference_(position),
element_response_(nullptr) {
SetResponseModel(model);
vector3r_t ncp = {{0.0, 0.0, 1.0}};
ncp_.reset(new coords::ITRFDirection(ncp));
vector3r_t ncppol0 = {{1.0, 0.0, 0.0}};
ncp_pol0_.reset(new coords::ITRFDirection(ncppol0));
}
void Station::SetResponseModel(const ElementResponseModel model) {
switch (model) {
case kHamaker:
element_response_.set(HamakerElementResponse::GetInstance(name_));
break;
case kOSKARDipole:
element_response_.set(OSKARElementResponseDipole::GetInstance());
break;
case kOSKARSphericalWave:
element_response_.set(OSKARElementResponseSphericalWave::GetInstance());
break;
case kLOBES:
element_response_.set(LOBESElementResponse::GetInstance(name_));
break;
default:
std::stringstream message;
message << "The requested element response model '" << model
<< "' is not implemented.";
throw std::runtime_error(message.str());
}
}
void Station::SetResponse(std::shared_ptr<ElementResponse> element_response) {
element_response_.set(element_response);
}
const std::string &Station::GetName() const { return name_; }
const vector3r_t &Station::GetPosition() const { return position_; }
void Station::SetPhaseReference(const vector3r_t &reference) {
phase_reference_ = reference;
}
const vector3r_t &Station::GetPhaseReference() const {
return phase_reference_;
}
void Station::SetAntenna(Antenna::Ptr antenna) {
antenna_ = antenna;
// The antenna can be either an Element or a BeamFormer
// If it is a BeamFormer we recursively extract the first antenna
// until we have an Element.
// The extraction returns copies so antenna_ remains unchanged.
// The element that is found is used in ComputeElementResponse to
// compute the element response.
while (auto beamformer = std::dynamic_pointer_cast<BeamFormer>(antenna)) {
antenna = beamformer->ExtractAntenna(0);
}
element_ = std::dynamic_pointer_cast<Element>(antenna);
}
// ========================================================
matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
const vector3r_t &direction,
size_t id,
const bool rotate) const {
Antenna::Options options;
options.rotate = rotate;
if (rotate) {
vector3r_t ncp_ = NCP(time);
vector3r_t east = normalize(cross(ncp_, direction));
vector3r_t north = cross(direction, east);
options.east = east;
options.north = north;
}
return element_->LocalResponse(time, freq, direction, id, options);
}
matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
const vector3r_t &direction,
const bool rotate) const {
Antenna::Options options;
options.rotate = rotate;
if (options.rotate) {
vector3r_t ncp_ = NCP(time);
vector3r_t east = normalize(cross(ncp_, direction));
vector3r_t north = cross(direction, east);
options.east = east;
options.north = north;
}
matrix22c_t response;
response = element_->Response(time, freq, direction, options);
return response;
}
matrix22c_t Station::Response(real_t time, real_t freq,
const vector3r_t &direction, real_t freq0,
const vector3r_t &station0,
const vector3r_t &tile0,
const bool rotate) const {
Antenna::Options options = {
.freq0 = freq0, .station0 = station0, .tile0 = tile0, .rotate = rotate};
if (rotate) {
vector3r_t ncp_ = NCP(time);
vector3r_t east = normalize(cross(ncp_, direction));
vector3r_t north = cross(direction, east);
options.east = east;
options.north = north;
}
matrix22c_t response = antenna_->Response(time, freq, direction, options);
return response;
}
diag22c_t Station::ArrayFactor(real_t time, real_t freq,
const vector3r_t &direction, real_t freq0,
const vector3r_t &station0,
const vector3r_t &tile0) const {
Antenna::Options options = {
.freq0 = freq0, .station0 = station0, .tile0 = tile0};
return antenna_->ArrayFactor(time, freq, direction, options);
}
matrix22r_t Station::Rotation(real_t time, const vector3r_t &direction) const {
// rotation needs to be optional, normally you only want to rotate your
// coordinatesytem for the center of your (mosaiced) image
// Compute the cross product of the NCP and the target direction. This
// yields a vector tangent to the celestial sphere at the target
// direction, pointing towards the East (the direction of +Y in the IAU
// definition, or positive right ascension).
// Test if the direction is equal to the NCP. If it is, take a random
// vector orthogonal to v1 (the east is not defined here).
vector3r_t v1;
if (std::abs(NCP(time)[0] - direction[0]) < 1e-9 &&
std::abs(NCP(time)[1] - direction[1]) < 1e-9 &&
std::abs(NCP(time)[2] - direction[2]) < 1e-9) {
// Make sure v1 is orthogonal to NCP(time). In the direction of the meridian
v1 = normalize(NCPPol0(time));
} else {
v1 = normalize(cross(NCP(time), direction));
}
// Compute the cross product of the antenna field normal (R) and the
// target direction. This yields a vector tangent to the topocentric
// spherical coordinate system at the target direction, pointing towards
// the direction of positive phi (which runs East over North around the
// pseudo zenith).
// Test if the normal is equal to the target direction. If it is, take
// a random vector orthogonal to the normal.
vector3r_t v2;
if (std::abs(antenna_->coordinate_system_.axes.r[0] - direction[0]) < 1e-9 &&
std::abs(antenna_->coordinate_system_.axes.r[1] - direction[1]) < 1e-9 &&
std::abs(antenna_->coordinate_system_.axes.r[2] - direction[2]) < 1e-9) {
// Nothing to be rotated if the direction is equal to zenith
v2 = v1;
} else {
v2 = normalize(cross(antenna_->coordinate_system_.axes.r, direction));
}
// Compute the cosine and sine of the parallactic angle, i.e. the angle
// between v1 and v2, both tangent to a latitude circle of their
// respective spherical coordinate systems.
real_t coschi = dot(v1, v2);
real_t sinchi;
if (coschi == 1.0)
sinchi = 0.0;
else
sinchi = dot(cross(v1, v2), direction);
// The input coordinate system is a right handed system with its third
// axis along the direction of propagation (IAU +Z). The output
// coordinate system is right handed as well, but its third axis points
// in the direction of arrival (i.e. exactly opposite).
//
// Because the electromagnetic field is always perpendicular to the
// direction of propagation, we only need to relate the (X, Y) axes of
// the input system to the corresponding (theta, phi) axes of the output
// system.
//
// To this end, we first rotate the input system around its third axis
// to align the Y axis with the phi axis. The X and theta axis are
// parallel after this rotation, but point in opposite directions. To
// align the X axis with the theta axis, we flip it.
//
// The Jones matrix to align the Y axis with the phi axis when these are
// separated by an angle phi (measured counter-clockwise around the
// direction of propagation, looking towards the origin), is given by:
//
// [ cos(phi) sin(phi)]
// [-sin(phi) cos(phi)]
//
// Here, cos(phi) and sin(phi) can be computed directly, without having
// to compute phi first (see the computation of coschi and sinchi
// above).
//
// Now, sinchi as computed above is opposite to sin(phi), because the
// direction used in the computation is the direction of arrival instead
// of the direction of propagation. Therefore, the sign of sinchi needs
// to be reversed. Furthermore, as explained above, the X axis has to be
// flipped to align with the theta axis. The Jones matrix returned from
// this function is therefore given by:
//
// [-coschi sinchi]
// [ sinchi coschi]
matrix22r_t rotation = {{{{-coschi, sinchi}}, {{sinchi, coschi}}}};
return rotation;
}
vector3r_t Station::NCP(real_t time) const { return ncp_->at(time); }
vector3r_t Station::NCPPol0(real_t time) const { return ncp_pol0_->at(time); }