AntennaType.cc

Go to the documentation of this file.

#include <rdet/AntennaType.h>
#include <det/VManager.h>
#include <det/Detector.h>

#include <utl/AugerException.h>

#include <cmath>
#include <sstream>

using namespace det;
using namespace std;
using namespace utl;


namespace rdet {

  std::map<std::string, AntennaPattern> AntennaType::fgAntennaPattern;  // global variable that holds the antenna pattern

  AntennaType::AntennaType(const string& antennaType)
  {
    fAntennaType = antennaType;
  }


  void
  AntennaType::BufferAntennaPattern()
  {
    if (fgAntennaPattern.find(fAntennaType) == fgAntennaPattern.end()) {  // check if AntennaPattern has been buffered before
      AntennaPattern antennaPattern;
      std::map<ResponseKey, VectorEffectiveLength> antennaResponse;
      // User Information which a lot is going to appear on the screen:
      ostringstream message;
      message << "Buffering Antenna Pattern for " << fAntennaType;
      INFO(message);

      utl::Validated<vector<string>> frequenciesS;          //fetch the frequencies as string
      GetAntennaData(frequenciesS, "frequency");  //cause they are used as identifier

      utl::Validated<vector<double>> frequencies;          //fetch the frequencies as double
      GetAntennaData(frequencies, "frequency");  //-- don't want to cheat the unit system ;)
      antennaPattern.frequency_lower_bound = frequencies.Get().front();
      antennaPattern.frequency_upper_bound = frequencies.Get().back();
      antennaPattern.frequencies = frequencies.Get();
      antennaPattern.nFrequency = antennaPattern.frequencies.size();
      message.str("");
      message << "buffering frequencies:";
      for (std::vector<double>::iterator it = antennaPattern.frequencies.begin();
           it != antennaPattern.frequencies.end(); ++it) {
        message << ' ' << *it / utl::MHz;
      }
      message << " MHz";
      INFO(message);

      // import zenith angles
      /* the import of the zenith angles is a complicated because in the antennamodel xml file the phi
       * and theta list define pairs. Thus, the theta list contains all theta angles nPhi times.
       * With the following algorithm the number of theta samples is determined and a std::vector
       * containing all theta angles is filled.
       * */
      utl::Validated<vector<double>> vTheta;    //fetch the according theta angles
      GetAntennaData(vTheta, "theta");
      antennaPattern.theta_lower_bound = vTheta.Get().front();
      antennaPattern.theta_upper_bound = vTheta.Get().back();
      unsigned int nTheta = 0;
      bool first_iteration = true;
      double theta_0 = vTheta.Get().at(0);
      for (unsigned int i = 0; i < vTheta.Get().size(); ++i) {
        const double theta = vTheta.Get().at(i);
        if (first_iteration && theta == theta_0 && i) {
          nTheta = i;
          first_iteration = false;
        }
        if (first_iteration) {
          antennaPattern.thetaAngles.push_back(theta);
        } else {  // check with remaining theta angles if theta angles were imported correctly
          if (antennaPattern.thetaAngles.at(i % nTheta) != theta) {
            ERROR("error in importing zenith angles");
            throw utl::AugerException("error in importing zenith angles");
          }
        }
      }
      antennaPattern.nTheta = nTheta;

      message.str("");
      message << "buffering theta angles:";
      for (std::vector<double>::iterator it = antennaPattern.thetaAngles.begin();
           it != antennaPattern.thetaAngles.end(); ++it) {
        message << ' ' << *it / utl::degree;
      }
      message << " degree";
      INFO(message);

      // determine nPhi
      /* The size of the phi and theta list of the antenna-model xml file is the product of the
       * number of theta samples and the number of phi samples. Thus, knowing nTheta, nPhi can be
       * calculated.
       * */
      const int nTotal = vTheta.Get().size();
      antennaPattern.nPhi = nTotal / nTheta;
      if (antennaPattern.nPhi * antennaPattern.nTheta != nTotal) {
        ERROR("error in determining number of azimuth samples");
        throw utl::AugerException("error in determining number of azimuth samples");
      }
      message.str("");
      message << "nPhi has been determined to " << antennaPattern.nPhi
              << "\t nTheta has been determined to " << antennaPattern.nTheta;
      INFO(message);

      // read in phi angles
      utl::Validated<vector<double>> vPhi;    //fetch the phi angles where the spheres
      GetAntennaData(vPhi, "phi");  //have been simulated
      antennaPattern.phi_lower_bound = vPhi.Get().front();
      antennaPattern.phi_upper_bound = vPhi.Get().back();
      for (int iPhi = 0; iPhi < antennaPattern.nPhi; ++iPhi)
        antennaPattern.phiAngles.push_back(vPhi.Get().at(iPhi * antennaPattern.nTheta));

      message.str("");
      message << "buffering phi angles:";
      for (std::vector<double>::iterator it = antennaPattern.phiAngles.begin();
           it != antennaPattern.phiAngles.end(); ++it) {
        message << ' ' << *it / utl::degree;
      }
      message << " degree";
      INFO(message);

      for (unsigned int iFrequency = 0; iFrequency < frequenciesS.Get().size(); ++iFrequency) {
        utl::Validated<vector<double>> sphereEAHTheta_amp;    //...
        utl::Validated<vector<double>> sphereEAHTheta_phase;  //fetching the
        utl::Validated<vector<double>> sphereEAHPhi_amp;      //sphere data with
        utl::Validated<vector<double>> sphereEAHPhi_phase;    //frequency string as id
        GetAntennaData(sphereEAHTheta_amp, "EAHTheta_amp", frequenciesS.Get().at(iFrequency));
        GetAntennaData(sphereEAHTheta_phase, "EAHTheta_phase", frequenciesS.Get().at(iFrequency));
        GetAntennaData(sphereEAHPhi_amp, "EAHPhi_amp", frequenciesS.Get().at(iFrequency));
        GetAntennaData(sphereEAHPhi_phase, "EAHPhi_phase", frequenciesS.Get().at(iFrequency));

        for (int iPhi = 0; iPhi < antennaPattern.nPhi; ++iPhi) {
          const double current_phi = antennaPattern.phiAngles.at(iPhi);
          //message.str("");
          //message << "buffering antenna pattern for phi = " << current_phi / utl::degree;
          //INFO(message.str().c_str());
          for (int iTheta = 0; iTheta < antennaPattern.nTheta; ++iTheta) {
            // check if we are still at the correct phi angle
            //message.str("");
            //message << "theta = " << fThetaAngles.at(iTheta) / utl::degree;
            //message << "\t(index = " << iPhi * fNTheta + iTheta << " phi = "
            //        << phi.Get().at(iPhi * fNTheta + iTheta) << " theta = "
            //        << theta.Get().at(iPhi * fNTheta + iTheta) << ")";
            //INFO(message);
            if (current_phi != vPhi.Get().at(iPhi * antennaPattern.nTheta + iTheta)) {
              throw utl::AugerException("phi angle has changed during theta loop");
            }

            const float Phi_amp = sphereEAHPhi_amp.Get().at(iPhi * antennaPattern.nTheta + iTheta);
            const float Phi_phase = sphereEAHPhi_phase.Get().at(iPhi * antennaPattern.nTheta + iTheta);
            const float Theta_amp = sphereEAHTheta_amp.Get().at(iPhi * antennaPattern.nTheta + iTheta);
            const float Theta_phase = sphereEAHTheta_phase.Get().at(iPhi * antennaPattern.nTheta + iTheta);

            VectorEffectiveLength vel;
            vel.Phi = std::polar(Phi_amp, Phi_phase);
            vel.Theta = std::polar(Theta_amp, Theta_phase);
            ResponseKey key(iFrequency, iTheta, iPhi);
            //cout << "(" << iFrequency << ", " << iTheta << ", " << iPhi << ") = " << VEL.Phi_amp
            //     << ", " << VEL.Phi_phase << ", " << VEL.Theta_amp << ", " << VEL.Theta_phase << endl;
            antennaResponse[key] = vel;
          }
        }
      }
      antennaPattern.AntennaResponse = antennaResponse;

      // calculate integrated antenna response
      for (unsigned int iFrequency = 0; iFrequency < frequenciesS.Get().size(); ++iFrequency) {
        const double integratedVEL =
          CalculateIntegratedEffectiveAntennaHeight(antennaResponse, iFrequency, antennaPattern.nTheta, antennaPattern.nPhi);
        antennaPattern.meanTransfer.push_back(integratedVEL);
      }

      //utl::Validated<vector<double>> MeanTransfer;  // fetch the mean transfer
      //GetAntennaData(MeanTransfer, "MeanTransfer");
      //antennaPattern.meanTransfer = MeanTransfer.Get();

      fgAntennaPattern[fAntennaType] = antennaPattern;
      message.str("");
      message << "finished buffering antenna " << fAntennaType;
      INFO(message);
    }
  }


  std::pair<std::complex<double>, std::complex<double>>
  AntennaType::GetElectricFieldResponse(const double theta, const double phi, const double freq, const std::string& interpolationMode)
  {
    //The antenna patterns go from 0 -> 360 deg in azimuth. Make sure we are in that range
    double matchedPhi = phi;
    while (matchedPhi < 0)
      matchedPhi += 360. * degree;
    while (matchedPhi >= 360. * degree)
      matchedPhi -= 360. * degree;

    if (interpolationMode == "Lookup") {
      return GetElectricFieldResponse_lookup(theta, matchedPhi, freq);
    } else if (interpolationMode == "LinearInterpolation") {
      return GetElectricFieldResponse_LinearInterpolation(theta, matchedPhi, freq);
    } else {
      ostringstream message;
      message << "interpolation mode " << interpolationMode << " unknown";
      throw utl::AugerException(message.str().c_str());
    }
  }


  std::pair<std::complex<double>, std::complex<double>>
  AntennaType::GetComplexRepresentationOfVectorEffectiveLength(const VectorEffectiveLength& vel)
  {
    return std::pair<std::complex<double>, std::complex<double>>(vel.Theta, vel.Phi);
  }


  double
  AntennaType::InterpolateLinear(const double x, const double x_low, const double x_up,
                                 const double y_low, const double y_up)
  {
    if (x_low == x_up)
      return y_low;

    const double result = y_low + (y_up - y_low) * (x - x_low) / (x_up - x_low);
    if (std::isnan(result)) {
      ostringstream msg;
      msg << "linear interpolation results in NaN:\n"
             "\tx = " << x << "\n"
             "\tx_low = " << x_low << "\n"
             "\tx_up = " << x_up << "\n"
             "\ty_low = " << y_low << "\n"
             "\ty_up = " << y_up << "\n"
             "\tresult = " << result << '\n';
      throw AugerException(msg.str());
    }
    return result;
  }


