initial commit of sergio's fvDOM model

src/thermophysicalModels/radiation/Make/files

@@ -1,13 +1,16 @@
+
 /* Radiation constants */
 radiationConstants/radiationConstants.C
 
-
 /* Radiation model */
 radiationModel/radiationModel/radiationModel.C
 radiationModel/radiationModel/newRadiationModel.C
 radiationModel/noRadiation/noRadiation.C
 radiationModel/P1/P1.C
-
+radiationModel/fvDOM/fvDOM/fvDOM.C
+radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
+radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.C
 
 /* Scatter model */
 submodels/scatterModel/scatterModel/scatterModel.C
@@ -21,11 +24,13 @@ submodels/absorptionEmissionModel/absorptionEmissionModel/newAbsorptionEmissionM
 submodels/absorptionEmissionModel/noAbsorptionEmission/noAbsorptionEmission.C
 submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.C
 submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmission.C
+submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 
 
 /* Boundary conditions */
 derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
-
+derivedFvPatchFields/GreyDiffusiveRadiation/GreyDiffusiveRadiationMixedFvPatchScalarField.C
-
+derivedFvPatchFields/WideBandDiffusiveRadiation/WideBandDiffusiveRadiationMixedFvPatchScalarField.C
 LIB = $(FOAM_LIBBIN)/libradiation

src/thermophysicalModels/radiation/Make/options

@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude \
+    -I radiationModel/fvDOM/interpolationLookUpTable
 
 LIB_LIBS = \
     -lfiniteVolume

