MeasuredDBMieModel.cc

Go to the documentation of this file.

#include <string>
#include <sstream>
#include <vector>
#include <cmath>
#include <limits>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <atm/ParametricXMLMieModel.h>
#include <atm/ScatteringResult.h>
#include <atm/AttenuationResult.h>
#include <atm/MeasuredDBMieModel.h>

#include <fwk/CentralConfig.h>

#include <utl/Point.h>
#include <utl/UTMPoint.h>
#include <utl/Vector.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/MathConstants.h>
#include <utl/StringCompare.h>

#include <atm/AerosolDB.h>
#include <atm/AerosolZone.h>

using namespace atm;
using namespace fdet;
using namespace utl;
using namespace std;
using namespace fwk;


MeasuredDBMieModel::MeasuredDBMieModel() :
  fVAODVsHeightMap(NULL),
  fDbIsEmptyNow(false),
  fNSigma(0.),
  fSystematicShift(0),
  fZoneParam(eAll)
{
}


void
MeasuredDBMieModel::Init()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("MeasuredDBMieModel");

  Branch zoneB;
  if (bool(zoneB = topB.GetChild("systematicShiftVAOD")))
    zoneB.GetData(fSystematicShift);

  if (bool(zoneB = topB.GetChild("zoneForce"))) {
    string zoneName;
    zoneB.GetData(zoneName);
    fZoneParam = StringToZoneParam(zoneName);
  }
  else if (bool(zoneB = topB.GetChild("zoneSwap"))) {
    vector<string> zoneNames;
    zoneB.GetData(zoneNames);
    if (zoneNames.size() == 2) {
      fSwapZones.clear();
      fSwapZones.push_back(StringToZoneParam(zoneNames[0]));
      fSwapZones.push_back(StringToZoneParam(zoneNames[1]));
    }
    else
      throw InvalidConfigurationException("Zone swap requires two zone names");
  }
}


bool
MeasuredDBMieModel::HasData() const
{
  if (!CheckForUpdates())
    return false;
  else if (fZoneParam == eAll)
    return true;
  else {// zone force
    ZonePosition::iterator zPosIt =
      find_if(fZonePositions.begin(), fZonePositions.end(),
              [&](const auto& zone) { return ZoneMatch(zone, ZoneParamToString(fZoneParam)); });
    return zPosIt != fZonePositions.end();
  }
}


bool
MeasuredDBMieModel::HasZone(const string& zoneName) const
{
  return find_if(fZonePositions.begin(), fZonePositions.end(),
                 [&](const auto& zone) { return ZoneMatch(zone, zoneName); }) != fZonePositions.end();
}


ScatteringResult
MeasuredDBMieModel::EvaluateMieScattering(const utl::Point& xInit,
                                          const utl::Point& xFinal,
                                          const double angle,
                                          const double distance,
                                          const vector<double>& wLength) const
{
  const AttenuationResult mieAttenuation = EvaluateMieAttenuation(xInit, xFinal, wLength);
  return EvaluateMieScattering(xInit, xFinal, angle, distance, mieAttenuation);
}


ScatteringResult
MeasuredDBMieModel::EvaluateMieScattering(const utl::Point& xInit,
                                          const utl::Point& xFinal,
                                          const double angle,
                                          const double distance,
                                          const AttenuationResult& mieAttenuation) const
{
  const utl::TabulatedFunctionErrors& attenuation = mieAttenuation.GetTransmissionFactor();
  utl::TabulatedFunctionErrors frac;
  const double fractionError = 0;
  const double wError = 0;
  // The model implemented here is wavelength independent
  for (TabulatedFunctionErrors::ArrayConstIterator it = attenuation.XBegin();
       it != attenuation.XEnd(); ++it) {
    const double scat = EvaluateMieScattering(xInit, xFinal, angle, distance, *it, attenuation.Y(*it));
    frac.PushBack(*it, wError, scat, fractionError);
  }
  return ScatteringResult(frac, angle);
}


