absco.h

Go to the documentation of this file.

#ifndef ABSCO_H
#define ABSCO_H
#include "gas_absorption.h"
#include "array_ad_with_unit.h"
#include "array_with_unit.h"
#include "null_deleter.h"
#include <boost/any.hpp>
#include <vector>

namespace FullPhysics {
class Absco;

/****************************************************************/
class AbscoInterpolator : public Printable<AbscoInterpolator> {
public:
  AbscoInterpolator(const boost::shared_ptr<Absco>& A,
                    const ArrayWithUnit<double, 1>& Press, 
                    const ArrayAdWithUnit<double, 1>& Temp,
                    const ArrayAdWithUnit<double, 1>& Broadener_vmr);
  virtual ~AbscoInterpolator() {}
  virtual void print(std::ostream& Os) const {Os << "AbscoInterpolator";}

  blitz::Array<double, 1> absorption_cross_section_noderiv(double wn) const;
  ArrayAd<double, 1> absorption_cross_section_deriv(double wn) const;
private:
  double interpol(double X, const std::vector<double>& Xv, 
                  int& i, double& df_dx) const;
  template<class T> blitz::Array<double, 1> 
  absorption_cross_section_noderiv_calc(double wn) const;
  template<class T> ArrayAd<double, 1> 
  absorption_cross_section_deriv_calc(double wn) const;
  boost::shared_ptr<Absco> absco;
  blitz::Array<double, 1> p;     // This is in Pascals
  ArrayAd<double, 1> t;         // This is in K
  ArrayAd<double, 1> b;         // This is dimensionless
  blitz::Array<double, 2> t_jac; // t.jacobian(), in C order and contiguous
  blitz::Array<double, 2> b_jac; // b.jacobian(), in C order and contiguous
  // Various indexes used in the interpolation
  blitz::Array<int, 1> ip, itp1, itp2, ib, ib2;
  blitz::Array<double, 1> dftp1_dt, dftp2_dt, dfb_db,
    fp, fb, ftp1, ftp2;
  mutable ArrayAd<double, 1> res; // We return the same results array
                                // each time, so we don't need to keep
                                // creating and destroying this.
};

/****************************************************************/
class Absco: public GasAbsorption {
public:
  virtual ~Absco() {}

  virtual std::string broadener_name() const = 0;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual double table_scale(double wn) const = 0;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual int number_broadener_vmr() const 
  { return broadener_vmr_grid().rows(); }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual int number_layer() const 
  {return pressure_grid().rows(); }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual int number_temperature() const 
  {return temperature_grid().cols(); }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 1> broadener_vmr_grid() const = 0;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 1> pressure_grid() const = 0;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 2> temperature_grid() const = 0;


  virtual DoubleWithUnit absorption_cross_section
  (double Wn, const DoubleWithUnit& Press, const DoubleWithUnit& Temp,
   const DoubleWithUnit& Broadener_vmr) const;
  virtual AutoDerivativeWithUnit<double>
  absorption_cross_section(double wn, 
    const DoubleWithUnit& press, 
    const AutoDerivativeWithUnit<double>& temp,
    const AutoDerivativeWithUnit<double>& broadener_vmr) const;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual bool is_float() const = 0;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------
  template <class T> blitz::Array<T, 3> read(double wn) const;
  friend class AbscoInterpolator;
protected:

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  virtual blitz::Array<double, 3> read_double(double wn) const = 0;
  virtual blitz::Array<float, 3> read_float(double wn) const = 0;
private:
  double interpol(double X, const std::vector<double>& Xv, 
                  int& i, double& df_dx) const;
  mutable std::vector<double> pgrid;
  mutable std::vector<std::vector<double> > tgrid;
  mutable std::vector<double> bgrid;
  void fill_pgrid_tgrid_and_bgrid() const;
};

template <> inline blitz::Array<double, 3> Absco::read<double>(double wn) const
{
  return read_double(wn);
}

template <> inline blitz::Array<float, 3> Absco::read<float>(double wn) const
{
  return read_float(wn);
}

}

#endif