rad_planck_mean.cc

Go to the documentation of this file.

 
 
#include "rad_planck_mean.h"
#include <iostream>
#include <cstdlib>
 
using namespace std;
 
 
const double rad_planck_mean::pmCoefs_H2O[6]   = {-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840E-5};
const double rad_planck_mean::pmCoefs_CO2[6]   = {18.741, -121.310, 273.500, -194.050, 56.310, -5.8169};
const double rad_planck_mean::pmCoefs_CO_lo[5] = {4.7869, -0.06953, 2.95775E-4, -4.25732E-7, 2.02894E-10};
const double rad_planck_mean::pmCoefs_CO_hi[5] = {10.09, -0.01183, 4.7753E-6, -5.87209E-10, -2.5334E-14};
const double rad_planck_mean::pmCoefs_CH4[5]   = {6.6334, -0.0035686, 1.6682E-08, 2.5611E-10, -2.6558E-14};
 
 
void rad_planck_mean::get_k_a_oneband(double       &kabs,
                                      double       &awts,
                                      const int    iband,
                                      const double T,
                                      const double P,
                                      const double fvsoot,
                                      const double xH2O,
                                      const double xCO2,
                                      const double xCO,
                                      const double xCH4){
 
    if(iband != 0) {
        cerr << "\n\n***** rad_planck_mean::get_k_a_oneband: iband should be zero since there is only one band in this model *****\n" << endl; 
        exit(0); 
    }
 
    //----- the curve fits are in the range 300<T<2500. If T is out of bounds either quit,
    //----- or just take T to be the limiting value. 
    //----- taking the limiting value is used below when the next few lines are commented.
 
    // if(T < 300.0 || T > 2500.0) {
    //     cerr << "\n\n***** rad_planck_mean::get_k_a_oneband: T is out of range 300-2500 K *****\n" << endl; 
    //     exit(0); 
    // }
 
    kabs = 0.0;
    awts = 1.0;
 
    double K;
    double Ti = (T<300) ? 1000.0/300 : ((T>2500) ? 1000.0/2500 : 1000.0 / T);
 
    //------------- H2O
 
    if(xH2O != 0.0){
        K = pmCoefs_H2O[0] + Ti*(pmCoefs_H2O[1] + Ti*(pmCoefs_H2O[2] + Ti*(pmCoefs_H2O[3] + Ti*(pmCoefs_H2O[4] + Ti*(pmCoefs_H2O[5])))));
        kabs += xH2O*P/101325.0*K;
    }
 
    //------------- CO2
 
    if(xCO2 != 0.0){
        K = pmCoefs_CO2[0] + Ti*(pmCoefs_CO2[1] + Ti*(pmCoefs_CO2[2] + Ti*(pmCoefs_CO2[3] + Ti*(pmCoefs_CO2[4] + Ti*(pmCoefs_CO2[5])))));
        kabs += xCO2*P/101325.0*K;
    }
 
    //------------- CO
 
    Ti = (T<300) ? 300 : ((T>2500) ? 2500 : T);
 
    if(xCO != 0.0) {
        if(T <= 750.0)
            K = pmCoefs_CO_lo[0] + Ti*(pmCoefs_CO_lo[1] + Ti*(pmCoefs_CO_lo[2] + Ti*(pmCoefs_CO_lo[3] + Ti*(pmCoefs_CO_lo[4]))));
        else
            K = pmCoefs_CO_hi[0] + Ti*(pmCoefs_CO_hi[1] + Ti*(pmCoefs_CO_hi[2] + Ti*(pmCoefs_CO_hi[3] + Ti*(pmCoefs_CO_hi[4]))));
        kabs += xCO*P/101325.0*K;
    }
 
    //------------- CH4
 
    if(xCH4 != 0.0){
        K = pmCoefs_CH4[0] + Ti*(pmCoefs_CH4[1] + Ti*(pmCoefs_CH4[2] + Ti*(pmCoefs_CH4[3] + Ti*(pmCoefs_CH4[4]))));
        kabs += xCH4*P/101325.0*K;
    }
 
    //------------- soot
 
    if(fvsoot != 0.0){
        kabs += 1817 * fvsoot*Ti;       // 1817 = 3.72*csoot/C2.
    }
}
 
 
void rad_planck_mean::get_k_a(vector<double> &kabs,
                              vector<double> &awts,
                              const double   T,
                              const double   P,
                              const double   fvsoot,
                              const double   xH2O,
                              const double   xCO2,
                              const double   xCO,
                              const double   xCH4){
 
    double k;
    double a;
 
    get_k_a_oneband(k, a, 0, T, P, fvsoot, xH2O, xCO2, xCO, xCH4);
 
    kabs.resize(1);  kabs[0] = k;
    awts.resize(1);  awts[0] = a; 
}