beam2016implementation.h

#ifndef _GET_FF_H__
#define _GET_FF_H__

/**
 * C++ implementation of Full Embeded Element beam model for MWA based on
 * beam_full_EE.py script and Sokolowski et al (2017) paper Implemented by
 * Marcin Sokolowski (May 2017) - marcin.sokolowski@curtin.edu.au
 */

#include <algorithm>
#include <complex>
#include <cstdlib>
#include <cstring>
#include <string>
#include <vector>
#include <map>
#include <memory>
#include <mutex>

#include <boost/math/special_functions/legendre.hpp>

#include <H5Cpp.h>

#include "factorialtable.h"
#include "recursivelock.h"

namespace everybeam {
namespace mwabeam {

/*
  Structure for Jones matrix :
  j00 j01
  j10 j11
*/
class JonesMatrix {
 public:
  std::complex<double> j00, j01, j10, j11;

  JonesMatrix(double j00_r = 0.00, double j01_r = 0.00, double j10_r = 0.00,
              double j11_r = 0.00)
      : j00(j00_r, 0), j01(j01_r, 0), j10(j10_r, 0), j11(j11_r, 0) {}

  static void zeros(std::vector<std::vector<JonesMatrix>>& jones, size_t x_size,
                    size_t y_size);
};

class Beam2016Implementation {
 public:
  struct DataSetIndex {
    DataSetIndex(char pol_, int antenna_, int freqHz_)
        : pol(pol_), antenna(antenna_), frequencyHz(freqHz_) {}

    bool operator<(const DataSetIndex& rhs) const {
      return std::make_tuple(pol, antenna, frequencyHz) <
             std::make_tuple(rhs.pol, rhs.antenna, rhs.frequencyHz);
    }

    std::string Name() const {
      return pol + std::to_string(antenna + 1) + '_' +
             std::to_string(frequencyHz);
    }

    char pol;
    int antenna, frequencyHz;
  };

  Beam2016Implementation(const double* delays, const double* amps,
                         const std::string& searchPath);

  //--------------------------------------------------------------------
  // Calculation of Jones matrix
  //------------------------------------------------------------------------------------------------------------------
  // Calculate jones matrix in a specified direction for a given frequency,
  // delays and amplitudes. Zenith normalisation can also be disabled - but by
  // default is enabled : This function will  :
  //    - calculate coefficients of spherical waves (SPW) if required (meaning
  //    if frequency, delays or amplitudes are different than last set of
  //    coefficients was calculated for)
  //    - calculate electric fields (equations 4 and 5 in Sokolowski et al
  //    paper)
  //    - normalise Jones matrix to one at zenith at the same frequency (if
  //    required by parameter):
  //
  // Only thse the functions should be used for external calls. The rest is
  // protected from external use.
  //
  // INPUT : (az_deg,za_deg) - azimuth and zenith angles [in degrees]
  //         freq_hz_param   - frequency in Hz
  //         delays          - delays in beamformer steps
  //         amps            - amplitudes
  //         bZenithNorm     - normalise to zenith (>0) or not (<=0)
  // OUTPUT : Jones matrix (normalised or not - depending on the parameter
  // bZenithNorm )
  JonesMatrix CalcJones(double az_deg, double za_deg, int freq_hz_param,
                        bool bZenithNorm = true);
  JonesMatrix CalcJones(double az_deg, double za_deg, int freq_hz_param,
                        const double* delays, const double* amps,
                        RecursiveLock<std::mutex>& lock,
                        bool bZenithNorm = true);

  // Calculation of Jones matrix for an image passed in the arrays azimuth and
  // zenith angles maps. This function calls the single direction one (above)
  // for all pixel in the input image. INPUT :
  //       2D array of azimuth angles in degrees
  //       2D array of zenith angles in degrees
  //       freq_hz_param   - frequency in Hz
  //       delays          - delays in beamformer steps
  //       amps            - amplutudes
  //       bZenithNorm     - normalise to zenith (>0) or not (<=0)
  // OUTPUT :
  //       2D array of JonesMatrix for each pixel
  void CalcJonesArray(std::vector<std::vector<double>>& azim_arr,
                      std::vector<std::vector<double>>& za_arr,
                      std::vector<std::vector<JonesMatrix>>& jones,
                      int freq_hz_param, bool bZenithNorm = true);

 private:
  //---------------------------------------------------- Calculations and
  // variables for spherical waves coefficients (see equations 3-6 in the
  // Sokolowski et al paper)
  // ----------------------------------------------------
  // Coefficients of spherical waves (SPW) - see equations 3-6 in the Sokolowski
  // et al paper 1 refers to s=1 (transverse electric - TE modes) - eq.5 2
  // refers to s=2 (transverse magnetic - TM modes) - eq.5 These coefficients
  // are calculated ones for a given frequency, delays, amplitudes so these
  // parameters are stored to only calculate new coefficients when they change.
  // X polarisation :
  struct Coefficients {
    std::vector<std::complex<double>> Q1_accum;
    std::vector<std::complex<double>> Q2_accum;
    std::vector<double> M_accum;
    std::vector<double> N_accum;
    std::vector<double>
        MabsM;    // precalculated m/abs(m) to make it once for all pointings
    double Nmax;  // maximum N coefficient for Y (=max(N_accum_X)) - to avoid
                  // relaculations
    std::vector<double>
        Cmn;  // coefficient under sumation in equation 3 for X pol.
  } _coefX, _coefY;

  // Calculation of Jones matrix for a single pointing direction (internal
  // function): INPUT : (az_rad,za_rad) - azimuth and zenith in [radians]
  JonesMatrix CalcJonesDirect(double az_rad, double za_rad,
                              const Coefficients& coefsX,
                              const Coefficients& coefsY);

  //--------------------------------------------------------------------
  // Calculation of Jones matrix components for given polarisation (eq. 4 and 5
  // in the Sokolowski et al (2017) paper --------------------------------
  // Internal function to calculate Jones matrix componenets for a given
  // polarisation (pol). The are calculated as electric field vectors E_theta_mn
  // (eq.4) and E_phi_mn (eq.5 in the Sokolowski et al 2017 paper). INPUT :
  //    phi,theta - FEKO convention coordinates (phi=90-azim), theta=za in
  //    radians already Q1_accum, Q2_accum, M_accum, N_accum, MabsM, Nmax -
  //    coefficients calculated earlier in CalcModes function for given
  //    frequency, delays and amplitudes
  // OUTPUT :
  //    Jones matrix filled with components for given polarisation (pol input
  //    parameter)
  void CalcSigmas(double phi, double theta, const Coefficients& coefficients,
                  char pol, JonesMatrix& jones_matrix) const;

  std::complex<double> JPower(size_t i) const { return _jPowerTable[i % 4]; }
  std::complex<double> _jPowerTable[4];

  // Information on last modes parameters - not to recalculate the same again
  // and again !
  int _calcModesLastFreqHz;
  std::vector<double> _calcModesLastDelays;
  std::vector<double> _calcModesLastAmps;

  // function comparing current parameters : frequency, delays and amplitudes
  // with those previously used to calculate spherical waves coefficients
  // (stored in the 3 variables above)
  bool IsCalcModesRequired(int freq_hz, int n_ant, const double* delays,
                           const double* amps);

  // Calculation of modes Q1 and Q2 and coefficients N and M and some derived
  // variables (MabsM_X,MabsM_Y,Nmax_X and Nmax_Y) to make it once for a given
  // pointing and then re-use for many different (azim,za) angles:

  // function calculating coefficients for X and Y and storing parameters
  // frequency, delays and amplitudes
  void GetModes(int freq_hz, size_t n_ant, const double* delays,
                const double* amps, Coefficients& coefsX, Coefficients& coefsY,
                RecursiveLock<std::mutex>& lock);

  // function calculating all coefficients Q1, Q2, N, M and derived MabsM, Nmax
  // for a given polarisation ("X" or "Y") - perhaps enum should be used here
  double CalcModes(int freq_hz, size_t n_ant, const double* delays,
                   const double* amp, char pol, Coefficients& coefs,
                   RecursiveLock<std::mutex>& lock);

  // Calculation of normalisation matrix :
  JonesMatrix CalcZenithNormMatrix(int freq_hz,
                                   RecursiveLock<std::mutex>& lock);

  std::map<int, JonesMatrix> _normJonesCache;

  double _delays[16], _amps[16];

  // ----------------------------------------------------------------
  // Checking frequences included in the H5 file (stored in vector m_freq_list)
  // :
  bool has_freq(int freq_hz) const;
  int find_closest_freq(int freq_hz) const;

  // HDF5 File interface and data structures for H5 data
  // Interface to HDF5 file format and structures to store H5 data
  // Functions for reading H5 file and its datasets :
  // Read data from H5 file, file name is specified in the object constructor
  void Read();

  const std::vector<std::vector<double>>& GetDataSet(
      const DataSetIndex& index, RecursiveLock<std::mutex>& lock);

  // Read dataset_name from H5 file
  void ReadDataSet(const std::string& dataset_name,
                   std::vector<std::vector<double>>& out_vector,
                   H5::H5File& h5File);

  // function for iteration call to H5Ovisit function :
  static herr_t list_obj_iterate(hid_t loc_id, const char* name,
                                 const H5O_info_t* info, void* operator_data);

  std::unique_ptr<H5::H5File> _h5File;
  std::string _searchPath;

  // Data structures for H5 file data :
  std::vector<std::string> m_obj_list;       // list of datasets in H5 file
  std::vector<int> m_freq_list;              // list of simulated frequencies
  std::vector<std::vector<double>> m_Modes;  // data in Modes DataSet

  FactorialTable _factorial;

  // Calculations of Legendre polynomials :
  static void lpmv(std::vector<double>& output, int n, double x);

  // OUTPUT : returns list of Legendre polynomial values calculated up to order
  // nmax :
  static int P1sin(int nmax, double theta, std::vector<double>& p1sin_out,
                   std::vector<double>& p1_out);

  std::map<DataSetIndex, std::vector<std::vector<double>>> _dataSetCache;

  std::mutex mutex_;
};
}  // namespace mwabeam
}  // namespace everybeam
#endif