Quadrature.h

Go to the documentation of this file.

#ifndef QUADRATURE_H_
#define QUADRATURE_H_

#ifdef __cplusplus
#ifdef SWIG
#include "Python.h"
#endif
#include "constants.h"
#include "log.h"
#include <sstream>
#endif
#include <vector>
#include <iostream>
#include <cmath>
#include <algorithm>


enum QuadratureType {
  TABUCHI_YAMAMOTO,
  LEONARD,
  GAUSS_LEGENDRE,
  EQUAL_WEIGHT,
  EQUAL_ANGLE
};

typedef std::vector<double> DoubleVec;
typedef DoubleVec::const_iterator DVCI;

template <typename T>
  std::ostream& operator<<(std::ostream& os, const std::vector<T>& vec) {
  if (vec.size() == 0) {
    os << "EMPTY";
  }

  else {

    using CI = typename std::vector<T>::const_iterator;

    os << "[";
    for (CI cit = vec.begin(); cit != vec.end(); ++cit) {
      if (cit != vec.begin()) {
        os << ", ";
      }
      os << *cit;
    }
    os << "]";
  }
  return os;
}



class Quadrature {

protected:

  QuadratureType _quad_type;

  size_t _num_azim;

  size_t _num_polar;

  std::vector<DoubleVec> _sin_thetas;

  std::vector<DoubleVec> _thetas;

  DoubleVec _phis;

  DoubleVec _azim_spacings;

  DoubleVec _azim_weights;

  std::vector<DoubleVec> _polar_spacings;

  std::vector<DoubleVec> _polar_weights;

  std::vector<DoubleVec> _total_weights;

  /* Templates for setting the same values to complimentary and supplementary
   * angles */
  template <typename T>
  void setPolarValues(std::vector< std::vector<T> >& vec, size_t azim_index,
                      size_t polar_index, T value) {
    size_t c_azim_index  = _num_azim/2 - azim_index - 1;
    size_t s_polar_index = _num_polar - polar_index - 1;
    vec.at(azim_index).at(polar_index)     = value;
    vec.at(c_azim_index).at(polar_index)   = value;
    vec.at(azim_index).at(s_polar_index)   = value;
    vec.at(c_azim_index).at(s_polar_index) = value;
  }

  template <typename T>
  void setAzimuthalValues(std::vector<T>& vec, size_t azim_index, T value) {
    vec.at(azim_index)                   = value;
    vec.at(_num_azim/2 - azim_index - 1) = value;
  }

  template <typename T>
  static void resize2D(std::vector< std::vector<T> >& vec, size_t dim1,
                       size_t dim2) {
    vec.resize(dim1);
    for (size_t i = 0; i < dim1; ++i)
    {
        vec.at(i).resize(dim2);
    }
  }

public:

  Quadrature();

  size_t getNumPolarAngles() const;
  size_t getNumAzimAngles() const;
  double getSinTheta(size_t azim, size_t polar) const;
  double getSinThetaInline(size_t azim, size_t polar) const;
  double getTheta(size_t azim, size_t polar) const;
  double getPhi(size_t azim) const;
  double getAzimWeight(size_t azim) const;
  double getPolarWeight(size_t azim, size_t polar) const;
  double getWeight(size_t azim, size_t polar) const;
  double getWeightInline(size_t azim, size_t polar) const;
  const std::vector<DoubleVec>& getSinThetas() const;
  const std::vector<DoubleVec>& getThetas() const;
  const DoubleVec& getPhis() const;
  const DoubleVec& getAzimWeights() const;
  const std::vector<DoubleVec>& getPolarWeights() const;
  QuadratureType getQuadratureType() const;
  const DoubleVec& getAzimSpacings() const;
  double getAzimSpacing(size_t azim) const;
  const std::vector<DoubleVec>& getPolarSpacings() const;
  double getPolarSpacing(size_t azim, size_t polar) const;

  virtual void setNumAzimAngles(size_t num_azim);
  virtual void setNumPolarAngles(size_t num_polar);
  void setThetas(const DoubleVec& thetas);
  void setPolarWeights(const DoubleVec& weights);
  void setTheta(double theta, size_t azim, size_t polar);
  void setPhi(double phi, size_t azim);
  void setAzimSpacing(double spacing, size_t azim);
  void setPolarSpacing(double spacing, size_t azim, size_t polar);
  void setAzimWeight(double weight, size_t azim);
  void setPolarWeight(double weight, size_t azim, size_t polar);

  virtual void initialize();
  virtual void precomputeWeights(bool solve_3D);

  std::string toString() const;

  friend std::ostream& operator<<(std::ostream& os, const Quadrature& quad);
};



class TYPolarQuad: public Quadrature {

private:

public:
  TYPolarQuad();
  void setNumPolarAngles(size_t num_polar);
  void initialize();
  void precomputeWeights(bool solve_3D);
};


class LeonardPolarQuad: public Quadrature {

private:

public:
  LeonardPolarQuad();
  void setNumPolarAngles(size_t num_polar);
  void initialize();
  void precomputeWeights(bool solve_3D);
};



class GLPolarQuad: public Quadrature {

private:

  std::vector <double> _roots;

  std::vector <double> _adjusted_roots;

  bool _correct_weights;

public:
  GLPolarQuad();
  void setNumPolarAngles(size_t num_polar);
  void useCorrectedWeights(bool use_corrected_weights);
  void initialize();
  void precomputeWeights(bool solve_3D);
  DoubleVec getCorrectedWeights(size_t azim) const;

  static double legendrePolynomial(size_t n, double x);
  static double logDerivLegendre(size_t n, double x);
  static double secondLogDerivLegendre(size_t n, double x);
  static double getSingleWeight(double root, size_t n);
  static DoubleVec getLegendreRoots(size_t n);
  static DoubleVec getGLWeights(const DoubleVec& roots, size_t n);
};



class EqualWeightPolarQuad: public Quadrature {

private:

public:
  EqualWeightPolarQuad();
  void setNumPolarAngles(size_t num_polar);
  void initialize();
  void precomputeWeights(bool solve_3D);
};



class EqualAnglePolarQuad: public Quadrature {

private:

public:
  EqualAnglePolarQuad();
  void setNumPolarAngles(size_t num_polar);
  void initialize();
  void precomputeWeights(bool solve_3D);
};


inline double Quadrature::getSinThetaInline(size_t azim, size_t polar) const {

  if (azim >= _num_azim/2)
    azim = _num_azim - azim - 1;

  return _sin_thetas[azim][polar];
}


inline double Quadrature::getWeightInline(size_t azim, size_t polar) const {
  return _total_weights[azim][polar];
}

#endif /* QUADRATURE_H_ */