double MeasuredDBMieModel::EvaluateMieScattering(const utl::Point& xInit,
                                                 const utl::Point& xFinal,
                                                 const double angle,
                                                 const double distance,
                                                 double wLength) const
{
  const double mAtt = EvaluateMieAttenuation(xInit, xFinal, wLength);
  return EvaluateMieScattering(xInit, xFinal, angle, distance, wLength, mAtt);
}


double MeasuredDBMieModel::EvaluateMieScattering(const utl::Point& xInit,
                                                 const utl::Point& /*xFinal*/,
                                                 const double angle,
                                                 const double distance,
                                                 double wLength,
                                                 const double mieAttenuation) const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("Mie model could not find data to compute scattering.");

  // TO BE DONE - CASE OF ZONE CROSSING
  // -- REFACTOR

  const double pf = EvaluateScatteringAngle(xInit, angle, wLength);
  return (1.0 - mieAttenuation) * pf / pow(distance, 2);
}


AttenuationResult
MeasuredDBMieModel::EvaluateMieAttenuation(const utl::Point& xInit,
                                           const utl::Point& xFinal,
                                           const vector<double>&  wLength) const
{
  // Not yet implemented errors !!
  //
  TabulatedFunctionErrors atten;

  for (std::vector<double>::const_iterator wIt = wLength.begin();
       wIt != wLength.end(); ++wIt) {

    const double value = EvaluateMieAttenuation(xInit, xFinal, *wIt);

    // wavelength dependence computed as
    // alpha(wav)=alpha(355nm) * (wav/355nm)^gamma   (gamma positive)
    //
    // This assumes gamma is the same for all heights.  Need to either
    // verify that this is true when reading DB, or else provide a
    // more general method
    //
    atten.PushBack(*wIt, 0., value, 0.);
  }

  return AttenuationResult(atten);
}

double
MeasuredDBMieModel::EvaluateMieAttenuation(const utl::Point& xInit,
                                           const utl::Point& xFinal,
                                           double wLength) const
{
  if  (!CheckForUpdates())
    throw NoDataForModelException("Mie model could not find data to compute attenuation");

  const string zoneName = GetZoneName(xInit);

  // TO BE DONE - CASE OF ZONE CROSSING
  // -- REFACTOR

  // Get height above reference ellipsoid for the specified points
  //
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  double hInit  = (wgs84.PointToLatitudeLongitudeHeight(xInit)).get<2>();
  double hFinal = (wgs84.PointToLatitudeLongitudeHeight(xFinal)).get<2>();

  // THIS CONDITION WILL BE REMOVED WHEN PROPER BOUNDS CHECKING EXISTS
  // IN TABULATEDFUNCTION
  const TabulatedFunction& zoneVAODVsHeight = (*fVAODVsHeightMap)[zoneName];
  const double hMin = zoneVAODVsHeight.XFront();
  const double hMax = zoneVAODVsHeight.XBack();

  const double distance = (xFinal - xInit).GetMag();
  const double deltaHeight = abs(hFinal - hInit);
  const double cosTheta = deltaHeight / distance;
  double averageHeight = min(hInit, hFinal) + 0.5 * deltaHeight;

  //BRD 24/2/06 --> use closest measured VAOD point
  if (hInit < hMin)
    hInit = hMin;
  if (hInit > hMax)
    hInit = hMax;
  if (hFinal < hMin)
    hFinal = hMin;
  if (hFinal > hMax)
    hFinal = hMax;
  if (averageHeight < hMin)
    averageHeight = hMin;
  if (averageHeight > hMax)
    averageHeight = hMax;

  double vaodSlant;

  if (deltaHeight > 1 * m) {
    // Accuracy could probably be improved by using attenuation length
    // for track segments in the top and bottom slices.  Here we are
    // relying on interpolation in VAOD table, for ease of implementation.
    // ALSO, IS 3rd ORDER INTERPOLATION TOO MUCH?
    const double vaodInit = zoneVAODVsHeight.InterpolateY(hInit,3);
    const double vaodFinal = zoneVAODVsHeight.InterpolateY(hFinal,3);
    vaodSlant = abs(vaodFinal - vaodInit) / cosTheta;
  } else {
    const double alpha = (*fAttVsHeightMap)[zoneName].Y(hInit);
    vaodSlant = distance*alpha;
  }

  // Not yet implemented errors !!
  //
  const double gamma = (*fAttLambdaVsHeightMap)[zoneName].Y(averageHeight);

  // This assumes gamma is the same for all heights.  Need to either
  // verify that this is true when reading DB, or else provide a
  // more general method.
  const double arg = -vaodSlant * pow(wLength/355./nanometer, -gamma);

  // Numerical sanity check, though argument should always be <= 0...
  if (arg < std::numeric_limits<double>::max_exponent10 * kLn10)
    return exp(arg);

  return std::numeric_limits<double>::max();
}


double
MeasuredDBMieModel::GetVerticalAerosolOpticalDepth(const unsigned int eyeId, const double altitude)
  const
{
  if  (!CheckForUpdates())
    throw NoDataForModelException("Mie model could not find data to compute VAOD");

  string zoneName = GetZoneName(eyeId);
  const TabulatedFunction& zoneVAODVsHeight = (*fVAODVsHeightMap)[zoneName];

  const double hMax = zoneVAODVsHeight.XBack();
  double VAOD;

  if (altitude >= hMax) {
    ostringstream warn;
    warn << "Requested VAOD at " << altitude / km << " km at eye " << eyeId
         << ", but VAOD profile only goes to " << hMax / km << "km."
         << " Returning heighest measured VAOD.";
    WARNING(warn);

    VAOD = zoneVAODVsHeight.YBack();
  }
  else {
    // IS 3rd ORDER INTERPOLATION TOO MUCH?
    VAOD =  zoneVAODVsHeight.InterpolateY(altitude, 3);
  }

  return VAOD;
}

double
MeasuredDBMieModel::GetAttenuationLength(const utl::Point& p,
                                         const double wLength) const
{
  const string zoneName = GetZoneName(p);
  const TabulatedFunction& zoneAlphaVsHeight = (*fAttVsHeightMap)[zoneName];

  // TO BE DONE - CASE OF ZONE CROSSING
  // -- REFACTOR

  // THIS CONDITION WILL BE REMOVED WHEN PROPER BOUNDS CHECKING EXISTS
  // IN TABULATEDFUNCTION
  const double hMin = zoneAlphaVsHeight.XFront();
  const double hMax = zoneAlphaVsHeight.XBack();

  // Get height above reference ellipsoid for the specified points.
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  double height = (wgs84.PointToLatitudeLongitudeHeight(p)).get<2>();

  // VMH 2018-11-09: when h > hMax we give higher VAODs by repeating the VAOD at hMax,
  // so we must give 0 for alpha (inf for att length) in that case because no change
  // in integrated aerosols means no aerosols at all
  if (height < hMin)
    height = hMin;
  if (height > hMax)
    return numeric_limits<double>::max();

  const double alpha = zoneAlphaVsHeight.Y(height);

  // Not yet implemented errors !!
  const double gamma = (*fAttLambdaVsHeightMap)[zoneName].Y(height);

  // wavelength dependence computed as
  // alpha(wav)=alpha(355nm) * (wav/355nm)^gamma   (gamma positive)
  //
  // This assumes gamma is the same for all heights.  Need to either
  // verify that this is true when reading DB, or else provide a
  // more general method.
  return alpha>0 ? 1. / (alpha * pow(wLength / 355. / nanometer, -gamma)) : numeric_limits<double>::max();
}


double
MeasuredDBMieModel::EvaluateScatteringAngle(const Point& p,
                                            const double angle,
                                            const double /*wLength*/)
  const
{
  const string zoneName = GetZoneName(p);
  const PFParameters& pfParams = (*fPhaseFuncMap)[zoneName];

  const double g2 = pfParams.g * pfParams.g;
  const double mu = cos(angle);
  const double a = 0.25 * (1 - g2) / kPi;
  const double b = pow(1 + g2 - 2 * pfParams.g * mu, -1.5);
  const double c = pfParams.f * 0.5 * (3*mu*mu - 1) * pow(1 + g2, -1.5);

  return a * (b + c);
}


