Skip to content
Snippets Groups Projects
Select Git revision
  • gec-84-casacore-3.7.1-in-wheels
  • master default protected
  • extends_cmake_to_use_fetchcontent
  • rtd-ubuntu24-plucky
  • ast-1644-temp-fix
  • ast-1600-fix-beam-for-meerkat-ska-mid
  • readthedocs-c++17
  • fix_oskar
  • 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
  • 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
40 results

Station.cc

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    Station.cc 9.09 KiB
    // 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/MathUtil.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)
        : itsName(name),
          itsPosition(position),
          itsPhaseReference(position),
          itsElementResponse(nullptr) {
      setModel(model);
      vector3r_t ncp = {{0.0, 0.0, 1.0}};
      itsNCP.reset(new coords::ITRFDirection(ncp));
      vector3r_t ncppol0 = {{1.0, 0.0, 0.0}};
      itsNCPPol0.reset(new coords::ITRFDirection(ncppol0));
    }
    
    void Station::setModel(const ElementResponseModel model) {
      switch (model) {
        case Hamaker:
          itsElementResponse.set(HamakerElementResponse::getInstance(itsName));
          break;
        case OSKARDipole:
          itsElementResponse.set(OSKARElementResponseDipole::getInstance());
          break;
        case OSKARSphericalWave:
          itsElementResponse.set(OSKARElementResponseSphericalWave::getInstance());
          break;
        case LOBES:
          itsElementResponse.set(LOBESElementResponse::getInstance(itsName));
          break;
        default:
          std::stringstream message;
          message << "The requested element response model '" << model
                  << "' is not implemented.";
          throw std::runtime_error(message.str());
      }
    }
    
    const std::string &Station::name() const { return itsName; }
    
    const vector3r_t &Station::position() const { return itsPosition; }
    
    void Station::setPhaseReference(const vector3r_t &reference) {
      itsPhaseReference = reference;
    }
    
    const vector3r_t &Station::phaseReference() const { return itsPhaseReference; }
    
    // ========================================================
    matrix22c_t Station::elementResponse(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 itsElement->local_response(time, freq, direction, id, options);
    }
    
    matrix22c_t Station::elementResponse(real_t time, real_t freq,
                                         const vector3r_t &direction,
                                         const bool rotate) const {
      //     if (rotate)
      //       return itsElement->response(time, freq, direction)
      //           * rotation(time, direction);
      //     else
      //       return itsElement->response(time, freq, direction);
    
      return itsElement->response(time, freq, direction);
    }
    
    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 = itsAntenna->response(time, freq, direction, options);
    
      //     if (rotate) {
      //         std::cout << "rotate" << std::endl;
      //         auto r = rotation(time, direction);
      //         std::cout << r[0][0] << ", " << r[0][1] << std::endl;
      //         std::cout << r[1][0] << ", " << r[1][1] << std::endl;
      //         response = response * r;
      //         std::cout << response[0][0] << std::endl;
      //     }
      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 itsAntenna->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(itsAntenna->m_coordinate_system.axes.r[0] - direction[0]) <
              1e-9 &&
          std::abs(itsAntenna->m_coordinate_system.axes.r[1] - direction[1]) <
              1e-9 &&
          std::abs(itsAntenna->m_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(itsAntenna->m_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 itsNCP->at(time); }
    
    vector3r_t Station::ncppol0(real_t time) const { return itsNCPPol0->at(time); }