diff --git a/src/cdi_analytic/Analytic_Models.cc b/src/cdi_analytic/Analytic_Models.cc index 392781fe36..11fc71730b 100644 --- a/src/cdi_analytic/Analytic_Models.cc +++ b/src/cdi_analytic/Analytic_Models.cc @@ -185,9 +185,10 @@ Constant_Analytic_Opacity_Model::get_parameters() const { Polynomial_Analytic_Opacity_Model::Polynomial_Analytic_Opacity_Model( const sf_char &packed) - : a(0.0), b(0.0), c(0.0), d(0.0), e(0.0), f(1.0), g(1.0), h(1.0) { + : a(0.0), b(0.0), c(0.0), d(0.0), e(0.0), f(1.0), g(1.0), h(1.0), i(0.0), + j(0.0), k(0.0) { // size of stream - size_t size = sizeof(int) + 8 * sizeof(double); + size_t size = sizeof(int) + 11 * sizeof(double); Require(packed.size() == size); @@ -204,7 +205,7 @@ Polynomial_Analytic_Opacity_Model::Polynomial_Analytic_Opacity_Model( "Tried to unpack the wrong type in Polynomial_Analytic_Opacity_Model"); // unpack the data - unpacker >> a >> b >> c >> d >> e >> f >> g >> h; + unpacker >> a >> b >> c >> d >> e >> f >> g >> h >> i >> j >> k; Ensure(unpacker.get_ptr() == unpacker.end()); } @@ -217,8 +218,8 @@ Polynomial_Analytic_Opacity_Model::pack() const { // get the registered indicator int indicator = POLYNOMIAL_ANALYTIC_OPACITY_MODEL; - // caculate the size in bytes: indicator + 8 * double - int size = sizeof(int) + 8 * sizeof(double); + // caculate the size in bytes: indicator + 11 * double + int size = sizeof(int) + 11 * sizeof(double); // make a vector of the appropriate size sf_char pdata(size); @@ -241,6 +242,9 @@ Polynomial_Analytic_Opacity_Model::pack() const { packer << f; packer << g; packer << h; + packer << i; + packer << j; + packer << k; // Check the size Ensure(packer.get_ptr() == &pdata[0] + size); @@ -253,7 +257,7 @@ Polynomial_Analytic_Opacity_Model::pack() const { Analytic_Opacity_Model::sf_double Polynomial_Analytic_Opacity_Model::get_parameters() const { - sf_double p(8); + sf_double p(11); p[0] = a; p[1] = b; p[2] = c; @@ -262,105 +266,12 @@ Polynomial_Analytic_Opacity_Model::get_parameters() const { p[5] = f; p[6] = g; p[7] = h; + p[8] = i; + p[9] = j; + p[10] = k; return p; } -//===========================================================================// -// STIMULATED_EMISSION_ANALYTIC_OPACITY_MODEL DEFINITIONS -//===========================================================================// -/*! - * \brief Unpacking constructor for Stimulated_Emission_Analytic_Opacity_Model - * objects. - * \bug There are no unit tests for this function. - */ -Stimulated_Emission_Analytic_Opacity_Model:: - Stimulated_Emission_Analytic_Opacity_Model(const sf_char &packed) - : a(0.0), b(0.0), c(0.0), d(0.0), e(0.0), f(1.0), g(1.0), h(1.0) { - // size of stream - size_t size = sizeof(int) + 8 * sizeof(double); - - Require(packed.size() == size); - - // make an unpacker - rtt_dsxx::Unpacker unpacker; - - // set the unpacker - unpacker.set_buffer(size, &packed[0]); - - // unpack the indicator - int indicator; - unpacker >> indicator; - Insist(indicator == STIMULATED_EMISSION_ANALYTIC_OPACITY_MODEL, - "Tried to unpack the wrong type in " - "Stimulated_Emission_Analytic_Opacity_Model"); - - // unpack the data - unpacker >> a >> b >> c >> d >> e >> f >> g >> h; - - Ensure(unpacker.get_ptr() == unpacker.end()); -} - -//---------------------------------------------------------------------------// -/*! - * \brief Packing functionfor Stimulated_Emission_Analytic_Opacity_Model - * objects. - * \bug There are no unit tests for this function, so commenting it out. - */ -Analytic_Opacity_Model::sf_char -Stimulated_Emission_Analytic_Opacity_Model::pack() const { - // get the registered indicator - int indicator = STIMULATED_EMISSION_ANALYTIC_OPACITY_MODEL; - - // caculate the size in bytes: indicator + 8 * double - int size = sizeof(int) + 8 * sizeof(double); - - // make a vector of the appropriate size - sf_char pdata(size); - - // make a packer - rtt_dsxx::Packer packer; - - // set the packer buffer - packer.set_buffer(size, &pdata[0]); - - // pack the indicator - packer << indicator; - - // pack the data - packer << a; - packer << b; - packer << c; - packer << d; - packer << e; - packer << f; - packer << g; - packer << h; - - // Check the size - Ensure(packer.get_ptr() == &pdata[0] + size); - - return pdata; -} - -//---------------------------------------------------------------------------// -/*! - * \brief Return the model parameters for a - * Stimulated_Emission_Analytic_Opacity_Model - * \bug There are no unit tests for this, so commenting it out. - */ -Analytic_Opacity_Model::sf_double -Stimulated_Emission_Analytic_Opacity_Model::get_parameters() const { - sf_double p(8); - p[0] = a; - p[1] = b; - p[2] = c; - p[3] = d; - p[4] = e; - p[5] = f; - p[6] = g; - p[7] = h; - return p; -} //===========================================================================// // POLYNOMIAL_SPECIFIC_HEAT_ANALYTIC_EOS_MODEL DEFINITIONS diff --git a/src/cdi_analytic/Analytic_Models.hh b/src/cdi_analytic/Analytic_Models.hh index c3dbc9dd02..e41a68ffff 100644 --- a/src/cdi_analytic/Analytic_Models.hh +++ b/src/cdi_analytic/Analytic_Models.hh @@ -32,7 +32,6 @@ namespace rtt_cdi_analytic { enum Opacity_Models { CONSTANT_ANALYTIC_OPACITY_MODEL, POLYNOMIAL_ANALYTIC_OPACITY_MODEL, - STIMULATED_EMISSION_ANALYTIC_OPACITY_MODEL, }; //----------------------------------------------------------------------------// @@ -158,222 +157,109 @@ public: * * The opacity is defined: * - * \arg opacity = (a + b * T^c * nu^e) * rho^d + * \arg opacity = a + (T/f)^c * (rho/g)^d * (nu/h)^e + * (1 - i * exp(-nu/T)) * (b + j * H(nu - k)) * - * where the coefficients have the following units: + * where i <= 0 means no stimulated emission correction + * and H is the Heaviside function. + * + * The coefficients are unitless or have the following units: * - * \arg a = [cm^2/g * (cm^3/g)^d] - * \arg b = [keV^(-c) * cm^2/g * (cm^3/g)^d] + * \arg a = [cm^2/g] + * \arg b = [cm^2/g] + * \arg f = [keV] + * \arg g = [g/cm^3] + * \arg h = [keV] + * \arg j = [cm^2/g] + * \arg k = [keV] */ class Polynomial_Analytic_Opacity_Model : public Analytic_Opacity_Model { private: // Coefficients - double a; //!< constant [cm^2/g * (cm^3/g)^d] - double b; //!< temperature multiplier [keV^(-c) * cm^2/g * (cm^3/g)^d] + double a; //!< constant [cm^2/g] + double b; //!< temperature multiplier [cm^2/g] double c; //!< temperature power double d; //!< density power double e; //!< frequency power - double f; //!< reference temperature - double g; //!< reference density - double h; //!< reference frequency + double f; //!< reference temperature [keV] + double g; //!< reference density [g/cm^3] + double h; //!< reference frequency [keV] + double i; //!< stimulated emission [0 or 1] + double j; //!< edge strength [cm^2/g] + double k; //!< edge location [keV] public: /*! * \brief Constructor. - * \param[in] a_ constant [cm^2/g (cm^3/g)^d] - * \param[in] b_ temperature multiplier [keV^(-c) cm^2/g (cm^3/g)^d] + * \param[in] a_ constant [cm^2/g] + * \param[in] b_ temperature multiplier [cm^2/g] * \param[in] c_ temperature power * \param[in] d_ density power * \param[in] e_ frequency power (default = 0) - * \param[in] f_ reference temperature (default = 1) - * \param[in] g_ reference density (default = 1) - * \param[in] h_ reference frequency (default = 1) + * \param[in] f_ reference temperature (default = 1 [keV]) + * \param[in] g_ reference density (default = 1 [g/cm^3]) + * \param[in] h_ reference frequency (default = 1 [keV]) + * \param[in] i_ stimulated emission (default = 0 [off]) + * \param[in] j_ edge strength (default = 0 [cm^2/g]) + * \param[in] k_ edge location (default = 0 [keV]) */ Polynomial_Analytic_Opacity_Model(double a_, double b_, double c_, double d_, double e_ = 0, double f_ = 1, double g_ = 1, - double h_ = 1) - : a(a_), b(b_), c(c_), d(d_), e(e_), f(f_), g(g_), h(h_) { + double h_ = 1, double i_ = 0, double j_ = 0, + double k_ = 0) + : a(a_), b(b_), c(c_), d(d_), e(e_), f(f_), g(g_), h(h_), i(i_), j(j_), + k(k_) { /*...*/ } //! Constructor for packed state. explicit Polynomial_Analytic_Opacity_Model(const sf_char &packed); - //! Calculate the opacity in units of cm^2/g - double calculate_opacity(double T, double rho, double nu0, double nu1) const { - using std::pow; - Require(c < 0.0 ? T > 0.0 : T >= 0.0); - Require(rho >= 0.0); - Require(nu1 > nu0); - Require(f > 0.0); - Require(g > 0.0); - Require(h > 0.0); - - double opacity(-1.0); - - //double nu = 0.5*(nu0+nu1); - double nu = sqrt(nu0 * nu1); - opacity = (a + b * pow(T / f, c) * pow(nu / h, e)) * pow(rho / g, d); - - Ensure(opacity >= 0.0); - return opacity; - } - //! Calculate the opacity in units of cm^2/g double calculate_opacity(double T, double rho, double nu) const { using std::pow; Require(c < 0.0 ? T > 0.0 : T >= 0.0); + Require(i > 0.0 ? T > 0.0 : T >= 0.0); Require(rho >= 0.0); Require(nu >= 0.0); Require(f > 0.0); Require(g > 0.0); Require(h > 0.0); - double opacity = (a + b * pow(T / f, c) * pow(nu / h, e)) * pow(rho / g, d); - - Ensure(opacity >= 0.0); - return opacity; - } - - //! Calculate the opacity in units of cm^2/g - double calculate_opacity(double T, double rho) const { - using std::pow; - Require(c < 0.0 ? T > 0.0 : T >= 0.0); - Require(rho >= 0.0); - - double opacity = (a + b * pow(T / f, c)) * pow(rho / g, d); + const double pows = pow(T / f, c) * pow(nu / h, e) * pow(rho / g, d); + const double stim = (i <= 0.0) ? 1.0 : 1.0 - exp(-nu / T); + const double jH = (nu >= k) ? j : 0.0; + const double opacity = a + pows * stim * (b + jH); Ensure(opacity >= 0.0); return opacity; } - //! Return the model parameters. - sf_double get_parameters() const; - - //! Pack up the class for persistence. - sf_char pack() const; -}; - -//---------------------------------------------------------------------------// -/*! - * \class Stimulated_Emission_Analytic_Opacity_Model - * \brief Derived Analytic_Opacity_Model class that defines a polynomial - * function for the opacity that includes stimulated emission. - * - * The opacity is defined: - * - * \arg opacity = (a + b * T^c * nu^e [ 1 - exp(-nu/T) ] ) * rho^d, - * - * i.e., - * - * \f[ \sigma = \left[a+bT^c \nu^e \left(1-e^{-\nu/T}\right) \right] \rho^d \f] - * - * where the coefficients have the following units: - * - * \arg a = [cm^2/g * (cm^3/g)^d] - * \arg b = [keV^(-c) * cm^2/g * (cm^3/g)^d] - * - * For a typical bound-free or free-free absorption model, one should set - * \arg a = 0 - * \arg b = constant > 0 - * \arg c = -1/2 ( inverse sqrt(T) ) - * \arg d = constant - * \arg e = -3 - * - * This produces a Planck or Rosseland opacity of the form \f[ \overline{\sigma} - * \propto T^{-7/2} \f] which is also known as Kramers' Opacity Law. - * - * \sa{ http://en.wikipedia.org/wiki/Kramers'_opacity_law } - */ -class Stimulated_Emission_Analytic_Opacity_Model - : public Analytic_Opacity_Model { -private: - // Coefficients - double a; //!< constant [cm^2/g * (cm^3/g)^d] - double b; //!< temperature multiplier [keV^(-c) * cm^2/g * (cm^3/g)^d] - double c; //!< temperature power - double d; //!< density power - double e; //!< frequency power - double f; //!< reference temperature - double g; //!< reference density - double h; //!< reference frequency - -public: - /*! - * \brief Constructor. - * \param a_ constant [cm^2/g (cm^3/g)^d] - * \param b_ temperature multiplier [keV^(-c) cm^2/g (cm^3/g)^d] - * \param c_ temperature power - * \param d_ density power - * \param e_ frequency power - * \param f_ reference temperature - * \param g_ reference density - * \param h_ reference frequency - * - * \bug no unit test exists. - */ - Stimulated_Emission_Analytic_Opacity_Model(double a_, double b_, double c_, - double d_, double e_ = 0, - double f_ = 1, double g_ = 1, - double h_ = 1) - : a(a_), b(b_), c(c_), d(d_), e(e_), f(f_), g(g_), h(h_) { - /*...*/ - } - - //! Constructor for packed state. - //! \bug no unit test exists. - explicit Stimulated_Emission_Analytic_Opacity_Model(const sf_char &packed); - //! Calculate the opacity in units of cm^2/g double calculate_opacity(double T, double rho, double nu0, double nu1) const { - using std::pow; - Require(T > 0.0); - Require(rho >= 0.0); Require(nu1 > nu0); - Require(f > 0.0); - Require(g > 0.0); - Require(h > 0.0); //double nu = 0.5*(nu0+nu1); double nu = sqrt(nu0 * nu1); - double opacity = (a + b * pow(T / f, c) * pow(nu / h, e) * - (1 - exp(-(nu / h) / (T / f)))) * - pow(rho / g, d); - - Ensure(opacity >= 0.0); - return opacity; + return calculate_opacity(T, rho, nu); } //! Calculate the opacity in units of cm^2/g - double calculate_opacity(double T, double rho, double nu) const { + double calculate_opacity(double T, double rho) const { using std::pow; - Require(T > 0.0); + Require(c < 0.0 ? T > 0.0 : T >= 0.0); Require(rho >= 0.0); - Require(nu >= 0.0); - Require(f > 0.0); - Require(g > 0.0); - Require(h > 0.0); - double opacity = (a + b * pow(T / f, c) * pow(nu / h, e) * - (1 - exp(-(nu / h) / (T / f)))) * - pow(rho / g, d); + const double opacity = a + b * pow(T / f, c) * pow(rho / g, d); Ensure(opacity >= 0.0); return opacity; } - //! Calculate the opacity in units of cm^2/g - double calculate_opacity(double, double) const { - Insist(false, "Stimatulated emission opacity model needs a frequency."); - return -1.0; - } - //! Return the model parameters. - //! \bug no unit test exists. sf_double get_parameters() const; //! Pack up the class for persistence. - //! \bug no unit test exists. sf_char pack() const; }; diff --git a/src/cdi_analytic/nGray_Analytic_Odfmg_Opacity.cc b/src/cdi_analytic/nGray_Analytic_Odfmg_Opacity.cc index ee6c3e63aa..233ffa036b 100644 --- a/src/cdi_analytic/nGray_Analytic_Odfmg_Opacity.cc +++ b/src/cdi_analytic/nGray_Analytic_Odfmg_Opacity.cc @@ -114,10 +114,6 @@ nGray_Analytic_Odfmg_Opacity::nGray_Analytic_Odfmg_Opacity( } else if (indicator == rtt_cdi_analytic::POLYNOMIAL_ANALYTIC_OPACITY_MODEL) { group_models[i].reset(new Polynomial_Analytic_Opacity_Model(models[i])); - } else if (indicator == - rtt_cdi_analytic::STIMULATED_EMISSION_ANALYTIC_OPACITY_MODEL) { - group_models[i].reset( - new Stimulated_Emission_Analytic_Opacity_Model(models[i])); } else { Insist(false, "Unregistered analytic opacity model!"); } diff --git a/src/cdi_analytic/test/tstAnalytic_Gray_Opacity.cc b/src/cdi_analytic/test/tstAnalytic_Gray_Opacity.cc index 30ff708042..1c98fec0aa 100644 --- a/src/cdi_analytic/test/tstAnalytic_Gray_Opacity.cc +++ b/src/cdi_analytic/test/tstAnalytic_Gray_Opacity.cc @@ -234,7 +234,7 @@ void CDI_test(rtt_dsxx::UnitTest &ut) { { std::vector params(amodel->get_parameters()); - std::vector expectedValue(8); + std::vector expectedValue(11); expectedValue[0] = 0.0; expectedValue[1] = 100.0; expectedValue[2] = -3.0; @@ -243,6 +243,9 @@ void CDI_test(rtt_dsxx::UnitTest &ut) { expectedValue[5] = 1.0; expectedValue[6] = 1.0; expectedValue[7] = 1.0; + expectedValue[8] = 0.0; + expectedValue[9] = 0.0; + expectedValue[10] = 0.0; double const tol(1.0e-12);