  std::complex<float>
  AntennaType::InterpolateLinear(const float x, const float x_low, const float x_up,
                                 const std::complex<float>& y_low, const std::complex<float>& y_up)
  {
    if (x_low == x_up)
      return y_low;

    const std::complex<float> result = y_low + (y_up - y_low) * ((x - x_low) / (x_up - x_low));

    if (std::isnan(std::abs(result))) {
      ostringstream msg;
      msg << "linear interpolation results in NaN:\n"
             "\tx = " << x << "\n"
             "\tx_low = " << x_low << "\n"
             "\tx_up = " << x_up << "\n"
             "\ty_low = " << y_low << "\n"
             "\ty_up = " << y_up << "\n"
             "\tresult = " << result << '\n';
      throw AugerException(msg.str());
    }

    return result;
  }


  VectorEffectiveLength
  AntennaType::InterpolateLinear(const double x, const double x_low, const double x_up,
                                 const VectorEffectiveLength& vel_low, const VectorEffectiveLength& vel_up)
  {
    const auto theta = InterpolateLinear(x, x_low, x_up, vel_low.Theta, vel_up.Theta);
    const auto phi = InterpolateLinear(x, x_low, x_up, vel_low.Phi, vel_up.Phi);
    return VectorEffectiveLength{theta, phi};
  }


  std::pair<std::complex<double>, std::complex<double>>
  AntennaType::GetElectricFieldResponse_lookup(const double theta, const double phi, const double freq)
  {
    BufferAntennaPattern();
    const double theta_lower_bound = fgAntennaPattern[fAntennaType].theta_lower_bound;
    const double theta_upper_bound = fgAntennaPattern[fAntennaType].theta_upper_bound;
    const double phi_lower_bound = fgAntennaPattern[fAntennaType].phi_lower_bound;
    const double phi_upper_bound = fgAntennaPattern[fAntennaType].phi_upper_bound;
    const double frequency_lower_bound = fgAntennaPattern[fAntennaType].frequency_lower_bound;
    const double frequency_upper_bound = fgAntennaPattern[fAntennaType].frequency_upper_bound;
    const int nTheta = fgAntennaPattern[fAntennaType].nTheta;
    const int nPhi = fgAntennaPattern[fAntennaType].nPhi;
    const int nFrequency = fgAntennaPattern[fAntennaType].nFrequency;

    const int iTheta =
      int(round((theta - theta_lower_bound) / (theta_upper_bound - theta_lower_bound) * (nTheta - 1)));
    const int iPhi = int(round((phi - phi_lower_bound) / (phi_upper_bound - phi_lower_bound) * (nPhi - 1)));
    const int iFrequency =
      int(round( (freq - frequency_lower_bound) / (frequency_upper_bound - frequency_lower_bound) * (nFrequency - 1)));
    ResponseKey key(iFrequency, iTheta, iPhi);
    VectorEffectiveLength& vel = fgAntennaPattern[fAntennaType].AntennaResponse.at(key);
    //cout << "f=" << freq / utl::MHz << "MHz theta=" << theta / utl::deg << " phi=" << phi / utl::deg
    //     << " :" << VEL.Phi_amp << ", " << VEL.Phi_phase / utl::deg << ", " << VEL.Theta_amp << ", "
    //     << VEL.Theta_phase / utl::deg << endl;
    return GetComplexRepresentationOfVectorEffectiveLength(vel);
  }


  std::pair<std::complex<double>, std::complex<double>>
  AntennaType::GetElectricFieldResponse_LinearInterpolation(const double theta, const double phi, const double freq)
  {
    BufferAntennaPattern();

    const double theta_lower_bound = fgAntennaPattern[fAntennaType].theta_lower_bound;
    const double theta_upper_bound = fgAntennaPattern[fAntennaType].theta_upper_bound;
    const double phi_lower_bound = fgAntennaPattern[fAntennaType].phi_lower_bound;
    const double phi_upper_bound = fgAntennaPattern[fAntennaType].phi_upper_bound;
    const double frequency_lower_bound = fgAntennaPattern[fAntennaType].frequency_lower_bound;
    const double frequency_upper_bound = fgAntennaPattern[fAntennaType].frequency_upper_bound;

    if ((freq < frequency_lower_bound || freq > frequency_upper_bound) ||
        (phi < phi_lower_bound || phi > phi_upper_bound) ||
        (theta < theta_lower_bound || theta > theta_upper_bound)) {
      ostringstream msg;
      msg << "frequency (" << freq / utl::MHz << "MHz) or zenith (" << theta / utl::degree
          << "deg) or azimuth (" << phi / utl::degree
          << "deg) angle requested that is not contained in the antenna pattern, returning 0";
      WARNING(msg);
      return std::pair<std::complex<double>, std::complex<double>>(std::complex<double>(0, 0),
                                                                   std::complex<double>(0, 0));
    }

    const int nTheta = fgAntennaPattern[fAntennaType].nTheta;
    const int nPhi = fgAntennaPattern[fAntennaType].nPhi;
    const int nFrequency = fgAntennaPattern[fAntennaType].nFrequency;
    const std::vector<double>& frequencies = fgAntennaPattern[fAntennaType].frequencies;
    const std::vector<double>& thetaAngles = fgAntennaPattern[fAntennaType].thetaAngles;
    const std::vector<double>& phiAngles = fgAntennaPattern[fAntennaType].phiAngles;

    const int iTheta_lower =
      int(floor((theta - theta_lower_bound) / (theta_upper_bound - theta_lower_bound) * (nTheta - 1)));
    const double theta_lower = thetaAngles.at(iTheta_lower);
    const int iTheta_upper =
      int(ceil((theta - theta_lower_bound) / (theta_upper_bound - theta_lower_bound) * (nTheta - 1)));
    const double theta_upper = thetaAngles.at(iTheta_upper);
    const int iPhi_lower =
      int(floor((phi - phi_lower_bound) / (phi_upper_bound - phi_lower_bound) * (nPhi - 1)));
    const double phi_lower = phiAngles.at(iPhi_lower);
    const int iPhi_upper =
      int(ceil((phi - phi_lower_bound) / (phi_upper_bound - phi_lower_bound) * (nPhi - 1)));
    const double phi_upper = phiAngles.at(iPhi_upper);
    const int iFrequency_lower =
      int(floor((freq - frequency_lower_bound) / (frequency_upper_bound - frequency_lower_bound) * (nFrequency - 1)));
    const double frequency_lower = frequencies.at(iFrequency_lower);
    const int iFrequency_upper =
      int( ceil((freq - frequency_lower_bound) / (frequency_upper_bound - frequency_lower_bound) * (nFrequency - 1)));
    const double frequency_upper = frequencies.at(iFrequency_upper);
    const std::map<ResponseKey, VectorEffectiveLength>& antennaResponse =
      fgAntennaPattern[fAntennaType].AntennaResponse;

    // lower frequency bound
    // theta low
    const VectorEffectiveLength vel_freq_low_theta_low =
      InterpolateLinear(
        phi, phi_lower, phi_upper,
        antennaResponse.at(ResponseKey(iFrequency_lower, iTheta_lower, iPhi_lower)),
        antennaResponse.at(ResponseKey(iFrequency_lower, iTheta_lower, iPhi_upper))
      );

    // theta up
    const VectorEffectiveLength vel_freq_low_theta_up =
      InterpolateLinear(
        phi, phi_lower, phi_upper,
        antennaResponse.at(ResponseKey(iFrequency_lower, iTheta_upper, iPhi_lower)),
        antennaResponse.at(ResponseKey(iFrequency_lower, iTheta_upper, iPhi_upper))
      );

    const VectorEffectiveLength vel_freq_low =
      InterpolateLinear(theta, theta_lower, theta_upper, vel_freq_low_theta_low, vel_freq_low_theta_up);

    // upper frequency bound
    // theta low
    const VectorEffectiveLength vel_freq_up_theta_low =
      InterpolateLinear(
        phi, phi_lower, phi_upper,
        antennaResponse.at(ResponseKey(iFrequency_upper, iTheta_lower, iPhi_lower)),
        antennaResponse.at(ResponseKey(iFrequency_upper, iTheta_lower, iPhi_upper))
      );

    // theta up
    const VectorEffectiveLength vel_freq_up_theta_up =
      InterpolateLinear(
        phi, phi_lower, phi_upper,
        antennaResponse.at(ResponseKey(iFrequency_upper, iTheta_upper, iPhi_lower)),
        antennaResponse.at(ResponseKey(iFrequency_upper, iTheta_upper, iPhi_upper))
      );

    const VectorEffectiveLength vel_freq_up =
      InterpolateLinear(theta, theta_lower, theta_upper, vel_freq_up_theta_low, vel_freq_up_theta_up);

    const VectorEffectiveLength vel =
      InterpolateLinear(freq, frequency_lower, frequency_upper, vel_freq_low, vel_freq_up);
    //cout << "f=" << freq / utl::MHz << "MHz theta=" << theta / utl::deg << " phi=" << phi / utl::deg
    //     << " :" << VEL.Phi_amp << ", " << VEL.Phi_phase / utl::deg << ", " << VEL.Theta_amp << ", "
    //     << VEL.Theta_phase / utl::deg << endl;
    return GetComplexRepresentationOfVectorEffectiveLength(vel);
  }


  double
  AntennaType::GetIntegratedEffectiveAntennaHeight(const double freq)
  {
    BufferAntennaPattern();
    const double frequency_lower_bound = fgAntennaPattern[fAntennaType].frequency_lower_bound;
    const double frequency_upper_bound = fgAntennaPattern[fAntennaType].frequency_upper_bound;
    //ostringstream msg;
    if (freq < frequency_lower_bound || freq > frequency_upper_bound) {
      //msg << "frequency (" << freq / utl::MHz
      //    << "MHz) requested that is not contained in the antenna pattern, returning 0";
      //WARNING(msg);
      return 0;
    }

    const int nFrequency = fgAntennaPattern[fAntennaType].nFrequency;
    const std::vector<double>& frequencies = fgAntennaPattern[fAntennaType].frequencies;

    const int iFrequency_lower =
      int(floor((freq - frequency_lower_bound) / (frequency_upper_bound - frequency_lower_bound) * (nFrequency - 1)));
    const double frequency_lower = frequencies.at(iFrequency_lower);
    const int iFrequency_upper =
      int(ceil((freq - frequency_lower_bound) / (frequency_upper_bound - frequency_lower_bound) * (nFrequency - 1)));
    const double frequency_upper = frequencies.at(iFrequency_upper);

    //msg.str("");
    //msg << "requesting frequency " << freq / utl::MHz << " MHz -> frequency bin " << iFrequency_lower
    //    << " and " << iFrequency_upper;
    //INFO(msg.str());

    return InterpolateLinear(freq, frequency_lower, frequency_upper,
                             fgAntennaPattern[fAntennaType].meanTransfer.at(iFrequency_lower),
                             fgAntennaPattern[fAntennaType].meanTransfer.at(iFrequency_upper));
  }


  double
  AntennaType::CalculateIntegratedEffectiveAntennaHeight(const std::map<ResponseKey, VectorEffectiveLength>& antennaResponse,
                                                         const int iFreq, const int nTheta, const int nPhi)
  {
    double integratedVEL = 0;
    for (int iTheta = 0; iTheta < nTheta; ++iTheta) {
      for (int iPhi = 0; iPhi < nPhi; ++iPhi) {
        VectorEffectiveLength vel = antennaResponse.at(ResponseKey(iFreq, iTheta, iPhi));
        integratedVEL += sqrt(std::norm(vel.Phi) + std::norm(vel.Theta));
      }
    }
    integratedVEL *= 2 * M_PI / (nTheta * nPhi * sqrt(2));
    return integratedVEL;
  }


  void
  AntennaType::NotFoundAndExit(const string& msg)
    const
  {
    ostringstream err;
    err << "Did not find requested component: " << msg;
    ERROR(err);
    exit(EXIT_FAILURE);
  }
}