src/thermophysicalModels/radiation/derivedFvPatchFields/GreyDiffusiveRadiation/GreyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -0,0 +1,312 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "GreyDiffusiveRadiationMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+#include "radiationModel.H"
+#include "fvDOM.H"
+#include "radiationConstants.H"
+#include "mathematicalConstants.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::
+GreyDiffusiveRadiationMixedFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_("undefined"),
+    emissivity_(0.0),
+    myRayIndex_(0),
+    myWaveLengthIndex_(0),
+    myRayIsInit_(-1),
+    qr_(0)
+{
+    refValue() = 0.0;
+    refGrad() = 0.0;
+    valueFraction() = 1.0;
+    qr_.setSize(p.size());
+}
+
+
+Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::
+GreyDiffusiveRadiationMixedFvPatchField
+(
+    const GreyDiffusiveRadiationMixedFvPatchField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+{}
+
+
+Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::
+GreyDiffusiveRadiationMixedFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_(dict.lookup("T")),
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    myRayIndex_(0),
+    myWaveLengthIndex_(0),
+    myRayIsInit_(-1),
+    qr_(0)
+{
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    refValue() = emissivity_*4.0*radiation::sigmaSB.value()*pow4(Tp) /
+                 Foam::mathematicalConstant::pi;
+    refGrad() = 0.0;
+
+    qr_.setSize(p.size());
+
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(refValue());
+    }
+}
+
+
+Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::
+GreyDiffusiveRadiationMixedFvPatchField
+(
+    const GreyDiffusiveRadiationMixedFvPatchField& ptf
+)
+:
+    mixedFvPatchScalarField(ptf),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+{}
+
+
+Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::
+GreyDiffusiveRadiationMixedFvPatchField
+(
+    const GreyDiffusiveRadiationMixedFvPatchField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(ptf, iF),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    scalarField::autoMap(m);
+
+}
+
+
+void Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    mixedFvPatchScalarField::rmap(ptf, addr);
+
+//    const GreyDiffusiveRadiationMixedFvPatchField& mrptf =
+        refCast<const GreyDiffusiveRadiationMixedFvPatchField>(ptf);
+}
+
+
+void Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    const radiationModel& rad =
+            db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& Dom(refCast<const fvDOM>(rad));
+
+    const label patchi = patch().index();
+
+    if(Dom.lambdaj() == 1)
+    {
+        if (myRayIsInit_ == -1)
+        {
+            for(label i=0; i < Dom.Ni() ; i++)
+            {
+                 for(label j=0; j < Dom.lambdaj() ; j++)
+                {
+                    const volScalarField& radiationField =
+                                        Dom.RadIntRayiLambdaj(i,j);
+                    if (&(radiationField.internalField()) ==
+                        &dimensionedInternalField())
+                    {
+                        myRayIndex_ = i;
+                        myWaveLengthIndex_ = j;
+                        myRayIsInit_ = 0.;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::radiation::"
+            "GreyDiffusiveRadiationMixedFvPatchField::"
+            "updateCoeffs"
+        )   << " a grey boundary condition is used with a non-grey"
+            << "absorption model"
+            << exit(FatalError);
+    }
+
+    vectorField n = patch().Sf()/patch().magSf();
+
+    scalarField& Iw = *(this);
+
+    qr_ =  Iw *(-n & Dom.RadIntRay(myRayIndex_).Di());
+
+    Dom.RadIntRay(myRayIndex_).add(qr_,patchi);
+
+    forAll(Iw, faceI)
+    {
+
+        scalar Ir = 0.0;
+
+        for(label i=0; i < Dom.Ni() ; i++) //
+        {
+            const vector& si = Dom.RadIntRay(i).Si();
+
+            const scalarField&  Iface = Dom.RadIntRay(i).Ilambdaj
+            (
+                myWaveLengthIndex_
+            ).boundaryField()[patch().index()];
+
+            scalar InOut = -n[faceI] & si;
+
+            if (InOut < 0.) // qin into the wall
+            {
+                const vector& di = Dom.RadIntRay(i).Di();
+                Ir += Iface[faceI]*mag(n[faceI] & di);
+            }
+        }
+
+        const vector& mySi = Dom.RadIntRay(myRayIndex_).Si();
+
+        scalar InOut = -n[faceI] & mySi;
+
+        if (InOut > 0.) //direction out of the wall
+        {
+            refGrad()[faceI] = 0.0;
+            valueFraction()[faceI] = 1.0;
+            refValue()[faceI] = ((1. - emissivity_) * Ir +
+                    emissivity_*radiation::sigmaSB.value()*pow4(Tp[faceI])) /
+                    Foam::mathematicalConstant::pi;
+
+        }
+        else if (InOut < 0.) //direction into the wall
+        {
+            valueFraction()[faceI] = 0.0;
+            refGrad()[faceI] = 0.0;
+            refValue()[faceI] = 0.0; //not used
+        }
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::radiation::GreyDiffusiveRadiationMixedFvPatchField::write(Ostream&
+os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("T") << TName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
+
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        GreyDiffusiveRadiationMixedFvPatchField
+    );
+}
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/GreyDiffusiveRadiation/GreyDiffusiveRadiationMixedFvPatchScalarField.H

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::GreyDiffusiveRadiationMixedFvPatchField
+
+Description
+    Radiation temperature specified
+
+SourceFiles
+    GreyDiffusiveRadiationMixedFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef GreyDiffusiveRadiationMixedFvPatchField_H
+#define GreyDiffusiveRadiationMixedFvPatchField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+/*---------------------------------------------------------------------------*\
+        Class GreyDiffusiveRadiationMixedFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+class GreyDiffusiveRadiationMixedFvPatchField
+:
+    public mixedFvPatchScalarField
+{
+
+    // Private data
+
+        //- Name of temperature field
+        word TName_;
+
+        //- Emissivity
+        scalar emissivity_;
+
+        //- Direction index
+        label myRayIndex_;
+
+        //- Direction index
+        label myWaveLengthIndex_;
+
+        //- Init ray flag
+        label myRayIsInit_;
+
+        //- Radiative heat flux on walls.
+        scalarField qr_;
+
+public:
+
+    //- Runtime type information
+    TypeName("GreyDiffusiveRadiation");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        GreyDiffusiveRadiationMixedFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        GreyDiffusiveRadiationMixedFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given GreyDiffusiveRadiationMixedFvPatchField
+        //  onto a new patch
+        GreyDiffusiveRadiationMixedFvPatchField
+        (
+            const GreyDiffusiveRadiationMixedFvPatchField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        GreyDiffusiveRadiationMixedFvPatchField
+        (
+            const GreyDiffusiveRadiationMixedFvPatchField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new GreyDiffusiveRadiationMixedFvPatchField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        GreyDiffusiveRadiationMixedFvPatchField
+        (
+            const GreyDiffusiveRadiationMixedFvPatchField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new GreyDiffusiveRadiationMixedFvPatchField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+        //- Return the temperature field name
+            const word& TName() const
+            {
+                return TName_;
+            }
+
+        //- Return reference to the temperature field name to allow
+            //  adjustment
+            word& TName()
+            {
+                return TName_;
+            }
+
+            //- Return the emissivity
+            const scalar& emissivity() const
+            {
+                return emissivity_;
+            }
+
+            //- Return reference to the emissivity to allow adjustment
+            scalar& emissivity()
+            {
+                return emissivity_;
+            }
+
+
+            //- Return heat flux on the boundary
+            const scalarField& qr() const
+            {
+                return qr_;
+            }
+
+            // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/WideBandDiffusiveRadiation/WideBandDiffusiveRadiationMixedFvPatchScalarField.C

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WideBandDiffusiveRadiationMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+#include "radiationModel.H"
+#include "fvDOM.H"
+#include "wideBandAbsorptionEmission.H"
+#include "radiationConstants.H"
+#include "mathematicalConstants.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::
+WideBandDiffusiveRadiationMixedFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_("undefined"),
+    emissivity_(0.0),
+    myRayIndex_(0),
+    myWaveLengthIndex_(0),
+    myRayIsInit_(-1),
+    qr_(0)
+{
+    refValue() = 0.0;
+    refGrad() = 0.0;
+    valueFraction() = 1.0;
+    qr_.setSize(p.size());
+}
+
+
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::
+WideBandDiffusiveRadiationMixedFvPatchField
+(
+    const WideBandDiffusiveRadiationMixedFvPatchField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+{}
+
+
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::
+WideBandDiffusiveRadiationMixedFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    TName_(dict.lookup("T")),
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    myRayIndex_(0),
+    myWaveLengthIndex_(0),
+    myRayIsInit_(-1),
+    qr_(0)
+{
+    const scalarField& Tp =
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+    refValue() = emissivity_*4.0*radiation::sigmaSB.value()*pow4(Tp) /
+                 Foam::mathematicalConstant::pi;
+    refGrad() = 0.0;
+
+    qr_.setSize(p.size());
+
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(refValue());
+    }
+}
+
+
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::
+WideBandDiffusiveRadiationMixedFvPatchField
+(
+    const WideBandDiffusiveRadiationMixedFvPatchField& ptf
+)
+:
+    mixedFvPatchScalarField(ptf),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+{}
+
+
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::
+WideBandDiffusiveRadiationMixedFvPatchField
+(
+    const WideBandDiffusiveRadiationMixedFvPatchField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(ptf, iF),
+    TName_(ptf.TName_),
+    emissivity_(ptf.emissivity_),
+    myRayIndex_(ptf.myRayIndex_),
+    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
+    myRayIsInit_(ptf.myRayIsInit_),
+    qr_(ptf.qr_)
+
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    scalarField::autoMap(m);
+
+}
+
+
+void Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    mixedFvPatchScalarField::rmap(ptf, addr);
+
+//    const GreyDiffusiveRadiationMixedFvPatchField& mrptf =
+        refCast<const WideBandDiffusiveRadiationMixedFvPatchField>(ptf);
+}
+
+
+void
+Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::updateCoeffs
+()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const radiationModel& rad =
+            db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& Dom(refCast<const fvDOM>(rad));
+
+    const label patchi = patch().index();
+
+    if(Dom.lambdaj() > 1)
+    {
+        if (myRayIsInit_ == -1)
+        {
+            for(label i=0; i < Dom.Ni() ; i++)
+            {
+                for(label j=0; j < Dom.lambdaj() ; j++)
+                {
+                    const volScalarField& radiationField =
+                                        Dom.RadIntRayiLambdaj(i,j);
+                    if (&(radiationField.internalField()) ==
+                        &dimensionedInternalField())
+                    {
+                        myRayIndex_ = i;
+                        myWaveLengthIndex_ = j;
+                        myRayIsInit_ = 0.;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::radiation::"
+            "WideBandDiffusiveRadiationMixedFvPatchScalarField::"
+            "updateCoeffs"
+        )   << " a Non-grey boundary condition is used with a grey"
+            << "absorption model"
+            << exit(FatalError);
+    }
+
+    vectorField n = patch().Sf()/patch().magSf();
+
+    scalarField& Iw = *(this);
+
+    qr_ =  Iw *(-n & Dom.RadIntRay(myRayIndex_).Di());
+
+    Dom.RadIntRay(myRayIndex_).add(qr_,patchi);
+
+    const scalarField Eb =
+            Dom.blackBody().bj(myWaveLengthIndex_).boundaryField()[patchi];
+
+    forAll(Iw, faceI)
+    {
+
+        scalar Ir = 0.0;
+        for(label i=0; i < Dom.Ni() ; i++) //
+        {
+            const vector& si = Dom.RadIntRay(i).Si();
+
+            const scalarField& Iface = Dom.RadIntRay(i).Ilambdaj
+            (
+                myWaveLengthIndex_
+            ).boundaryField()[patch().index()];
+
+            scalar InOut = -n[faceI] & si;
+
+            if (InOut < 0.) // qin into the wall
+            {
+                const vector& di = Dom.RadIntRay(i).Di();
+                Ir = Ir + Iface[faceI]*mag(n[faceI] & di);
+            }
+        }
+
+
+        const vector& mySi = Dom.RadIntRay(myRayIndex_).Si();
+
+        scalar InOut = -n[faceI] & mySi;
+
+        if (InOut > 0.) //direction out of the wall
+        {
+            refGrad()[faceI] = 0.0;
+            valueFraction()[faceI] = 1.0;
+            refValue()[faceI] = ((1. - emissivity_) * Ir +
+                    emissivity_* Eb[faceI]) / Foam::mathematicalConstant::pi;
+        }
+        else if (InOut < 0.) //direction into the wall
+        {
+            valueFraction()[faceI] = 0.0;
+            refGrad()[faceI] = 0.0;
+            refValue()[faceI] = 0.0; //not used
+        }
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::radiation::WideBandDiffusiveRadiationMixedFvPatchField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("T") << TName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        WideBandDiffusiveRadiationMixedFvPatchField
+    );
+}
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/WideBandDiffusiveRadiation/WideBandDiffusiveRadiationMixedFvPatchScalarField.H

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::WideBandDiffusiveRadiationMixedFvPatchField
+
+Description
+    Radiation temperature specified
+
+SourceFiles
+    WideBandDiffusiveRadiationMixedFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WideBandDiffusiveRadiationMixedFvPatchField_H
+#define WideBandDiffusiveRadiationMixedFvPatchField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+/*---------------------------------------------------------------------------*\
+        Class WideBandDiffusiveRadiationMixedFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+class WideBandDiffusiveRadiationMixedFvPatchField
+:
+    public mixedFvPatchScalarField
+{
+
+    // Private data
+
+        //- Name of temperature field
+        word TName_;
+
+        //- Emissivity
+        scalar emissivity_;
+
+        //- Direction index
+        label myRayIndex_;
+
+        //- Direction index
+        label myWaveLengthIndex_;
+
+        //- Init ray flag
+        label myRayIsInit_;
+
+        //- Radiative heat flux on walls.
+        scalarField qr_;
+
+public:
+
+    //- Runtime type information
+    TypeName("WideBandDiffusiveRadiation");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        WideBandDiffusiveRadiationMixedFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        WideBandDiffusiveRadiationMixedFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given GreyDiffusiveRadiationMixedFvPatchField
+        //  onto a new patch
+        WideBandDiffusiveRadiationMixedFvPatchField
+        (
+            const WideBandDiffusiveRadiationMixedFvPatchField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        WideBandDiffusiveRadiationMixedFvPatchField
+        (
+            const WideBandDiffusiveRadiationMixedFvPatchField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new WideBandDiffusiveRadiationMixedFvPatchField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        WideBandDiffusiveRadiationMixedFvPatchField
+        (
+            const WideBandDiffusiveRadiationMixedFvPatchField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new WideBandDiffusiveRadiationMixedFvPatchField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+        //- Return the temperature field name
+            const word& TName() const
+            {
+                return TName_;
+            }
+
+        //- Return reference to the temperature field name to allow
+            //  adjustment
+            word& TName()
+            {
+                return TName_;
+            }
+
+            //- Return the emissivity
+            const scalar& emissivity() const
+            {
+                return emissivity_;
+            }
+
+            //- Return reference to the emissivity to allow adjustment
+            scalar& emissivity()
+            {
+                return emissivity_;
+            }
+
+
+            //- Return heat flux on the boundary
+            const scalarField& qr() const
+            {
+                return qr_;
+            }
+
+            // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/P1/P1.C

@@ -50,6 +50,9 @@ namespace Foam
     }
 }
 
+//  Radiation solver iteration
+Foam::label Foam::radiation::P1::iterRadId = pTraits<label>::one;
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -135,10 +138,12 @@ bool Foam::radiation::P1::read()
 
 void Foam::radiation::P1::correct()
 {
-    if (!radiation_)
+    if (!radiation_ || !(iterRadId == nFlowIterPerRadIter_))
     {
+        iterRadId++;
         return;
     }
+
     a_ = absorptionEmission_->a();
     e_ = absorptionEmission_->e();
     E_ = absorptionEmission_->E();
@@ -166,6 +171,8 @@ void Foam::radiation::P1::correct()
      ==
       - 4.0*(e_*radiation::sigmaSB*pow4(T_) + E_)
     );
+
+    iterRadId = pTraits<label>::one;
 }
 
 

src/thermophysicalModels/radiation/radiationModel/P1/P1.H

@@ -86,6 +86,10 @@ class P1
 
 public:
 
+   // Static data members
+
+        static label iterRadId;
+
     //- Runtime type information
     TypeName("P1");
 

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "absorptionCoeffs.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::radiation::absorptionCoeffs::~absorptionCoeffs()
+{}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::radiation::absorptionCoeffs::absorptionCoeffs(Istream& is)
+:
+   Tcommon_(readScalar(is)),
+   Tlow_(readScalar(is)),
+   Thigh_(readScalar(is)),
+   invTemp_(readBool(is))
+{
+    for
+    (
+        label coefLabel=0;
+        absorptionCoeffs::nCoeffs_;
+        coefLabel++
+    )
+    {
+        is >> highACoeffs_[coefLabel];
+    }
+
+    for
+    (
+        label coefLabel=0;
+        absorptionCoeffs::nCoeffs_;
+        coefLabel++
+    )
+    {
+        is >> lowACoeffs_[coefLabel];
+    }
+}
+
+void Foam::radiation::absorptionCoeffs::checkT(const scalar T) const
+{
+    if (T <  Tlow_ || T > Thigh_)
+    {
+        FatalErrorIn
+        (
+            "absCoeffs::checkT(const scalar T) const"
+        )   << "attempt to use absCoeff"
+               " out of temperature range "
+            << Tlow_ << " -> " << Thigh_ << ";  T = " << T
+            << abort(FatalError);
+    }
+}
+
+
+const Foam::radiation::absorptionCoeffs::coeffArray&
+Foam::radiation::absorptionCoeffs::coeffs
+(
+    const scalar T
+) const
+{
+    checkT(T);
+
+    if (T < Tcommon_)
+    {
+        return lowACoeffs_;
+    }
+    else
+    {
+        return highACoeffs_;
+    }
+}
+
+void Foam::radiation::absorptionCoeffs::init
+(
+    const dictionary& dict
+)
+{
+    Tcommon_ = readScalar(dict.lookup("Tcommon"));
+    Tlow_ = readScalar(dict.lookup("Tlow"));
+    Thigh_ = readScalar(dict.lookup("Thigh"));
+    invTemp_ = readBool(dict.lookup("invTemp"));
+
+    dict.lookup("loTcoeffs") >> lowACoeffs_;
+    dict.lookup("hiTcoeffs") >> highACoeffs_;
+}
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/absorptionCoeffs/absorptionCoeffs.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::janafThermo
+
+Description
+    Absorption coefficients class used in greyMeanAbsorptionEmission and
+    wideBandAbsorptionEmission
+
+SourceFiles
+    absorptionCoeffs.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef absorptionCoeffs_H
+#define absorptionCoeffs_H
+
+#include "List.H"
+#include "IOstreams.H"
+#include "IOdictionary.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class absorptionCoeffs Declaration
+\*---------------------------------------------------------------------------*/
+
+class absorptionCoeffs
+{
+
+public:
+
+        static const int nCoeffs_ = 6;
+        typedef FixedList<scalar, nCoeffs_> coeffArray;
+
+private:
+
+    // Private data
+
+        // Temperature limits of applicability of functions
+        scalar Tcommon_;
+
+        scalar Tlow_;
+
+        scalar Thigh_;
+
+        // Polynomio using inverse temperatures
+        bool invTemp_;
+
+        coeffArray highACoeffs_;
+        coeffArray lowACoeffs_;
+
+
+    // Private member functions
+
+        //- Check given temperature is within the range of the fitted coeffs
+        void checkT(const scalar T) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from Istream
+        absorptionCoeffs(Istream&);
+
+        // Null constructor
+        absorptionCoeffs()
+        {}
+
+    // Destructor
+
+        ~absorptionCoeffs();
+
+    // member functions
+
+        //- Return the coefficients corresponding to the given temperature
+        const coeffArray& coeffs(const scalar T) const;
+
+        // Init from a dictionary
+        void init(const dictionary&);
+
+        // Init from an Istram
+        inline void init(Istream&)
+        {
+             absorptionCoeffs(Istream);
+        }
+
+
+    // Acces Functions
+        inline bool invTemp() const
+        {
+            return  invTemp_;
+        }
+
+        inline scalar Tcommon() const
+        {
+            return  Tcommon_;
+        }
+
+        inline scalar Tlow() const
+        {
+            return  Tlow_;
+        }
+
+        inline scalar Thigh() const
+        {
+            return  Thigh_;
+        }
+
+        inline const coeffArray& highACoeffs() const
+        {
+            return  highACoeffs_;
+        }
+
+        inline const coeffArray& lowACoeffs() const
+        {
+            return  lowACoeffs_;
+        }
+
+};
+
+} // End namespace Foam
+
+} // End namespace radiation
+
+#endif
+

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "blackBodyEmission.H"
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::blackBodyEmission::blackBodyEmission
+(
+    const fileName& fn,
+    const word& instance,
+    label lambdaj,
+    const volScalarField& T
+)
+:blackBodyEmissiveTable_(fn, instance, T.mesh()),
+C1_("C1",dimensionSet(1,4,3,0,0,0,0),3.7419e-16),
+C2_("C2",dimensionSet(0,1,0,1,0,0,0),14.388e-6),
+bj_(0),
+T_(T)
+{
+    bj_.setSize(lambdaj);
+    for(label i=0; i < lambdaj; i++)
+    {
+        bj_.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "bj_" + Foam::name(i) ,
+                    T.mesh().time().timeName(),
+                    T.mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                radiation::sigmaSB*pow4(T)
+            )
+        );
+
+    }
+}
+
+// * * * * * * *  Destructor * * * * * * * * * * * * * * * * * * * * * * * * //
+Foam::radiation::blackBodyEmission::~blackBodyEmission()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+Foam::scalar Foam::radiation::blackBodyEmission::flambdaT
+(
+    const scalar lambdaT
+) const
+{
+    return  blackBodyEmissiveTable_.LookUp(lambdaT*1e6)[1];
+}
+
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::blackBodyEmission::EbDeltaLambdaT
+(
+    const volScalarField& T,
+    const Vector2D<scalar>& band
+) const
+{
+    tmp<volScalarField> Eb
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Eb",
+                T.mesh().time().timeName(),
+                T.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            radiation::sigmaSB*pow4(T)
+        )
+    );
+
+
+    if (band == Vector2D<scalar>::one)
+    {
+        return Eb;
+    }
+    else
+    {
+        forAll(T,i)
+        {
+            scalar T1 = flambdaT(band[1]*T[i]);
+            scalar T2 = flambdaT(band[0]*T[i]);
+            dimensionedScalar flambdaDelta
+            (
+                "flambdaDelta",
+                dimless,
+                T1 - T2
+            );
+            Eb()[i] = Eb()[i]*flambdaDelta.value();
+        }
+        return Eb;
+    }
+
+}
+
+
+void Foam::radiation::blackBodyEmission::correct
+(
+    label j,
+    const Vector2D<scalar>& band
+)
+{
+    bj_[j] = EbDeltaLambdaT(T_, band);
+}
+
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::blackBodyEmission
+
+Description
+    Class Black Body Emission Intensity.
+
+SourceFiles
+    blackBodyEmission.C
+
+
+\*---------------------------------------------------------------------------*/
+#ifndef blackModyEmission_H
+#define blackModyEmission_H
+
+#include "volFields.H"
+#include "dimensionedScalar.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+#include "interpolationLookUpTable.H"
+#include "Vector2D.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+/*---------------------------------------------------------------------------*\
+                           Class blackBodyEmission Declaration
+\*---------------------------------------------------------------------------*/
+class blackBodyEmission
+{
+    // Private data
+    mutable interpolationLookUpTable<scalar> blackBodyEmissiveTable_;
+
+    scalar flambdaT(const scalar lambdaT) const;
+
+    const dimensionedScalar C1_;
+
+    const dimensionedScalar C2_;
+
+    // Ptr List for enregy black body emission
+    PtrList<volScalarField> bj_;
+
+    // T
+    const volScalarField& T_;
+
+public:
+
+
+    // Constructors
+    blackBodyEmission
+    (
+        const fileName& fn,
+        const word& instance,
+        label lambdaj,
+        const volScalarField& T
+    );
+
+    // Destructor
+    ~blackBodyEmission();
+
+
+    // - Spectral emission for the black body at T and lambda
+    inline dimensionedScalar EblambdaT
+    (
+        const dimensionedScalar T,
+        const scalar lambda
+    ) const
+    {
+        return (C1_/(pow5(lambda)*(exp(C2_/(lambda*T))-1)));
+    }
+
+    // Integral Energy at T from lambda1 to lambda2
+
+    tmp<Foam::volScalarField> EbDeltaLambdaT //dimensionedScalar
+    (
+        const volScalarField& T,
+        const Vector2D<scalar>& band
+    ) const;
+
+    // Edit
+
+        // Update black body emission
+        void correct(label j, const Vector2D<scalar>& band);
+
+    // Acces members
+
+        //   Black body spectrum
+       inline const volScalarField& bj(label i) const
+       {
+            return bj_[i];
+       }
+
+};
+
+}
+}
+
+#endif
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmissivePower

@@ -0,0 +1,174 @@
+171
+(
+1000    0.00032
+1100    0.00091
+1200    0.00213
+1300    0.00432
+1400    0.00779
+1500    0.01285
+1600    0.01972
+1700    0.02853
+1800    0.03934
+1900    0.05210
+2000    0.06672
+2100    0.08305
+2200    0.10088
+2300    0.12002
+2400    0.14025
+2500    0.16135
+2600    0.18311
+2700    0.20535
+2800    0.22788
+2900    0.25055
+3000    0.27322
+3100    0.29576
+3200    0.31809
+3300    0.34009
+3400    0.36172
+3500    0.38290
+3600    0.40359
+3700    0.42375
+3800    0.44336
+3900    0.46240
+4000    0.48085
+4100    0.49872
+4200    0.51599
+4300    0.53267
+4400    0.54877
+4500    0.56429
+4600    0.57925
+4700    0.59366
+4800    0.60753
+4900    0.62088
+5000    0.63372
+5100    0.64606
+5200    0.65794
+5300    0.66935
+5400    0.68033
+5500    0.69087
+5600    0.70101
+5700    0.71076
+5800    0.72012
+5900    0.72913
+6000    0.73778
+6100    0.74610
+6200    0.75410
+6300    0.76180
+6400    0.76920
+6500    0.77631
+6600    0.78316
+6700    0.78975
+6800    0.79609
+6900    0.80219
+7000    0.80807
+7100    0.81373
+7200    0.81918
+7300    0.82443
+7400    0.82949
+7500    0.83436
+7600    0.83906
+7700    0.84359
+7800    0.84796
+7900    0.85218
+8000    0.85625
+8100    0.86017
+8200    0.86396
+8300    0.86762
+8400    0.87115
+8500    0.87456
+8600    0.87786
+8700    0.88105
+8800    0.88413
+8900    0.88711
+9000    0.88999
+9100   0.89277
+9200   0.89547
+9300   0.89807
+9400   0.90060
+9500   0.90304
+9600   0.90541
+9700   0.90770
+9800   0.90992
+9900   0.91207
+10000   0.91415
+10200   0.91813
+10400   0.92188
+10600   0.92540
+10800   0.92872
+11000   0.93184
+11200   0.93479
+11400   0.93758
+11600   0.94021
+11800   0.94270
+12000   0.94505
+12200   0.94728
+12400   0.94939
+12600   0.95139
+12800   0.95329
+13000   0.95509
+13200   0.95680
+13400   0.95843
+13600   0.95998
+13800   0.96145
+14000   0.96285
+14200   0.96418
+14400   0.96546
+14600   0.96667
+14800   0.96783
+15000   0.96893
+15200   0.96999
+15400   0.97100
+15600   0.97196
+15800   0.97288
+16000   0.97377
+16200   0.97461
+16400   0.97542
+16600   0.97620
+16800   0.97694
+17000   0.97765
+17200   0.97834
+17400   0.97899
+17600   0.97962
+17800   0.98023
+18000   0.98080
+18200   0.98137
+18400   0.98191
+18600   0.98243
+18900   0.98293
+19000   0.98340
+19200   0.98387
+19400   0.98431
+19600   0.98474
+19800   0.98515
+20000   0.98555
+21000   0.98735
+22000   0.98886
+23000   0.99014
+24000   0.99123
+25000   0.99217
+26000   0.99297
+27000   0.99367
+28000   0.99429
+29000   0.99482
+30000   0.99529
+31000   0.99571
+32000   0.99607
+33000   0.99640
+34000   0.99669
+35000   0.99695
+35000   0.99719
+36000   0.99740
+37000   0.99759
+38000   0.99776
+39000   0.99792
+40000   0.99806
+41000   0.99819
+42000   0.99831
+43000   0.99842
+44000   0.99851
+45000   0.99861
+46000   0.99869
+47000   0.99877
+48000   0.99884
+49000   0.99890
+)

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C

@@ -0,0 +1,353 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvDOM.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+
+#include "absorptionEmissionModel.H"
+#include "scatterModel.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(fvDOM, 0);
+
+        addToRunTimeSelectionTable
+        (
+            radiationModel,
+            fvDOM,
+            dictionary
+        );
+    }
+}
+
+//  Radiation solver iterator counter
+Foam::label Foam::radiation::fvDOM::iterRadId = pTraits<label>::one;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
+:
+    radiationModel(typeName, T),
+    G_
+    (
+        IOobject
+        (
+            "G",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qr_
+    (
+        IOobject
+        (
+            "Qr",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+    ),
+    a_
+    (
+        IOobject
+        (
+            "a",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("a", dimless/dimLength, 0.0)
+    ),
+    aj_(0),
+    e_
+    (
+        IOobject
+        (
+            "e",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("a", dimless/dimLength, 0.0)
+    ),
+    E_
+    (
+        IOobject
+        (
+            "E",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+    ),
+    Ntheta_(readLabel(radiationModelCoeffs_.lookup("Ntheta"))),
+    Nphi_(readLabel(radiationModelCoeffs_.lookup("Nphi"))),
+    Ni_(0),
+    lambdaj_(absorptionEmission_->nBands()),
+    blackBody_(fileName("blackBodyEmissivePower"), "constant", lambdaj_, T)
+{
+    aj_.setSize(lambdaj_);
+    if (mesh_.nSolutionD() == 3)    //3D
+    {
+        RadIntRay_.setSize(4.* Nphi_* Ntheta_);
+        Ni_ = 4. * Nphi_ * Ntheta_;
+        scalar deltaPhi = mathematicalConstant::pi / (2. * Nphi_);
+        scalar deltaTheta = mathematicalConstant::pi / Ntheta_;
+        label i = 0;
+        for(label n = 1 ; n <= Ntheta_ ; n++)
+        {
+            for(label m = 1 ; m <= 4*Nphi_ ; m++)
+            {
+                scalar thetai = (2.*n - 1.)*deltaTheta/2.;
+                scalar phii = (2.*m - 1.)*deltaPhi/2.;
+                RadIntRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        phii, thetai, deltaPhi,
+                        deltaTheta, lambdaj_, mesh_,
+                        absorptionEmission_, blackBody_
+                    )
+                );
+                i++;
+            }
+        }
+    }
+    else
+    {
+        if (mesh_.nSolutionD() == 2)    //2D (X & Y)
+        {
+            scalar thetai = mathematicalConstant::pi/2.;
+            scalar deltaTheta = mathematicalConstant::pi;
+            RadIntRay_.setSize(4.* Nphi_);
+            Ni_ = 4. * Nphi_;
+            scalar deltaPhi = mathematicalConstant::pi / (2. * Nphi_);
+            label i = 0;
+            for(label m = 1 ; m <= 4*Nphi_ ; m++)
+            {
+                scalar phii = (2.*m - 1.)*deltaPhi/2.;
+                RadIntRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        phii, thetai, deltaPhi,
+                        deltaTheta, lambdaj_, mesh_,
+                        absorptionEmission_, blackBody_
+                    )
+                );
+                i++;
+            }
+        }
+        else    //1D (X)
+        {
+            scalar thetai = mathematicalConstant::pi/2.;
+            scalar deltaTheta = mathematicalConstant::pi;
+            RadIntRay_.setSize(2);
+            Ni_ = 2.;
+            scalar deltaPhi = mathematicalConstant::pi;
+            label i = 0;
+            for(label m = 1 ; m <= 2 ; m++)
+            {
+                scalar phii = (2*m - 1.)*deltaPhi/2.;
+                RadIntRay_.set
+                (
+                    i,
+                    new radiativeIntensityRay
+                    (
+                        phii, thetai, deltaPhi,
+                        deltaTheta, lambdaj_, mesh_,
+                        absorptionEmission_, blackBody_
+                    )
+                );
+                i++;
+            }
+
+        }
+    }
+
+
+// Construct absorption for wave length
+    for(label i=0; i < lambdaj_; i++)
+    {
+        volScalarField* volPtr= new volScalarField
+        (
+            IOobject
+            (
+                "aj_" + Foam::name(i) ,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            a_
+        );
+
+        aj_.set(i, volPtr);
+    }
+
+}
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::fvDOM::~fvDOM()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::radiation::fvDOM::read()
+{
+    if (radiationModel::read())
+    {
+        // nothing to read
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::radiation::fvDOM::correct()
+{
+    if (!radiation_ || !(iterRadId == nFlowIterPerRadIter_))
+    {
+        iterRadId++;
+        return;
+    }
+
+    absorptionEmission_->correct(a_, aj_);
+
+    updateBlackBodyEmission();
+
+    scalar maxResidual = 0;
+    scalar convergenceCriterion = 0;
+    radiationModelCoeffs_.readIfPresent("convergence", convergenceCriterion);
+    label radIter = 0;
+    do
+    {
+         radIter ++;
+        for(label i = 0; i < Ni_; i++) //
+        {
+            maxResidual = 0;
+            scalar maxBandResidual = RadIntRay_[i].correct(this);
+            maxResidual = max(maxBandResidual, maxResidual);
+        }
+
+        Info << "Radiation solver Iter: " <<  radIter << endl;
+
+    }while(maxResidual > convergenceCriterion);
+
+    updateG();
+
+    iterRadId = pTraits<label>::one;
+}
+
+Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
+{
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Rp",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            4.0*a_*radiation::sigmaSB //absorptionEmission_->a()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::radiation::fvDOM::Ru() const
+{
+
+    const DimensionedField<scalar, volMesh>& G =
+        G_.dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> E =
+        absorptionEmission_->ECont()().dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> a =
+        a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
+
+    return  a*G - 4.0*E;
+}
+
+void Foam::radiation::fvDOM::updateBlackBodyEmission()
+{
+
+    for(label j=0; j < lambdaj_; j++)
+    {
+          blackBody_.correct(j, absorptionEmission_->bands(j));
+    }
+
+}
+
+void Foam::radiation::fvDOM::updateG()
+{
+
+    G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+    Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+
+    for(label i = 0; i < Ni_; i++)
+    {
+        RadIntRay_[i].addIntensity();
+        G_ +=  RadIntRay_[i].I()* RadIntRay_[i].omegai();
+        Qr_ += RadIntRay_[i].Qri();
+    }
+}
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H

@@ -0,0 +1,257 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::fvDOM
+
+Description
+
+    Finite Volume Discrete Ordinary Method. It solves the RTE equation for n
+    directions in a participating media. It does not consider scatter.
+
+    Available absorption models:
+                                greyMeanAbsoprtionEmission
+                                wideBandAbsorptionEmission
+
+    i.e. dictionary
+    fvDOMCoeffs
+    {
+        Nphi    1;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+        Ntheta  2;          // polar angles in P1 (from Z to X-Y plane)
+        convergence 1e-4;   //convergence criteria for radiation iteration
+    }
+
+    nFlowIterPerRadIter   1; // Number of flow iterations per radiation
+                                iteration
+
+    The total number of solid angles is  4 * Nphi * Ntheta.
+
+    In 1D the direction of the rays is X (Nphi and Ntheta are ignored)
+    In 2D the direction of the rays is on X-Y plane (only Nphi is considered)
+    In 3D (Nphi and Ntheta are considered)
+
+SourceFiles
+    fvDOM.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiationModelfvDOM_H
+#define radiationModelfvDOM_H
+
+#include "radiativeIntensityRay.H"
+#include "radiationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class fvDOM Declaration
+\*---------------------------------------------------------------------------*/
+
+class fvDOM
+:
+    public radiationModel
+{
+
+    // Private data
+
+        //- Incident radiation  [W/m2]
+        volScalarField G_;
+
+        //- Total Radiative heat flux [W/m2]
+        volScalarField Qr_;
+
+        //- Total Absorption coefficient  [1/m]
+        volScalarField a_;
+
+        //- wavelength total Absorption coefficient [1/m]
+        PtrList<volScalarField> aj_;
+
+        //- Total Emission coefficient [1/m]
+        volScalarField e_;
+
+        //- Emission contribution [Kg/m/s^3]
+        volScalarField E_;
+
+        //- Number of solid angles in theta
+        label Ntheta_;
+
+        //- Number of solid angles in phi
+        label Nphi_ ;
+
+        //- Total number of directions
+        label Ni_;
+
+        //- Number of wavelength bands
+        label lambdaj_;
+
+        //- Black Body
+        blackBodyEmission blackBody_;
+
+        //- List of Pointers to RadiativeIntensityRay
+        PtrList<radiativeIntensityRay> RadIntRay_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        fvDOM(const fvDOM&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fvDOM&);
+
+        //- Update Absorption Coefficients
+//        void updateAbsorptionCoeffs(void);
+
+        //- Update Black Body Emissiom
+        void updateBlackBodyEmission(void);
+
+public:
+
+    // Static data members
+
+        static label iterRadId;
+
+    //- Runtime type information
+        TypeName("fvDOM");
+
+
+    // Constructors
+
+        //- Construct from components
+        fvDOM(const volScalarField& T);
+
+
+    // Destructor
+
+        ~fvDOM();
+
+
+    // Member functions
+
+        // Edit
+
+        //- Update radiationSource varible
+        void correct();
+
+        //- Read radiationProperties dictionary
+        bool read();
+
+        //- Update G and calculate total heat flux on boundary
+        void updateG();
+
+        //- Source term component (for power of T^4)
+        virtual tmp<volScalarField> Rp() const;
+
+        //- Source term component (constant)
+        virtual tmp<DimensionedField<scalar, volMesh> > Ru() const;
+
+        // Access
+
+        //- Intensity Ray on i direction
+        inline const radiativeIntensityRay& RadIntRay(label i) const
+        {
+            return  RadIntRay_[i];
+        }
+
+        //- Intensity Ray on i direction and j band-width
+        inline const volScalarField& RadIntRayiLambdaj
+        (
+            const label i,
+            const label j
+        ) const
+        {
+            return  RadIntRay_[i].Ilambdaj(j);
+        }
+
+        //-Number of angles in theta
+        label Ntheta() const
+        {
+            return Ntheta_;
+        }
+
+        //- Number of angles in phi
+        label Nphi() const
+        {
+            return Nphi_;
+        }
+
+        //- Number of directions
+        label Ni() const
+        {
+            return Ni_;
+        }
+
+        //- Number of wavelengths
+        inline const label& lambdaj() const
+        {
+            return lambdaj_;
+        }
+
+        // Const access to a
+        inline const volScalarField& a() const
+        {
+            return a_;
+        }
+
+        // Const access to aj
+        inline const volScalarField& aj(label i) const
+        {
+            return aj_[i];
+        }
+
+        // Const access to G
+        inline const volScalarField& G() const
+        {
+            return G_;
+        }
+
+        //  Const access to Qr
+        inline const volScalarField& Qr() const
+        {
+            return Qr_;
+        }
+
+        //  Const access to blavkBody
+        virtual const blackBodyEmission& blackBody() const
+        {
+            return blackBody_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C

@@ -0,0 +1,540 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "IFstream.H"
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template <class Type>
+Foam::label Foam::interpolationLookUpTable <Type>::index
+(
+	const List<scalar>& indices,
+	const bool lastDim
+) const
+{
+    label totalindex = 0;
+
+    for(int i = 0;i < dim_.size()-1;i++)
+    {
+        label dim = 1;
+        for(int j = i + 1 ;j <  dim_.size() ; j++)
+		{
+			dim *=(dim_[j]+1);
+		}
+        totalindex += Foam::min
+		(
+			Foam::max(label((indices[i]-min_[i])/delta_[i]),0),
+				dim_[i]
+		)*dim;
+    }
+
+	if(lastDim)
+	{
+
+
+    	label iLastdim = dim_.size() -1;
+    	totalindex  += Foam::min
+		(
+			Foam::max
+			(
+				label((indices[iLastdim]-min_[iLastdim])/delta_[iLastdim]),
+				0
+			),
+			dim_[iLastdim]
+		);
+	}
+
+    return totalindex;
+}
+
+template <class Type>
+Foam::label Foam::interpolationLookUpTable <Type>::index
+(
+    const scalar indice
+) const
+{
+
+    label i = 0;
+    label totalindex =
+    Foam::min
+    (
+        Foam::max
+        (
+            label((indice-min_[i])/delta_[i]),
+            0
+        ),
+        dim_[i]
+    );
+
+    return totalindex;
+}
+
+template<class Type>
+bool Foam::interpolationLookUpTable<Type>::checkRange
+(
+	const scalar lookUpValue,
+	const label interfield
+) const
+{
+	if(lookUpValue >= min_[interfield] &&  lookUpValue <= max_[interfield])
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+}
+
+template<class Type>
+Foam::scalar Foam::interpolationLookUpTable<Type>::interpolate
+(
+	const label lo,
+	const label hi,
+	const scalar lookUpValue, //Xi
+	const label ofield, // Yo , Delta Y
+	const label interfield // Delta X
+) const
+{
+
+
+		if(	List<scalarField>::operator[](interfield).operator[](hi)
+          != List<scalarField>::operator[](interfield).operator[](lo))
+        {
+
+		scalar output
+    	(
+		List<scalarField>::operator[](ofield).operator[](lo)
+     +  (
+        	List<scalarField>::operator[](ofield).operator[](hi)
+       	  - List<scalarField>::operator[](ofield).operator[](lo)
+        )
+      * (
+        	lookUpValue
+          - List<scalarField>::operator[](interfield).operator[](lo)
+        )
+       /(
+        	List<scalarField>::operator[](interfield).operator[](hi)
+          - List<scalarField>::operator[](interfield).operator[](lo)
+        )
+        );
+        	return output;
+    	}
+    	else
+    	{
+    		return List<scalarField>::operator[](ofield).operator[](lo);
+    	}
+
+
+}
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::dimensionTable()
+{
+    min_.setSize(entries_.size());
+    dim_.setSize(entries_.size());
+	delta_.setSize(entries_.size());
+    max_.setSize(entries_.size());
+	entryIndices_.setSize(entries_.size());
+	outputIndices_.setSize(output_.size());
+	label index = 0;
+	label tableDim = 1;
+
+    forAll(entries_,i)
+    {
+        dim_[i] = readLabel(entries_[i].lookup("N"));
+	    max_[i] = readScalar(entries_[i].lookup("max"));
+    	min_[i] = readScalar(entries_[i].lookup("min"));
+		delta_[i] = (max_[i] - min_[i]) / dim_[i];
+        tableDim *= (dim_[i] + 1);
+		fieldIndices_.insert(entries_[i].lookup("name"),index);
+		entryIndices_[i] = index;
+		index ++;
+    }
+
+	forAll(output_,i)
+    {
+		fieldIndices_.insert(output_[i].lookup("name"),index);
+		outputIndices_[i] = index;
+		index++;
+	}
+
+    List<scalarField >& internal = *this;
+
+	internal.setSize(entries_.size()+output_.size());
+
+	interpOutput_.setSize(entries_.size() + output_.size());
+
+	forAll(internal, i)
+	{
+		internal[i].setSize(tableDim);
+	}
+}
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::readTable
+(
+    const word& instance,
+    const fvMesh& mesh
+)
+{
+	IOdictionary control
+    (
+        IOobject
+        (
+            fileName_, //,
+            instance, //i,
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    control.lookup("fields") >> entries_;
+    control.lookup("output") >> output_;
+	control.lookup("values") >> *this;
+
+	dimensionTable();
+
+    check();
+
+    if (this->size() == 0)
+    {
+        FatalErrorIn
+        (
+                "Foam::interpolationLookUpTable<Type>::readTable()"
+        )   << "table is empty" << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable()
+:
+	List<scalarField >(),
+    fileName_("fileNameIsUndefined")
+{}
+
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+	const fileName& fn, const word& instance, const fvMesh& mesh
+)
+:
+	List<scalarField >(),
+    fileName_(fn),
+	dim_(0),
+	min_(0),
+	delta_(0.),
+	max_(0.),
+	entries_(0),
+	output_(0),
+	entryIndices_(0),
+	outputIndices_(0),
+	interpOutput_(0)
+{
+    readTable(instance, mesh);
+}
+
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+     const interpolationLookUpTable& interpTable
+)
+:
+	List<scalarField >(interpTable),
+    fileName_(interpTable.fileName_),
+	entryIndices_(interpTable.entryIndices_),
+	outputIndices_(interpTable.outputIndices_),
+	dim_(interpTable.dim_),
+	min_(interpTable.min_),
+	delta_(interpTable.delta_),
+	max_(interpTable.max_),
+	entries_(0),
+	output_(0),
+	interpOutput_(interpTable.interpOutput_)
+{}
+
+template<class Type>
+Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
+(
+	const dictionary& dict
+)
+:
+	List<scalarField >(),
+    fileName_(fileName(dict.lookup("fileName")).expand()),
+    dim_(0),
+    min_(.0),
+    delta_(.0),
+    max_(.0),
+    entries_(dict.lookup("fields")),
+	output_(dict.lookup("output")),
+	entryIndices_(0),
+	outputIndices_(0),
+    fieldIndices_(0),
+	interpOutput_(0)
+{
+	dimensionTable();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::check() const
+{
+
+//  check order in the first dimension.
+
+    scalar prevValue = List<scalarField >::operator[](0).operator[](0);
+    label dim = 1 ;
+    for(int j = 1 ;j < dim_.size() ; j++)
+	{
+		dim *=(dim_[j]+1);
+	}
+
+    for (label i = 1; i < dim_[0]; i++)
+    {
+        label index = i*dim;
+        const scalar currValue =
+        List<scalarField >::operator[](0).operator[](index);
+        // avoid duplicate values (divide-by-zero error)
+        if (currValue <= prevValue)
+        {
+            FatalErrorIn
+            (
+                    "Foam::interpolationLookUpTable<Type>::checkOrder() const"
+            )   << "out-of-order value: "
+                << currValue << " at index " << index << nl
+                << exit(FatalError);
+        }
+        prevValue = currValue;
+    }
+}
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::write
+(
+    Ostream& os,
+    const fileName& fn,
+    const word& instance,
+    const fvMesh& mesh
+) const
+{
+    IOdictionary control
+    (
+        IOobject
+        (
+            fn, //,
+            instance, //i,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+   control.writeHeader(os) ;
+
+   os.writeKeyword("fields");
+   os << entries_ << token::END_STATEMENT << nl;
+
+   os.writeKeyword("output");
+   os << output_ << token::END_STATEMENT << nl;
+
+    if (this->size() == 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::interpolationTable<Type>::write()"
+        )   << "table is empty" << nl
+            << exit(FatalError);
+    }
+	os.writeKeyword("values");
+    os << *this << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::scalarField&
+Foam::interpolationLookUpTable<Type>::operator[]
+(
+	const label i
+)
+{
+    label ii = i;
+    label n  = this->size();
+
+    if (n <= 1)
+    {
+        FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			 "(const label) const"
+		)   << "table has (" << n << ") columns" << nl
+            << exit(FatalError);
+    }
+    else if (ii < 0)
+    {
+		FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			 "(const label) const"
+		)   << "index (" << ii << ") underflow" << nl
+            << exit(FatalError);
+	}
+
+    else if (ii > n)
+    {
+		FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			"(const label) const"
+		)   << "index (" << ii << ") overflow" << nl
+            << exit(FatalError);
+	}
+
+	return List<scalarField >::operator[](ii);
+}
+
+
+template<class Type>
+const Foam::scalarField&
+Foam::interpolationLookUpTable<Type>::operator[]
+(
+	const label i
+) const
+{
+    label ii = i;
+    label n  = this->size();
+
+    if (n <= 1)
+    {
+        FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			 "(const label) const"
+		)   << "table has (" << n << ") columns" << nl
+            << exit(FatalError);
+    }
+    else if (ii < 0)
+    {
+		FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			 "(const label) const"
+		)   << "index (" << ii << ") underflow" << nl
+            << exit(FatalError);
+	}
+
+    else if (ii > n)
+    {
+		FatalErrorIn
+        (
+			"Foam::interpolationLookUpTable<Type>::operator[]"
+			"(const label) const"
+		)   << "index (" << ii << ") overflow" << nl
+            << exit(FatalError);
+	}
+
+	return List<scalarField >::operator[](ii);
+}
+
+template<class Type>
+bool Foam::interpolationLookUpTable<Type>::found
+(
+    const word& key
+) const
+{
+    for
+    (
+        HashTable<label>::const_iterator iter = fieldIndices_.begin();
+        iter != fieldIndices_.end();
+        ++iter
+    )
+    {
+        if (fieldIndices_.found(key))
+        {
+            return true;
+        };
+    }
+    return false;
+}
+
+
+template<class Type>
+const Foam::scalarList&
+Foam::interpolationLookUpTable<Type>::LookUp
+(
+    const scalar retvals
+)
+{
+    const label lo = index(retvals);
+    findHi(lo, retvals);
+    return interpOutput_;
+}
+
+
+template<class Type>
+void Foam::interpolationLookUpTable<Type>::findHi
+(
+    const label lo,
+    const scalar retvals
+)
+{
+
+    forAll(outputIndices_,j)
+    {
+        scalar tmp = 0;
+        label ofield = outputIndices_[j];
+        scalar baseValue = List<scalarField>::operator[](ofield).operator[](lo);
+
+        forAll(entryIndices_,i)
+        {
+            if (checkRange(retvals,entryIndices_[i]))
+            {
+                label dim = 1;
+
+                label hi = Foam::min(lo + dim,(*this)[0].size()-1);
+
+                tmp +=(interpolate(lo,hi,retvals,ofield,entryIndices_[i])
+                    - baseValue);
+            }
+            interpOutput_[entryIndices_[i]] = retvals;
+        }
+
+        tmp += baseValue;
+        interpOutput_[outputIndices_[j]] = tmp;
+    }
+
+}

src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H

@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::interpolationLookUpTable
+
+Description
+    A list of lists. Interpolates based on the first dimension.
+    The values must be positive and monotonically increasing in each dimension
+
+
+Note
+    - Accessing an empty list results in an error.
+    - Accessing a list with a single element always returns the same value.
+
+SourceFiles
+    interpolationLookUpTable.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationLookUpTable_H
+#define interpolationLookUpTable_H
+
+#include "List.H"
+#include "ListOps.H"
+#include "scalarField.H"
+#include "HashTable.H"
+#include "IOdictionary.H"
+#include "fvCFD.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class interpolationLookUpTable Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolationLookUpTable
+:
+	  public List<scalarField>
+{
+public:
+
+    // Public data types
+
+private:
+
+    // Private data
+
+        //- File name
+        fileName fileName_;
+
+	//- table dimensions
+        List<label>  dim_;
+
+	//- min on each dimension
+        List<scalar>  min_;
+
+	//- deltas on each dimension
+		List<scalar>  delta_;
+
+	//- maximum on each dimension
+        List<scalar>  max_;
+
+        //- entries dictionaries
+        List<dictionary> entries_;
+
+		//- output dictionaries
+		List<dictionary> output_;
+
+		//- input Indeces from the Look Up Table
+		List<label> entryIndices_;
+
+		//- output Indeces from the Look Up Table
+		List<label> outputIndices_;
+
+		//- field names and indices
+		HashTable<label> fieldIndices_;
+
+		//- Output List containing input and interpolation values of outputs
+		List<scalar> interpOutput_;
+
+        // Private Member Functions
+
+		//- Read the table of data from file
+        void readTable(const word& instance, const fvMesh& mesh);
+
+		//- Dimension Table from dictionaries input and output
+		void dimensionTable();
+
+        // Index table finding using scalarList
+		label index(const List<scalar>&, const bool lastDim=true) const;
+
+        // Index table finding using scalar
+        label index(const scalar) const;
+
+		//- check range of lookUpValue
+		bool checkRange(const scalar, const label) const;
+
+		//- Interpolate function return an scalar
+		scalar interpolate
+		(
+			const label lo,
+			const label hi,
+			const scalar lookUpValue,
+			const label ofield,
+			const label interfield
+		) const;
+
+		// Check list is monotonically increasing
+        void check() const;
+
+        // find hi index, interpolate and populate interpOutput_
+        void findHi(const label lo, const scalar retvals);
+
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        interpolationLookUpTable();
+
+        //- Construct given the name of the file containing the table of data
+        interpolationLookUpTable
+        (
+            const fileName& fn,
+            const word& instance,
+            const fvMesh& mesh
+        );
+
+         //- Construct from dictionary
+        interpolationLookUpTable(const dictionary& dict);
+
+
+		//- Construct copy
+        interpolationLookUpTable(const interpolationLookUpTable& interpTable);
+
+    // Member Functions
+
+        bool found(const word& key) const;
+
+        //- Return the output list given a single input scalar
+        const List<scalar>& LookUp(const scalar);
+
+		//- Write Look Up Table to filename.
+        void write
+        (
+            Ostream& os,
+            const fileName& fn,
+            const word& instance,
+            const fvMesh& mesh
+        ) const;
+
+
+      // Acces
+
+		inline label findFieldIndex(const word& field) const
+		{
+			return fieldIndices_[field];
+		}
+
+		 inline const List<dictionary>& output() const
+        {
+            return output_;
+        }
+
+		inline const List<dictionary>& entries() const
+        {
+            return entries_;
+        }
+
+        inline const List<scalar>& min() const
+        {
+            return min_;
+        }
+
+        inline const List<label>& dim() const
+        {
+            return dim_;
+        }
+
+		inline const List<scalar>& delta() const
+		{
+	    	return delta_;
+		}
+
+        inline const List<scalar>& max() const
+        {
+            return max_;
+        }
+
+        inline const word tableName() const
+        {
+            return fileName_.name();
+        }
+     // Member Operators
+
+        //- Return an element of constant List<scalar, Type>
+        const scalarField& operator[](const label) const;
+
+		//- Return an element of List<scalar, Type>
+		scalarField & operator[](const label);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#ifdef NoRepository
+#   include "interpolationLookUpTable.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -0,0 +1,250 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "radiativeIntensityRay.H"
+#include "fvm.H"
+#include "fvDOM.H"
+
+#include "absorptionEmissionModel.H"
+#include "scatterModel.H"
+#include "mathematicalConstants.H"
+#include "radiationConstants.H"
+#include "radiationModel.H"
+#include "Vector2D.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+Foam::label Foam::radiation::radiativeIntensityRay::rayId = 0;
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
+(
+    scalar& phii,
+    scalar& thetai,
+    scalar& deltaPhi,
+    scalar& deltaTheta,
+    label& lambdaj,
+    const fvMesh& mesh,
+    const absorptionEmissionModel& absEmmModel,
+    const blackBodyEmission& blackBody
+)
+:
+    absEmmModel_(absEmmModel),
+    mesh_(mesh),
+    blackBody_(blackBody),
+    I_
+    (
+        IOobject
+        (
+            "I" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qri_
+    (
+        IOobject
+        (
+            "Qr" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+    )
+{
+    init(phii,thetai,deltaPhi,deltaTheta,lambdaj);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::radiation::radiativeIntensityRay::init
+(
+    scalar phii, scalar thetai, scalar deltaPhi, scalar deltaTheta,
+    scalar lambdaj
+)
+{
+    scalar sinTheta = Foam::sin(thetai);
+    scalar cosTheta = Foam::cos(thetai);
+    scalar sinPhi = Foam::sin(phii);
+    scalar cosPhi = Foam::cos(phii);
+    vector s = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
+    Si_ = (s);
+    thetai_ = thetai;
+    phii_ = phii;
+    omegai_ = 2. * Foam::sin(thetai)*Foam::sin(deltaTheta/2.)*deltaPhi;
+    Nlambdaj_ = (lambdaj);
+    Ilambdaj_.setSize(Nlambdaj_);
+    const vector& d = vector
+    (
+        sinPhi * Foam::sin(0.5*deltaPhi) * (deltaTheta - Foam::cos(2.*thetai)
+        * Foam::sin(deltaTheta)) ,
+        cosPhi * Foam::sin(0.5*deltaPhi) * (deltaTheta - Foam::cos(2.*thetai)
+        * Foam::sin(deltaTheta)) ,
+        0.5 * deltaPhi * Foam::sin(2.*thetai) * Foam::sin(deltaTheta)
+    );
+
+    Di_ = (d);
+
+    forAll(Ilambdaj_, i)
+    {
+        IOobject header
+        (
+            "Ilambda_" + name(rayId) + "_"+ name(i),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ
+        );
+        // check if field exists and can be read
+        if (header.headerOk())
+        {
+            Ilambdaj_.set
+            (
+                i,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "Ilambda_" + name(rayId) + "_" + name(i),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::MUST_READ,
+                        IOobject::AUTO_WRITE
+                    ),
+                    mesh_
+                )
+            );
+        }
+        else
+        {
+            volScalarField Idefault
+            (
+                IOobject
+                (
+                    "Idefault",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_
+            );
+
+            Ilambdaj_.set
+            (
+                i,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "Ilambda_" + name(rayId) + "_"+ name(i),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::AUTO_WRITE
+                    ),
+                    Idefault
+                )
+            );
+        }
+    }
+    rayId++;
+}
+
+Foam::scalar Foam::radiation::radiativeIntensityRay::correct
+(
+    fvDOM* DomPtr
+)
+{
+    Qri_ =  dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+
+    scalar maxResidual = 0.0;
+
+    for(label i = 0; i < Nlambdaj_; i++)
+    {
+        volScalarField k = DomPtr->aj(i);
+
+        volScalarField E = absEmmModel_.ECont(i)/Foam::mathematicalConstant::pi;
+
+        surfaceScalarField Ji = Di_ & mesh_.Sf();
+
+        volScalarField Ib = blackBody_.bj(i)/Foam::mathematicalConstant::pi;
+
+        fvScalarMatrix IiEq
+        (
+            fvm::div(Ji,Ilambdaj_[i],"div(Ji,Ii_h)")
+          + fvm::Sp(k*omegai_, Ilambdaj_[i])
+            == k*omegai_*Ib + E
+        );
+
+        IiEq.relax();
+
+        scalar eqnResidual = solve
+        (
+            IiEq,
+            mesh_.solver("Ii")
+        ).initialResidual();
+
+        maxResidual = max(eqnResidual, maxResidual);
+
+    }
+    return maxResidual;
+}
+
+void Foam::radiation::radiativeIntensityRay::addIntensity()
+{
+    I_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+
+    for(label i = 0; i < Nlambdaj_; i++)
+    {
+         I_ += absEmmModel_.addRadInt(i, Ilambdaj_[i]);
+    }
+}
+
+void Foam::radiation::radiativeIntensityRay::add
+(
+    const scalarField & qr,
+    label patchI
+) const
+{
+      Qri_.boundaryField()[patchI] += qr;
+}
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::radiativeIntensityRay
+
+Description
+
+SourceFiles
+    radiativeIntensityRay.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiativeIntensityRay_H
+#define radiativeIntensityRay_H
+
+#include "absorptionEmissionModel.H"
+#include "blackBodyEmission.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+
+/*---------------------------------------------------------------------------*\
+                           Class radiativeIntensityRay Declaration
+\*---------------------------------------------------------------------------*/
+class fvDOM;
+
+
+class radiativeIntensityRay
+{
+
+    // Private data
+
+        //- absorptionEmissionModel
+        const absorptionEmissionModel& absEmmModel_;
+
+        //- Temperature
+        const fvMesh& mesh_;
+
+        // black Body
+        const blackBodyEmission&  blackBody_;
+
+        //- Total Radiative Intensity on i direction  / [W/m2]
+        volScalarField I_;
+
+        // - Total  radiative heat flux on boundary on i direction
+        mutable volScalarField Qri_;
+
+        //- Direction i
+        vector Si_;
+
+        //- theta angle of direction i
+        scalar thetai_;
+
+        //- phi angle of direction i
+        scalar phii_;
+
+        //- solid angle
+        scalar omegai_;
+
+        //- Number of bands on i direction
+        label Nlambdaj_;
+
+        //- average vector inside the solid angle
+        vector Di_;
+
+        //- List of Pointers to Radiative Intensity wave-length
+        PtrList<volScalarField> Ilambdaj_;
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        radiativeIntensityRay(const radiativeIntensityRay&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const radiativeIntensityRay&);
+
+        //- Integrate Intensity on this direction
+//        void IntegrateRadiativeIntensity(const label i);
+public:
+
+    // Static data members
+        static label rayId;
+
+    // Constructors
+
+        //- Null constructor
+        radiativeIntensityRay
+        (
+            scalar& phii,
+            scalar& thetai,
+            scalar& deltaPhi,
+            scalar& deltaTheta,
+            label& lambdaj,
+            const fvMesh& mesh_,
+            const absorptionEmissionModel& absEmmModel_,
+            const blackBodyEmission& blackBody
+        );
+
+
+    // Destructor
+
+        ~radiativeIntensityRay();
+
+
+    // Member functions
+
+        // Edit
+
+            // Update radiative intensity on i direction
+            scalar correct(fvDOM*);
+
+            // init the ray on i direction
+            void init
+            (
+                scalar phii, scalar thetai, scalar deltaPhi,scalar
+                deltaTheta, scalar lambda
+            );
+
+            // add radiative heat flux on walls from the boundary patch
+            void add(const scalarField&, label) const;
+
+            // add Radiative intensities from all the bands
+            void addIntensity();
+
+        // Access
+
+            //- Return Intensity on i direction
+            inline const volScalarField& I() const
+            {
+                return I_;
+            }
+
+            //- Return heat flux on boundary on i direction
+            inline const volScalarField& Qri() const
+            {
+                return Qri_;
+            }
+
+            inline const vector Si() const
+            {
+                return Si_;
+            }
+
+            inline const vector Di() const
+            {
+                return Di_;
+            }
+
+            scalar lambdaj() const
+            {
+                return  Nlambdaj_;
+            }
+
+            scalar phii() const
+            {
+                return  phii_;
+            }
+
+            scalar thetai() const
+            {
+                return  thetai_;
+            }
+
+            scalar omegai() const
+            {
+                return  omegai_;
+            }
+
+
+            const volScalarField& Ilambdaj(label i) const
+            {
+                return  Ilambdaj_[i];
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C

@@ -52,8 +52,8 @@ autoPtr<radiationModel> radiationModel::New
             IOobject
             (
                 "radiationProperties",
-                T.time().constant(),
+                T.mesh().time().constant(),
-                T.db(),
+                T.mesh().db(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
@@ -86,7 +86,7 @@ autoPtr<radiationModel> radiationModel::New
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace radiation
+} // End radiation
 } // End namespace Foam
 
 // ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C

@@ -54,8 +54,8 @@ Foam::radiation::radiationModel::radiationModel
         IOobject
         (
             "radiationProperties",
-            T.time().constant(),
+            T.mesh().time().constant(),
-            T.db(),
+            T.mesh().db(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
@@ -64,6 +64,7 @@ Foam::radiation::radiationModel::radiationModel
     mesh_(T.mesh()),
     radiation_(lookup("radiation")),
     radiationModelCoeffs_(subDict(type + "Coeffs")),
+    nFlowIterPerRadIter_(readLabel(lookup("nFlowIterPerRadIter"))),
     absorptionEmission_(absorptionEmissionModel::New(*this, mesh_)),
     scatter_(scatterModel::New(*this, mesh_))
 {}

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H

@@ -49,6 +49,7 @@ SourceFiles
 #include "volFields.H"
 #include "basicThermo.H"
 #include "fvMatrices.H"
+#include "blackBodyEmission.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,9 +83,12 @@ protected:
         //- Model specific dictionary input parameters
         Switch radiation_;
 
+
         //- Radiation model dictionary
         dictionary radiationModelCoeffs_;
 
+        //- Number of iteration in the Flow solver per Radiative solver
+        label nFlowIterPerRadIter_;
 
         // References to the radiation sub-models
 
@@ -95,6 +99,7 @@ protected:
             autoPtr<scatterModel> scatter_;
 
 
+
 private:
 
     // Private Member Functions
@@ -106,6 +111,7 @@ private:
         void operator=(const radiationModel&);
 
 
+
 public:
 
     //- Runtime type information
@@ -174,6 +180,8 @@ public:
             (
                 basicThermo& thermo
             ) const;
+
+
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C

@@ -59,14 +59,14 @@ Foam::radiation::absorptionEmissionModel::~absorptionEmissionModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::a() const
+Foam::radiation::absorptionEmissionModel::a(label i) const
 {
-    return aDisp() + aCont();
+    return aDisp(i) + aCont(i);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aCont() const
+Foam::radiation::absorptionEmissionModel::aCont(label i) const
 {
     return tmp<volScalarField>
     (
@@ -89,7 +89,7 @@ Foam::radiation::absorptionEmissionModel::aCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aDisp() const
+Foam::radiation::absorptionEmissionModel::aDisp(label i) const
 {
     return tmp<volScalarField>
     (
@@ -112,14 +112,14 @@ Foam::radiation::absorptionEmissionModel::aDisp() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::e() const
+Foam::radiation::absorptionEmissionModel::e(label i) const
 {
-    return eDisp() + eCont();
+    return eDisp(i) + eCont(i);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eCont() const
+Foam::radiation::absorptionEmissionModel::eCont(label i) const
 {
     return tmp<volScalarField>
     (
@@ -142,7 +142,7 @@ Foam::radiation::absorptionEmissionModel::eCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eDisp() const
+Foam::radiation::absorptionEmissionModel::eDisp(label i) const
 {
     return tmp<volScalarField>
     (
@@ -165,14 +165,14 @@ Foam::radiation::absorptionEmissionModel::eDisp() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::E() const
+Foam::radiation::absorptionEmissionModel::E(label i) const
 {
-    return EDisp() + ECont();
+    return EDisp(i) + ECont(i);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::ECont() const
+Foam::radiation::absorptionEmissionModel::ECont(label i) const
 {
     return tmp<volScalarField>
     (
@@ -195,7 +195,7 @@ Foam::radiation::absorptionEmissionModel::ECont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::EDisp() const
+Foam::radiation::absorptionEmissionModel::EDisp(label i) const
 {
     return tmp<volScalarField>
     (
@@ -216,5 +216,41 @@ Foam::radiation::absorptionEmissionModel::EDisp() const
     );
 }
 
+Foam::label
+Foam::radiation::absorptionEmissionModel::nBands() const
+{
+    return pTraits<label>::one;
+}
 
+const Foam::Vector2D<Foam::scalar>&
+Foam::radiation::absorptionEmissionModel::bands(label n) const
+{
+
+    return Vector2D<scalar>::one;
+}
+
+bool Foam::radiation::absorptionEmissionModel::isGrey(void) const
+{
+    return false;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::addRadInt
+(
+    const label i,
+    const volScalarField& Ilambdaj
+) const
+{
+    return Ilambdaj;
+}
+
+void Foam::radiation::absorptionEmissionModel::correct
+(
+    volScalarField& a,
+    PtrList<volScalarField>& aj
+) const
+{
+       a.internalField()  =  this->a();
+       aj[0].internalField() =  a.internalField();
+}
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H

@@ -38,6 +38,7 @@ Description
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 #include "volFields.H"
+#include "Vector2D.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -115,37 +116,69 @@ public:
             // Absorption coefficient
 
                 //- Absorption coefficient (net)
-                virtual tmp<volScalarField> a() const;
+                virtual tmp<volScalarField> a(const label i = 0) const;
 
                 //- Absorption coefficient for continuous phase
-                virtual tmp<volScalarField> aCont() const;
+                virtual tmp<volScalarField> aCont(const label i = 0) const;
 
                 //- Absorption coefficient for dispersed phase
-                virtual tmp<volScalarField> aDisp() const;
+                virtual tmp<volScalarField> aDisp(const label i = 0) const;
 
 
             // Emission coefficient
 
                 //- Emission coefficient (net)
-                virtual tmp<volScalarField> e() const;
+                virtual tmp<volScalarField> e(const label i = 0) const;
 
                 //- Return emission coefficient for continuous phase
-                virtual tmp<volScalarField> eCont() const;
+                virtual tmp<volScalarField> eCont(const label i = 0) const;
 
                 //- Return emission coefficient for dispersed phase
-                virtual tmp<volScalarField> eDisp() const;
+                virtual tmp<volScalarField> eDisp(const label i = 0) const;
 
 
             // Emission contribution
 
                 //- Emission contribution (net)
-                virtual tmp<volScalarField> E() const;
+                virtual tmp<volScalarField> E(const label i = 0) const;
 
                 //- Emission contribution for continuous phase
-                virtual tmp<volScalarField> ECont() const;
+                virtual tmp<volScalarField> ECont(const label i = 0) const;
 
                 //- Emission contribution for dispersed phase
-                virtual tmp<volScalarField> EDisp() const;
+                virtual tmp<volScalarField> EDisp(const label i = 0) const;
+
+                inline const fvMesh& mesh() const
+                {
+                    return mesh_;
+                }
+
+        // Default for grey absorptionEmission model return 1
+        virtual label nBands() const;
+
+        // Default for grey absorptionEmission model return Vector2D::one
+        virtual const Vector2D<scalar>& bands
+        (
+            label n
+        ) const;
+
+        // Is the absorptionEmission grey
+        virtual bool isGrey(void) const;
+
+
+        // Add radiative intensity fir ray i
+        virtual tmp<volScalarField> addRadInt
+        (
+            const label i,
+            const volScalarField& Ilambdaj
+        ) const;
+
+        // Correct absorption coefficients
+        virtual void correct
+        (
+            volScalarField& a_,
+            PtrList<volScalarField>& aj_
+        ) const ;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.C

@@ -70,7 +70,7 @@ Foam::radiation::constantAbsorptionEmission::~constantAbsorptionEmission()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::aCont() const
+Foam::radiation::constantAbsorptionEmission::aCont(label i) const
 {
     tmp<volScalarField> ta
     (
@@ -95,7 +95,7 @@ Foam::radiation::constantAbsorptionEmission::aCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::eCont() const
+Foam::radiation::constantAbsorptionEmission::eCont(label i) const
 {
     tmp<volScalarField> te
     (
@@ -120,7 +120,7 @@ Foam::radiation::constantAbsorptionEmission::eCont() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::ECont() const
+Foam::radiation::constantAbsorptionEmission::ECont(label i) const
 {
     tmp<volScalarField> tE
     (

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.H

@@ -98,19 +98,26 @@ public:
             // Absorption coefficient
 
                 //- Absorption coefficient for continuous phase
-                tmp<volScalarField> aCont() const;
+                tmp<volScalarField> aCont(const label i = 0) const;
 
 
             // Emission coefficient
 
                 //- Emission coefficient for continuous phase
-                tmp<volScalarField> eCont() const;
+                tmp<volScalarField> eCont(const label i = 0) const;
 
 
             // Emission contribution
 
                 //- Emission contribution for continuous phase
-                tmp<volScalarField> ECont() const;
+                tmp<volScalarField> ECont(const label i = 0) const;
+
+    // Member Functions
+
+        inline bool isGrey(void) const
+        {
+            return true;
+        }
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C

@@ -0,0 +1,273 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "greyMeanAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(greyMeanAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            greyMeanAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+    speciesNames_(0),
+    specieIndex_(0),
+    LookUpTable_
+    (
+        fileName(coeffsDict_.lookup("LookUpTableFileName")),
+        "constant",
+        mesh
+    ),
+    thermo_(mesh.db().lookupObject<basicThermo>("thermophysicalProperties")),
+    EhrrCoeff_(readScalar(coeffsDict_.lookup("EhrrCoeff"))),
+    Yj_(0)
+{
+    Yj_.setSize(nSpecies_);
+    label nFunc = 0;
+    const dictionary& functionDicts = dict.subDict(typeName+"Coeffs");
+
+    forAllConstIter(dictionary, functionDicts, iter)
+    {
+        // safety:
+        if (!iter().isDict())
+        {
+            continue;
+        }
+        const word& key = iter().keyword();
+        speciesNames_.insert(key, nFunc);
+        const dictionary& dict = iter().dict();
+        coeffs_[nFunc].init(dict);
+        nFunc++;
+    }
+
+//  Check that all the species on the dictionary are present in the LookupTable
+//  and save the corresponding indexes of the LookupTable
+
+    label j = 0;
+    forAllConstIter(HashTable<label>, speciesNames_, iter)
+    {
+        if(mesh.db().foundObject<volScalarField>("ft"))
+        {
+
+            if(LookUpTable_.found(iter.key()))
+            {
+                label index = LookUpTable_.findFieldIndex(iter.key());
+
+                Info << "specie: " << iter.key() << " found on LookUpTable"
+                << " with index: " << index << endl;
+
+                specieIndex_[iter()] = index;
+            }
+            else if(mesh.db().foundObject<volScalarField>(iter.key()))
+            {
+                volScalarField& Y = const_cast<volScalarField&>
+                        (mesh.db().lookupObject<volScalarField>(iter.key()));
+                Yj_.set
+                (
+                    j,
+                    &Y
+                );
+                specieIndex_[iter()] = 0.;
+                j++;
+                Info << "specie : " << iter.key() << " is being solved"
+                     << endl;
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "Foam::radiation::greyMeanAbsorptionEmission(const"
+                    "dictionary& dict, const fvMesh& mesh)"
+                )   << "specie : " << iter.key()
+                    << " is neither in Look Up Table : "
+                    << LookUpTable_.tableName()
+                    << " nor is being solved" <<nl
+                    << exit(FatalError);
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "Foam::radiation::greyMeanAbsorptionEmission(const"
+                "dictionary& dict, const fvMesh& mesh)"
+            )   << "specie ft is not present " <<nl
+                << exit(FatalError);
+
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::greyMeanAbsorptionEmission::~greyMeanAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::aCont(label iband) const
+{
+    const volScalarField& T = thermo_.T();
+    const volScalarField& p = thermo_.p();
+    const volScalarField& ft =
+            this->mesh().db().lookupObject<volScalarField>("ft");
+
+    label nSpecies = speciesNames_.size();
+
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "a",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a",dimless/dimLength, 0.0)
+        )
+    );
+
+    scalarField& a = ta().internalField();
+
+    forAll(a, i)
+    {
+        //
+        const List<scalar>& species = LookUpTable_.LookUp(ft[i]);
+
+        for(label n=0; n<nSpecies; n++)
+        {
+            label l = 0;
+            scalar Yipi = 0;
+            if(specieIndex_[n] != 0)
+            {
+                //moles x pressure [atm]
+                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+            }
+            else
+            {
+                Yipi = Yj_[l][i]; // mass fraction
+                l++;
+            }
+
+            const absorptionCoeffs::coeffArray& b = coeffs_[n].coeffs(T[i]);
+
+            scalar Ti = T[i];
+            if (coeffs_[n].invTemp())        //negative Temp exponents
+            {
+                Ti = 1./T[i];
+            }
+            a[i]+= Yipi*
+            (((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti + b[0]);
+        }
+    }
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::eCont(label iband) const
+{
+    tmp<volScalarField> e
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "e",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("e",dimless/dimLength, 0.0)
+        )
+    );
+
+    return e;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::greyMeanAbsorptionEmission::ECont(label iband) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "E",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+
+    if (mesh().db().foundObject<volScalarField>("hrr"))
+    {
+        const volScalarField& hrr =
+                        mesh().db().lookupObject<volScalarField>("hrr");
+        E().internalField() =  EhrrCoeff_ *  hrr.internalField();
+    }
+
+    return E;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::greyMeanAbsorptionEmission
+
+Description
+    greyMeanAbsorptionEmission radiation absorption and emission coefficients
+    for continuous  phase
+
+    The coefficients for the species in the LookUpTable have to be specified
+    for use in moles x P [atm].i.e. (k[i] = species[i]*p*9.869231e-6).
+
+    The coefficients for CO and soot or any other added are multiplied by the
+    respective mass fraction being solved.
+
+    All the species in the dictionary need either to be in the Look up Table or
+    being solved. Conversely, all the species solved do not need to be included
+    in the calculation of the absorption coefficient.
+
+    The names of the species in the absorption dictionary must match exactly the
+    name in the Look Up table or the name of the field being solved.
+
+    The look Up table ("SpeciesTable") file should be in constant
+
+    i.e. dictionary
+
+    LookUpTableFileName     "SpeciesTable";
+
+    EhrrCoeff                0.0;
+
+     CO2
+    {
+        Tcommon         300.;   //Common Temp
+        invTemp         true;   //Is the polynomio using inverse temperature.
+        Tlow            300.;   //Low Temp
+        Thigh           2500.;  //Hight Temp
+
+        loTcoeffs       //coefss for T < Tcommon
+        (
+            0           //  a0            +
+            0           //  a1*T          +
+            0           //  a2*T^(+/-)2   +
+            0           //  a3*T^(+/-)3   +
+            0           //  a4*T^(+/-)4   +
+            0           //  a5*T^(+/-)5   +
+        );
+        hiTcoeffs       //coefss for T > Tcommon
+        (
+            18.741
+            -121.31e3
+            273.5e6
+            -194.05e9
+            56.31e12
+            -5.8169e15
+        );
+
+    }
+
+SourceFiles
+    greyMeanAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef greyMeanAbsorptionEmission_H
+#define greyMeanAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "basicThermo.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class constantAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class greyMeanAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+
+public:
+
+
+        // - Maximum number of species considered for absorptivity
+        static const int nSpecies_ = 5;
+
+        //  Absorption Coefficients
+        absorptionCoeffs coeffs_[nSpecies_];
+
+private:
+        //- Absorption model dictionary
+        dictionary coeffsDict_;
+
+        // Hash table with species names
+        HashTable<label> speciesNames_;
+
+        // Indexes of species in the LookUpTable
+        FixedList<label,nSpecies_> specieIndex_;
+
+
+        // Look Up table of species related to ft
+        mutable interpolationLookUpTable<scalar> LookUpTable_;
+
+        // Thermo
+        const basicThermo& thermo_;
+
+        //- Emission constant coefficient
+        const scalar EhrrCoeff_;
+
+        //- Pointer list of species in the registry involved in the absorption
+        UPtrList<volScalarField> Yj_;
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("greyMeanAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        greyMeanAbsorptionEmission
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Destructor
+
+        ~greyMeanAbsorptionEmission();
+
+
+    // Member Operators
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label i = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label i = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label i = 0) const;
+
+    // Member Functions
+
+        inline bool isGrey(void) const
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C

@@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "wideBandAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(wideBandAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            wideBandAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+    speciesNames_(0),
+    specieIndex_(0),
+    LookUpTable_
+    (
+        fileName(coeffsDict_.lookup("LookUpTableFileName")),
+        "constant",
+        mesh
+    ),
+    thermo_(mesh.db().lookupObject<basicThermo>("thermophysicalProperties")),
+    Yj_(0),
+    totalWaveLength_(0)
+{
+    Yj_.setSize(nSpecies_);
+    label nBand = 0;
+    const dictionary& functionDicts = dict.subDict(typeName +"Coeffs");
+    forAllConstIter(dictionary, functionDicts, iter)
+    {
+        // safety:
+        if (!iter().isDict())
+        {
+            continue;
+        }
+
+        const dictionary& dict = iter().dict();
+        dict.lookup("bandLimits") >> iBands_[nBand];
+        dict.lookup("EhrrCoeff") >> iEhrrCoeffs_[nBand];
+        totalWaveLength_ += (iBands_[nBand][1] - iBands_[nBand][0]);
+
+        label nSpec = 0;
+        const dictionary& SpecDicts = dict.subDict("species");
+        forAllConstIter(dictionary, SpecDicts, iter)
+        {
+            const word& key = iter().keyword();
+            if (nBand == 0)
+            {
+                speciesNames_.insert(key, nSpec);
+            }
+            else
+            {
+                if (!speciesNames_.found(key))
+                {
+                    FatalErrorIn
+                    (
+                        "Foam::radiation::wideBandAbsorptionEmission(const"
+                        "dictionary& dict, const fvMesh& mesh)"
+                    )   << "specie : " << key << "is not in all the bands"
+                        << nl
+                        << exit(FatalError);
+                }
+            }
+            coeffs_[nSpec][nBand].init(SpecDicts.subDict(key));
+            nSpec++;
+        }
+        nBand++;
+    }
+    nBands_ = nBand;
+
+/*
+Check that all the species on the dictionary are present in the
+LookupTable  and save the corresponding indexes of the LookupTable
+*/
+
+    label j = 0;
+    forAllConstIter(HashTable<label>, speciesNames_, iter)
+    {
+
+        if(LookUpTable_.found(iter.key()))
+        {
+            label index = LookUpTable_.findFieldIndex(iter.key());
+            Info << "specie: " << iter.key() << " found in LookUpTable"
+            << " with index: " << index << endl;
+            specieIndex_[iter()] = index;
+        }
+        else if
+        (mesh.db().foundObject<volScalarField>(iter.key()))
+        {
+            volScalarField& Y = const_cast<volScalarField&>
+                    (mesh.db().lookupObject<volScalarField>(iter.key()));
+
+            Yj_.set(j, &Y);
+
+            specieIndex_[iter()] = 0.;
+            j++;
+            Info << "specie : " << iter.key() << " is being solved" << endl;
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "radiation::wideBandAbsorptionEmission(const"
+                "dictionary& dict, const fvMesh& mesh)"
+             )  << "specie : " << iter.key()
+                << " is neither in Look Up Table : "
+                << LookUpTable_.tableName() <<" nor is being slved"
+                << exit(FatalError);
+        }
+    }
+}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::wideBandAbsorptionEmission::~wideBandAbsorptionEmission()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::aCont(label iband) const
+{
+    const volScalarField& T = thermo_.T();
+    const volScalarField& p = thermo_.p();
+    const volScalarField& ft =
+            this->mesh().db().lookupObject<volScalarField>("ft");
+
+    label nSpecies = speciesNames_.size();
+
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "a",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a",dimless/dimLength, 0.0)
+        )
+    );
+
+    scalarField& a = ta().internalField();
+
+    forAll(a, i)
+    {
+
+        const List<scalar>& species = LookUpTable_.LookUp(ft[i]);
+
+        for(label n=0; n<nSpecies; n++)
+        {
+            label l = 0;
+            scalar Yipi = 0;
+            if(specieIndex_[n] != 0) //moles x pressure [atm]
+            {
+                Yipi = species[specieIndex_[n]]*p[i]*9.869231e-6;
+            }
+            else
+            {
+                Yipi = Yj_[l][i]; // mass fraction from species being solved
+                l++;
+            }
+
+            scalar Ti = T[i];
+
+            const absorptionCoeffs::coeffArray& b =
+                                            coeffs_[n][iband].coeffs(T[i]);
+
+            if (coeffs_[n][iband].invTemp())
+            {
+                Ti = 1./T[i];
+            }
+
+            a[i]+=Yipi*
+            (((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti + b[0]);
+        }
+    }
+
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::eCont(label iband) const
+{
+    tmp<volScalarField> e
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "e",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("e",dimless/dimLength, 0.0)
+        )
+    );
+
+    return e;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::ECont(label iband) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "E",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+
+    if (mesh().db().foundObject<volScalarField>("hrr"))
+    {
+        const volScalarField& hrr =
+                        mesh().db().lookupObject<volScalarField>("hrr");
+        E().internalField() =  iEhrrCoeffs_[iband] *  hrr.internalField() *
+                    (iBands_[iband][1] - iBands_[iband][0])/totalWaveLength_;
+    }
+
+
+    return E;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wideBandAbsorptionEmission::addRadInt
+(
+    const label i,
+    const volScalarField& Ilambdaj
+) const
+{
+     return Ilambdaj*(iBands_[i][1] - iBands_[i][0])/totalWaveLength_;
+}
+
+
+void Foam::radiation::wideBandAbsorptionEmission::correct
+(
+    volScalarField& a,
+    PtrList<volScalarField>& aj
+) const
+{
+    a = dimensionedScalar("zero",dimless/dimLength, 0.0);
+
+    for(label j = 0; j < nBands_; j++)
+        {
+            Info << "Calculating... absorption in band : " << j <<endl;
+            aj[j].internalField() = this->a(j);
+            Info << "Calculated absorption in band : " << j <<endl;
+            a.internalField() += aj[j].internalField() *
+                            (iBands_[j][1] - iBands_[j][0])/totalWaveLength_;
+        }
+
+}
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H

@@ -0,0 +1,264 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::radiation::greyMeanAbsorptionEmission
+
+Description
+
+    wideBandAbsorptionEmission radiation absorption and emission coefficients
+    for continuous phase.
+
+    All the bands should have the same number of species and have to be entered
+    in the same order.
+
+    There is no check of continuity of the bands. They should not ovelap or
+    have gaps.
+
+    The black body emission power table(constant/blackBodyEmissivePower) range
+    of lambda * T = [1000; 10000] x 10E-6 (90% of the total emission).
+
+    The emission constant proportionality is specified per band (EhrrCoeff).
+
+    The coefficients for the species in the LookUpTable have to be specified
+    for use in moles x P [atm].i.e. (k[i] = species[i]*p*9.869231e-6).
+
+    The coefficients for CO and soot or any other added are multiplied by the
+    respective mass fraction being solved.
+
+    The look Up table file should be in the constant directory.
+
+    band dictionary:
+
+    band0
+    {
+        bandLimits (1.0e-6 2.63e-6);
+        EhrrCoeff       0.0;
+        species
+        {
+            CH4
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+            CO2
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+
+            H2O
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+
+            Ysoot
+            {
+                Tcommon         300.;
+                Tlow            300.;
+                Thigh           2500.;
+                invTemp         false;
+                loTcoeffs (0 0 0 0 0 0) ;
+                hiTcoeffs (.1 0 0 0 0 0);
+            }
+        }
+    }
+
+
+
+
+
+SourceFiles
+    wideBandAbsorptionEmission.C
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef wideBandAbsorptionEmission_H
+#define wideBandAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "basicThermo.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class wideBandAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class wideBandAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+
+public:
+
+
+        // - Maximum number of species considered for absorptivity
+        static const int nSpecies_ = 5;
+
+        // - Maximum number of bands
+        static const int maxBands_ = 10;
+
+        //  Absorption Coefficients
+        FixedList< FixedList<absorptionCoeffs, nSpecies_> , maxBands_> coeffs_;
+
+private:
+        //- Absorption model dictionary
+        dictionary coeffsDict_;
+
+        //- Hash table with species names
+        HashTable<label> speciesNames_;
+
+        //- Indexes of species in the LookUpTable
+        FixedList<label,nSpecies_> specieIndex_;
+
+        //- Bands
+        FixedList<Vector2D<scalar>, maxBands_ > iBands_;
+
+        //- Proportion of the heat released rate emitted
+        FixedList<scalar, maxBands_ > iEhrrCoeffs_;
+
+        //- Look Up table of species related to ft
+        mutable interpolationLookUpTable<scalar> LookUpTable_;
+
+        //- Thermo
+        const basicThermo& thermo_;
+
+        //- Bands
+        label nBands_ ;
+
+        //- Pointer list of species being solved involved in the absorption
+        UPtrList<volScalarField> Yj_;
+
+        // Total wave lenght cover by the bands
+        scalar totalWaveLength_;
+
+public:
+
+
+    //- Runtime type information
+    TypeName("wideBandAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        wideBandAbsorptionEmission
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Destructor
+
+        ~wideBandAbsorptionEmission();
+
+
+    // Member Operators
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label i = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label i = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label i = 0) const;
+
+    // Member Functions
+
+        inline bool isGrey(void) const
+        {
+            return false;
+        }
+
+        // Number of bands
+        label nBands() const
+        {
+            return nBands_;
+        }
+
+        // lower and upper limit of band i
+        inline const Vector2D<scalar>& bands(label i) const
+        {
+            return iBands_[i];
+        }
+
+        //  Add contribution of Ilambdaj to the total Radiative Intensity on
+        //direction i
+        tmp<volScalarField> addRadInt
+        (
+            const label i,
+            const volScalarField& Ilambdaj
+        ) const;
+
+        void correct(volScalarField& a_, PtrList<volScalarField>& aj) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //