Simulator.h

Go to the documentation of this file.

// Simulator.h: Compute visibilities for different model components types
// (implementation of ModelComponentVisitor).
//
// Copyright (C) 2023 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef DP3_BASE_SIMULATOR_H_
#define DP3_BASE_SIMULATOR_H_
 
#include <vector>
 
#include <casacore/casa/Arrays/Cube.h>
 
#include <xtensor/containers/xtensor.hpp>
 
#include "Baseline.h"
#include "ModelComponent.h"
#include "ModelComponentVisitor.h"
#include "Direction.h"
 
namespace dp3 {
namespace base {
 
inline void radec2lmn(const Direction& reference, const Direction& direction,
                      double* lmn) {
  const double delta_ra = direction.ra - reference.ra;
  const double sin_delta_ra = std::sin(delta_ra);
  const double cos_delta_ra = std::cos(delta_ra);
  const double sin_dec = std::sin(direction.dec);
  const double cos_dec = std::cos(direction.dec);
  const double sin_dec0 = std::sin(reference.dec);
  const double cos_dec0 = std::cos(reference.dec);
 
  lmn[0] = cos_dec * sin_delta_ra;
  lmn[1] = sin_dec * cos_dec0 - cos_dec * sin_dec0 * cos_delta_ra;
  // Normally, n is calculated using n = std::sqrt(1.0 - l * l - m * m).
  // However the sign of n is lost that way, so a calculation is used that
  // avoids losing the sign. This formula can be found in Perley (1989).
  // Be aware that we asserted that the sign is wrong in Perley (1989),
  // so this is corrected.
  lmn[2] = sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
}
 
 
typedef std::complex<double> dcomplex;
 
class Simulator : public ModelComponentVisitor {
 public:
  Simulator(const Direction& reference, size_t nStation,
            const std::vector<Baseline>& baselines,
            const std::vector<double>& freq,
            const std::vector<double>& chanWidths,
            const xt::xtensor<double, 2>& stationUVW,
            casacore::Cube<dcomplex>& buffer, bool correctFreqSmearing,
            bool stokesIOnly);
 
  Simulator(const Direction& reference, size_t nStation,
            const std::vector<Baseline>& baselines,
            const std::vector<double>& freq,
            const std::vector<double>& chanWidths,
            const std::vector<double>& scaled_ncp_uvw,
            const xt::xtensor<double, 2>& stationUVW,
            casacore::Cube<dcomplex>& buffer, bool correctTimeSmearing,
            bool correctFreqSmearing, bool stokesIOnly);
 
  // Note DuoMatrix is actually two T matrices
  // T: floating point type, ideally float, double, or long double.
  template <typename T>
  class DuoMatrix {
   public:
    DuoMatrix() = default;
 
    DuoMatrix(size_t nrows, size_t ncols) { resize(nrows, ncols); }
 
    void resize(size_t nrows, size_t ncols) {
      itsNRows = nrows;
      itsData_real.resize(nrows * ncols);
      itsData_imag.resize(nrows * ncols);
    }
 
    size_t nRows() { return itsNRows; }
    size_t nCols() { return itsData_real.size() / itsNRows; }
 
    T& real(size_t row, size_t col) {
      return itsData_real[col * itsNRows + row];
    }
    T& imag(size_t row, size_t col) {
      return itsData_imag[col * itsNRows + row];
    }
    T* realdata() { return &itsData_real[0]; }
    T* imagdata() { return &itsData_imag[0]; }
 
   private:
    // Use separate storage for real/imag parts
    std::vector<T> itsData_real;
    std::vector<T> itsData_imag;
    size_t itsNRows{0};
  };
 
  void simulate(const std::shared_ptr<const ModelComponent>& component);
 
 private:
  void visit(const PointSource& component) override;
  void visit(const GaussianSource& component) override;
 
 private:
  Direction itsReference;
  size_t itsNStation, itsNBaseline, itsNChannel;
  bool itsCorrectTimeSmearing;
  bool itsCorrectFreqSmearing;
  bool itsStokesIOnly;
  std::vector<Baseline> itsBaselines;
  std::vector<double> itsFreq;
  std::vector<double> itsChanWidths;
  std::vector<double> itsScaledNcpUvw;
  const xt::xtensor<double, 2>* itsStationUVW;
  casacore::Cube<dcomplex> itsBuffer;
  std::vector<double> itsStationPhases;
  std::vector<double> itsStationEarthRotation;
  DuoMatrix<double> itsShiftBuffer;
  DuoMatrix<double> itsSpectrumBuffer;
};
 
 
}  // namespace base
}  // namespace dp3
 
#endif