// Returns true if data are found and cached or if data already exist
// for the current detector time.
bool
MeasuredDBMieModel::CheckForUpdates() const
{
  if (fCurrentTime == det::Detector::GetInstance().GetTime() && fCurrentTime)
    return !fDbIsEmptyNow;

  fCurrentTime = det::Detector::GetInstance().GetTime();
  fDbIsEmptyNow = false;

  // Note: these tables are shadow pointers.  No need to call delete.
  fAttVsHeightMap = new ZoneFunction;
  fCentralVAODVsHeightMap = new ZoneFunction;
  fMinVAODVsHeightMap = new ZoneFunction;
  fMaxVAODVsHeightMap = new ZoneFunction;
  fAttLambdaVsHeightMap = new ZoneFunction;
  fPhaseFuncMap = new ZonePhaseFunction;

  fZonePositions.clear();

  const AerosolDB& aerosol1 =
    det::Detector::GetInstance().GetAtmosphere().GetAerosolDB();

  unsigned int zoneCount = 0;
  for (AerosolDB::ZoneIterator zIt = aerosol1.ZonesBegin();
       zIt != aerosol1.ZonesEnd(); ++zIt, ++zoneCount)
  {
    // fill zone name and position
    // height should not be relevant for these calculations
    fZonePositions[zIt->GetName()] =
      UTMPoint(
        zIt->GetNorthing(),
        zIt->GetEasting(),
        0., 19, 'H',
        ReferenceEllipsoid::GetEllipsoidIDFromString("WGS84")).GetPoint();

    TabulatedFunction attVsH;
    TabulatedFunction centralVaodVsH;
    TabulatedFunction minVaodVsH;
    TabulatedFunction maxVaodVsH;
    TabulatedFunction attLambdaVsH;

    // Fill VAOD slices; extrapolate lower end of VAOD table down to VAOD = 0
    // (for more accurate interpolations).
    centralVaodVsH.PushBack(zIt->AttSlicesBegin()->GetMinHeight(), 0.);
    minVaodVsH.PushBack(zIt->AttSlicesBegin()->GetMinHeight(), 0.);
    maxVaodVsH.PushBack(zIt->AttSlicesBegin()->GetMinHeight(), 0.);

    double sliceRoof = 0.;
    double att = 0.;
    double attLambda = 0.;
    for (AerosolZone::AttSliceIterator attIt = zIt->AttSlicesBegin();
         attIt != zIt->AttSlicesEnd(); ++attIt)
    {
      const double hMin = attIt->GetMinHeight();
      const double hMax = attIt->GetMaxHeight();  // DB stores VAOD @ top of each slice

      att       = attIt->GetAttAlpha();
      attLambda = attIt->GetAttLambda();

      try {
        if (fSystematicShift == 0) {
          minVaodVsH.PushBack(hMax, attIt->GetMinVAODUncor());
          maxVaodVsH.PushBack(hMax, attIt->GetMaxVAODUncor());
          centralVaodVsH.PushBack(hMax, attIt->GetVAOD());
        }
        else if (fSystematicShift > 0) {
          const double dVAOD = fSystematicShift*(attIt->GetMaxVAODCor()-attIt->GetVAOD());
          minVaodVsH.PushBack(hMax, attIt->GetMinVAODUncor() + dVAOD);
          maxVaodVsH.PushBack(hMax, attIt->GetMaxVAODUncor() + dVAOD);
          centralVaodVsH.PushBack(hMax, attIt->GetVAOD() + dVAOD);
        }
        else if (fSystematicShift < 0) {
          const double minVaod = attIt->GetMinVAODUncor();
          const double maxVaod = attIt->GetMaxVAODUncor();
          const double vaod = attIt->GetVAOD();
          const double dVAOD =
            fabs(fSystematicShift) * (attIt->GetVAOD()-attIt->GetMinVAODCor());

          if(vaod > dVAOD)
            centralVaodVsH.PushBack(hMax, vaod - dVAOD);
          else
            centralVaodVsH.PushBack(hMax, 0.);

          if(minVaod > dVAOD)
            minVaodVsH.PushBack(hMax, minVaod - dVAOD);
          else
            minVaodVsH.PushBack(hMax, 0.);

          if(maxVaod > dVAOD)
            maxVaodVsH.PushBack(hMax, maxVaod - dVAOD);
          else
            maxVaodVsH.PushBack(hMax, 0.);
        }
      }
      catch (utl::DataNotFoundInDBException& err) {
        ostringstream warn;
        warn << "New Aerosol DB structure (Atm_Aerosol_1_A) not found, trying old DB structure (Atm_Aerosol_0_A)";
        WARNING(warn);

        minVaodVsH.PushBack(hMax, attIt->GetMinVAOD());
        maxVaodVsH.PushBack(hMax, attIt->GetMaxVAOD());
        centralVaodVsH.PushBack(hMax, attIt->GetVAOD());
      }

      attVsH.PushBack(hMin, att);
      attLambdaVsH.PushBack(hMin, attLambda);

      sliceRoof = attIt->GetMaxHeight();
    }

    // Set value at top of highest slice to be the same as at
    // the bottom of same slice (not sure if this is the best thing to do).
    attVsH.PushBack(sliceRoof, att);
    attLambdaVsH.PushBack(sliceRoof, attLambda);

    (*fAttVsHeightMap)[zIt->GetName()] = attVsH;
    (*fCentralVAODVsHeightMap)[zIt->GetName()] = centralVaodVsH;
    (*fMaxVAODVsHeightMap)[zIt->GetName()] = maxVaodVsH;
    (*fMinVAODVsHeightMap)[zIt->GetName()] = minVaodVsH;

    // set VAOD pointer according to user request
    ToggleVAODPtr();

    (*fAttLambdaVsHeightMap)[zIt->GetName()] = attLambdaVsH;

    // Fill PF slices
    for (AerosolZone::PFSliceIterator pfIt = zIt->PFSlicesBegin();
         pfIt != zIt->PFSlicesEnd() ; ++pfIt) {

      PFParameters pfParams;
      pfParams.f          = pfIt->GetF();
      pfParams.fErr       = pfIt->GetFError();
      pfParams.g          = pfIt->GetG();
      pfParams.gErr       = pfIt->GetGError();
      pfParams.fLambda    = pfIt->GetFLambda();
      pfParams.fLambdaErr = pfIt->GetFLambdaError();
      pfParams.gLambda    = pfIt->GetGLambda();
      pfParams.gLambdaErr = pfIt->GetGLambdaError();

      (*fPhaseFuncMap)[zIt->GetName()] = pfParams;
    }
  }
  if (zoneCount == 0) {
    fDbIsEmptyNow = true;

     ostringstream msg;
     msg << "Mie model requested data from database, but no data found for time ";
     msg << det::Detector::GetInstance().GetTime();
     DEBUGLOG(msg);

    return false;
  }
  return true;
}


// First attempt at zone selection.  Pick the zone with coordinates closest to
// the given point.
string
MeasuredDBMieModel::GetZoneName(const utl::Point& pt) const
{
  string zoneName("");

  if (fZoneParam == eAll) {
    double smallestDist = numeric_limits<float>::max(); // was HUGE
    for (map<string, Point>::iterator zIt = fZonePositions.begin();
         zIt != fZonePositions.end(); ++zIt) {

      const double dist = (pt - zIt->second).GetMag();

      if (dist <= smallestDist) {
        smallestDist = dist;
        zoneName = zIt->first;
      }
    }

    // Zone swapping: see if zone is in the swap list, and swap the zone name
    if (fSwapZones.size()) {
      string zone1 = ZoneParamToString(fSwapZones[0]);
      string zone2 = ZoneParamToString(fSwapZones[1]);
      bool isZone1 = false;
      bool isZone2 = false;

      if ((isZone1 = ContainsSubstringEquivalent(zoneName, zone1)) ||
          (isZone2 = ContainsSubstringEquivalent(zoneName, zone2)))
      {
        ZonePosition::iterator zPosIt;
        if (isZone1 && HasZone(zone2))
          zPosIt = find_if(fZonePositions.begin(), fZonePositions.end(),
                           [&](const auto& zone) { return ZoneMatch(zone, zone2); });
        else if (isZone2 && HasZone(zone1))
          zPosIt = find_if(fZonePositions.begin(), fZonePositions.end(),
                           [&](const auto& zone) { return ZoneMatch(zone, zone1); });
        else {
          ostringstream error;
          error << "Zone swap: missing zone for \""
                << (isZone1 ? zone2 : (isZone2 ? zone1 : "UNKNOWN")) << "\"";
          throw DataNotFoundInDBException(error.str());
        }
        zoneName = zPosIt->first;
      }
    }
    return zoneName;
   }
  // Zone forcing enabled
  else {
    ZonePosition::iterator zPosIt =
      find_if(fZonePositions.begin(), fZonePositions.end(),
              [&](const auto& zone) { return ZoneMatch(zone, ZoneParamToString(fZoneParam)); });

    if (zPosIt != fZonePositions.end())
      return zPosIt->first;
    else {
      ostringstream error;
      error << "Zone force: missing zone for \""
            << ZoneParamToString(fZoneParam) << "\"";
      throw DataNotFoundInDBException(error.str());
    }
  }
}

// Pick the aerosol zone nearest to the given eye.
string
MeasuredDBMieModel::GetZoneName(const unsigned int eyeId) const
{
  const utl::Point& eyePos =
    det::Detector::GetInstance().GetFDetector().GetEye(eyeId).GetPosition();

  return GetZoneName(eyePos);
}


bool
MeasuredDBMieModel::CrossesZone(const utl::Point& /*pt1*/, const utl::Point& /*p2*/) const
{
  // To be implemented. Return true if the two points are not in the same zone.
  return false;
}


void
MeasuredDBMieModel::SetUncertaintyBound(double nSigma) const
{
  fNSigma = nSigma;
  ToggleVAODPtr();
}


void
MeasuredDBMieModel::ToggleVAODPtr() const
{
  if (fNSigma == 0.)
    fVAODVsHeightMap = fCentralVAODVsHeightMap.Get();
  else if (fNSigma == 1.)
    fVAODVsHeightMap = fMaxVAODVsHeightMap.Get();
  else if (fNSigma == -1.)
    fVAODVsHeightMap = fMinVAODVsHeightMap.Get();
  else {
    ostringstream err;
    err << " Mie model could not find data for " << fNSigma
        << " sigma uncertainty bound";
    throw NoDataForModelException(err.str());
  }
}


MeasuredDBMieModel::ZoneParam
MeasuredDBMieModel::StringToZoneParam(const std::string& zoneString) const
{
  if (ContainsSubstringEquivalent(zoneString, "Leones"))
    return eLosLeones;
  else if (ContainsSubstringEquivalent(zoneString, "Morados"))
    return eLosMorados;
  else if (ContainsSubstringEquivalent(zoneString, "Amarilla"))
    return eLomaAmarilla;
  else if (ContainsSubstringEquivalent(zoneString, "Coihueco"))
    return eCoihueco;
  else
    return eAll;
}


std::string
MeasuredDBMieModel::ZoneParamToString(const ZoneParam& zone) const
{
  string zoneName("");

  switch (zone) {
    case eLosLeones:
      zoneName = "Leones";
      break;
    case eLosMorados:
      zoneName = "Morados";
      break;
    case eLomaAmarilla:
      zoneName = "Amarilla";
      break;
    case eCoihueco:
      zoneName = "Coihueco";
      break;
    default:
      break;
  }

  return zoneName;
}


bool
MeasuredDBMieModel::ZoneMatch(const ZonePoint& zonePt, const string& zoneName) const
{
  return ContainsSubstringEquivalent(zonePt.first, zoneName);
}