ParametricXMLProfileModel.cc

Go to the documentation of this file.

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TabulatedFunction.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/CentralConfig.h>

#include <atm/ProfileResult.h>
#include <atm/ParametricXMLProfileModel.h>

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

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

ParametricXMLProfileModel::ParametricXMLProfileModel() : VProfileModel() {
}

ParametricXMLProfileModel::~ParametricXMLProfileModel() {
}

void ParametricXMLProfileModel::Init(){
  
  CentralConfig* cc = CentralConfig::GetInstance();

  Branch topB =
    cc->GetTopBranch("ParametricXMLProfileModel");

  vector<double> pressure;
  vector<double> height;
  vector<double> density;
  vector<double> temperature;
  vector<double> humidity;

  topB.GetChild("height").GetData(height);
  topB.GetChild("pressure").GetData(pressure);
  topB.GetChild("density").GetData(density);
  topB.GetChild("temperature").GetData(temperature);
  topB.GetChild("humidity").GetData(humidity);

  if ((height.size()!=pressure.size())||
      (height.size()!=density.size())||
      (height.size()!=temperature.size())||
      (height.size()!=humidity.size()))
  {
    ERROR("Dimensions of the atmosphere profile vectors are different. ");
    throw utl::NonExistentComponentException();
  }

  // make tabulated functions
  TabulatedFunctionErrors tabRIVsHeight;
  TabulatedFunctionErrors tabLogPressureVsHeight;
  TabulatedFunctionErrors tabLogVaporPressureVsHeight;
  TabulatedFunctionErrors tabTemperatureVsHeight;
  TabulatedFunctionErrors tabLogDepthVsHeight;
  TabulatedFunctionErrors tabLogDensityVsHeight;
  TabulatedFunctionErrors tabHeightVsLogDepth;

  const double minLog = kLn10 * numeric_limits<double>::min_exponent10;

  for (unsigned int i = 0 ; i < height.size() ; i++) {

    tabLogDepthVsHeight.PushBack(
          height[i], 0.,
          log(pressure[i]*(kOverburdenSeaLevel/atmosphere)), minLog);

    // Move calculation of refraction index to base class once
    // TabulatedFunction is fixed
    const double refracIndex = utl::RefractionIndex::LorentzLorentz(pressure[i]*(kOverburdenSeaLevel/atmosphere));
    tabRIVsHeight.PushBack(height[i], 0., refracIndex, 0.);

    tabLogPressureVsHeight.PushBack(height[i], 0., log(pressure[i]), minLog);

    tabLogVaporPressureVsHeight.PushBack(height[i], 0.,
        (humidity[i] > 0.) ? log(humidity[i]) : minLog, minLog);

    tabTemperatureVsHeight.PushBack(height[i], 0., temperature[i], 0.);

    tabLogDensityVsHeight.PushBack(height[i], 0., log(density[i]), minLog);
  }

  // have to fill tabulated functions in ascending order or ordinate (tabulated function should be fixed..)
  for (int i = height.size()-1 ; i >= 0 ; i--) {
    tabHeightVsLogDepth.PushBack(
      log(pressure[i]*(kOverburdenSeaLevel/atmosphere)), minLog, height[i], 0.);
  }

  INFO("CREATING PROFILE RESULTS");

  // create profile results
  *fTabLogPressureVsHeight = ProfileResult(
      tabLogPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogVaporPressureVsHeight = ProfileResult(
      tabLogVaporPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabTemperatureVsHeight = ProfileResult(
      tabTemperatureVsHeight, ProfileResult::eLinear, ProfileResult::eLinear);

  *fTabLogDensityVsHeight = ProfileResult(
      tabLogDensityVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogDepthVsHeight = ProfileResult(
      tabLogDepthVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabHeightVsLogDepth = ProfileResult(
      tabHeightVsLogDepth, ProfileResult::eLog, ProfileResult::eLinear);

  *fTabRIVsHeight = ProfileResult(
      tabRIVsHeight, ProfileResult::eLinear, ProfileResult::eLinear);

  INFO("DONE CREATING PROFILE RESULTS");

}

const ProfileResult&
ParametricXMLProfileModel::EvaluatePressureVsHeight() const
{
  return *fTabLogPressureVsHeight;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateVaporPressureVsHeight() const
{
  return *fTabLogVaporPressureVsHeight;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateTemperatureVsHeight() const
{
  return *fTabTemperatureVsHeight;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateDensityVsHeight() const
{
  return *fTabLogDensityVsHeight;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateDepthVsHeight() const
{
  return *fTabLogDepthVsHeight;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateHeightVsDepth() const
{
  return *fTabHeightVsLogDepth;
}

const ProfileResult&
ParametricXMLProfileModel::EvaluateRefractionIndexVsHeight() const
{
  return *fTabRIVsHeight;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: