diff --git a/src/thermophysicalModels/radiation/Make/files b/src/thermophysicalModels/radiation/Make/files
index cfce31241011a8f5df4eb006fef0c4ba04fd8683..4e7fb1b11cbaab0879b316bb76a2bde95eb4fa0d 100644
--- a/src/thermophysicalModels/radiation/Make/files
+++ b/src/thermophysicalModels/radiation/Make/files
@@ -1,4 +1,3 @@
-
 /* Radiation constants */
 radiationConstants/radiationConstants.C
 
@@ -31,6 +30,7 @@ submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionE
 /* Boundary conditions */
 derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
-derivedFvPatchFields/GreyDiffusiveRadiation/GreyDiffusiveRadiationMixedFvPatchScalarField.C
-derivedFvPatchFields/WideBandDiffusiveRadiation/WideBandDiffusiveRadiationMixedFvPatchScalarField.C
+derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+
 LIB = $(FOAM_LIBBIN)/libradiation
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index 93d64a1cbdb8f8d37e93510016cb3d7172c8857c..7d4e3f748ba93f1e5a3c7783006767b6b4bf5f9b 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -29,7 +29,6 @@ License
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 
-#include "radiationModel.H"
 #include "fvDOM.H"
 #include "radiationConstants.H"
 #include "mathematicalConstants.H"
@@ -47,9 +46,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(p, iF),
     TName_("undefinedT"),
     emissivity_(0.0),
-    myRayIndex_(0),
-    myWaveLengthIndex_(0),
-    myRayIsInit_(-1),
+    rayId_(0),
+    wavelengthId_(0),
     qr_(0)
 {
     refValue() = 0.0;
@@ -71,9 +69,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 {}
 
@@ -89,9 +86,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
     emissivity_(readScalar(dict.lookup("emissivity"))),
-    myRayIndex_(0),
-    myWaveLengthIndex_(0),
-    myRayIsInit_(-1),
+    rayId_(-1),
+    wavelengthId_(-1),
     qr_(0)
 {
     const scalarField& Tp =
@@ -126,9 +122,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 {}
 
@@ -143,9 +138,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 
 {}
@@ -158,9 +152,8 @@ void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    scalarField::autoMap(m);
-
-    qr_.automap(m);
+    mixedFvPatchScalarField::autoMap(m);
+    qr_.autoMap(m);
 }
 
 
@@ -190,32 +183,28 @@ updateCoeffs()
     const scalarField& Tp =
         patch().lookupPatchField<volScalarField, scalar>(TName_);
 
-    const radiationModel& rad =
-        db().lookupObject<radiationModel>("radiationProperties");
-
-    const fvDOM& dom = refCast<const fvDOM>(rad);
+    const fvDOM& dom = db().lookupObject<fvDOM>("radiationProperties");
 
-    const label patchi = patch().index();
+    const label patchI = patch().index();
 
-    if (dom.lambdaj() == 1)
+    if (dom.nLambda() == 1)
     {
-        if (myRayIsInit_ == -1)
+        if (rayId_ == -1)
         {
-            for (label i=0; i < Dom.Ni() ; i++)
+            for (label rayI=0; rayI < dom.nRay(); rayI++)
             {
-                for (label j=0; j < dom.lambdaj() ; j++)
+                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
                 {
                     const volScalarField& radiationField =
-                        dom.RadIntRayiLambdaj(i,j);
+                        dom.IRayWave(rayI, lambdaI);
                     if
                     (
                         &(radiationField.internalField())
                      == &dimensionedInternalField()
                     )
                     {
-                        myRayIndex_ = i;
-                        myWaveLengthIndex_ = j;
-                        myRayIsInit_ = 0.0;
+                        rayId_ = rayI;
+                        wavelengthId_ = lambdaI;
                         break;
                     }
                 }
@@ -234,41 +223,38 @@ updateCoeffs()
             << exit(FatalError);
     }
 
-    vectorField n = patch().Sf()/patch().magSf();
-
     scalarField& Iw = *this;
+    vectorField n = patch().Sf()/patch().magSf();
 
-    qr_ =  Iw*(-n & dom.RadIntRay(myRayIndex_).Di());
+    radiativeIntensityRay& ray =
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
 
-    dom.RadIntRay(myRayIndex_).add(qr_,patchi);
+    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
 
     forAll(Iw, faceI)
     {
         scalar Ir = 0.0;
 
-        for (label i=0; i < dom.Ni(); i++)
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
         {
-            const vector& si = dom.RadIntRay(i).Si();
-
-            const scalarField& Iface = dom.RadIntRay(i).Ilambdaj
-            (
-                myWaveLengthIndex_
-            ).boundaryField()[patch().index()];
+            const vector& d = dom.IRay(rayI).d();
 
-            scalar InOut = -n[faceI] & si;
+            const scalarField& Iface =
+                dom.IRay(rayI).IWave
+                (
+                    wavelengthId_
+                ).boundaryField()[patchI];
 
-            if (InOut < 0.) // qin into the wall
+            if ((-n[faceI] & d) < 0.0) // qin into the wall
             {
-                const vector& di = dom.RadIntRay(i).Di();
-                Ir += Iface[faceI]*mag(n[faceI] & di);
+                const vector& dAve = dom.IRay(rayI).dAve();
+                Ir += Iface[faceI]*mag(n[faceI] & dAve);
             }
         }
 
-        const vector& mySi = dom.RadIntRay(myRayIndex_).Si();
-
-        scalar InOut = -n[faceI] & mySi;
+        const vector& d = dom.IRay(rayId_).d();
 
-        if (InOut > 0.) //direction out of the wall
+        if ((-n[faceI] & d) > 0.) //direction out of the wall
         {
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
@@ -280,7 +266,7 @@ updateCoeffs()
                /mathematicalConstant::pi;
 
         }
-        else if (InOut < 0.) //direction into the wall
+        else //direction into the wall
         {
             valueFraction()[faceI] = 0.0;
             refGrad()[faceI] = 0.0;
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
index 0c22d52caa4480d1bc78f7329ed7a9736204b9c3..d1a98d67fb349b137d7522113ead5eaca258c6d5 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
@@ -60,14 +60,11 @@ class greyDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
-        //- Direction index
-        label myRayIndex_;
+        //- Ray index
+        label rayId_;
 
-        //- Direction index
-        label myWaveLengthIndex_;
-
-        //- Initialise ray flag
-        label myRayIsInit_;
+        //- Wavelength index
+        label wavelengthId_;
 
         //- Radiative heat flux on walls
         scalarField qr_;
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
index a760ada47c8e988ce5630cc9048ac98a2fe92063..7937da15352120c4ca08e8f4f89265492903db04 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -24,12 +24,11 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "WideBandDiffusiveRadiationMixedFvPatchScalarField.H"
+#include "wideBandDiffusiveRadiationMixedFvPatchScalarField.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 
-#include "radiationModel.H"
 #include "fvDOM.H"
 #include "wideBandAbsorptionEmission.H"
 #include "radiationConstants.H"
@@ -45,11 +44,10 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    TName_("undefined"),
+    TName_("undefinedT"),
     emissivity_(0.0),
-    myRayIndex_(0),
-    myWaveLengthIndex_(0),
-    myRayIsInit_(-1),
+    rayId_(0),
+    wavelengthId_(0),
     qr_(p.size(), 0.0)
 {
     refValue() = 0.0;
@@ -71,9 +69,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 {}
 
@@ -89,9 +86,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
     emissivity_(readScalar(dict.lookup("emissivity"))),
-    myRayIndex_(0),
-    myWaveLengthIndex_(0),
-    myRayIsInit_(-1),
+    rayId_(0),
+    wavelengthId_(0),
     qr_(p.size(), 0.0)
 {
     const scalarField& Tp =
@@ -125,9 +121,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 {}
 
@@ -142,9 +137,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
     emissivity_(ptf.emissivity_),
-    myRayIndex_(ptf.myRayIndex_),
-    myWaveLengthIndex_(ptf.myWaveLengthIndex_),
-    myRayIsInit_(ptf.myRayIsInit_),
+    rayId_(ptf.rayId_),
+    wavelengthId_(ptf.wavelengthId_),
     qr_(ptf.qr_)
 {}
 
@@ -156,9 +150,8 @@ void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    scalarField::autoMap(m);
-
-    qr_.automap(m);
+    mixedFvPatchScalarField::autoMap(m);
+    qr_.autoMap(m);
 }
 
 
@@ -185,32 +178,28 @@ updateCoeffs()
         return;
     }
 
-    const radiationModel& rad =
-        db().lookupObject<radiationModel>("radiationProperties");
-
-    const fvDOM& dom(refCast<const fvDOM>(rad));
+    const fvDOM& dom = db().lookupObject<fvDOM>("radiationProperties");
 
-    const label patchi = patch().index();
+    const label patchI = patch().index();
 
-    if (dom.lambdaj() > 1)
+    if (dom.nLambda() > 1)
     {
-        if (myRayIsInit_ == -1)
+        if (rayId_ == -1)
         {
-            for (label i=0; i < dom.Ni() ; i++)
+            for (label rayI=0; rayI < dom.nRay(); rayI++)
             {
-                for (label j=0; j < dom.lambdaj() ; j++)
+                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
                 {
                     const volScalarField& radiationField =
-                        dom.RadIntRayiLambdaj(i,j);
+                        dom.IRayWave(rayI, lambdaI);
                     if
                     (
                         &(radiationField.internalField())
                       ==&dimensionedInternalField()
                     )
                     {
-                        myRayIndex_ = i;
-                        myWaveLengthIndex_ = j;
-                        myRayIsInit_ = 0.;
+                        rayId_ = rayI;
+                        wavelengthId_ = lambdaI;
                         break;
                     }
                 }
@@ -225,58 +214,55 @@ updateCoeffs()
             "wideBandDiffusiveRadiationMixedFvPatchScalarField::"
             "updateCoeffs"
         )   << " a Non-grey boundary condition is used with a grey"
-            << "absorption model"
-            << exit(FatalError);
+            << "absorption model" << nl << exit(FatalError);
     }
 
-    vectorField n = patch().Sf()/patch().magSf();
-
     scalarField& Iw = *this;
+    vectorField n = patch().Sf()/patch().magSf();
 
-    qr_ =  Iw*(-n & dom.RadIntRay(myRayIndex_).Di());
+    radiativeIntensityRay& ray =
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
 
-    dom.RadIntRay(myRayIndex_).add(qr_,patchi);
+    ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
 
     const scalarField Eb =
-        dom.blackBody().bj(myWaveLengthIndex_).boundaryField()[patchi];
+        dom.blackBody().bj(wavelengthId_).boundaryField()[patchI];
 
     forAll(Iw, faceI)
     {
         scalar Ir = 0.0;
-        for(label i=0; i < dom.Ni(); i++)
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
         {
-            const vector& si = dom.RadIntRay(i).Si();
+            const vector& d = dom.IRay(rayI).d();
 
-            const scalarField& Iface = dom.RadIntRay(i).Ilambdaj
-            (
-                myWaveLengthIndex_
-            ).boundaryField()[patch().index()];
+            const scalarField& Iface =
+                dom.IRay(rayI).IWave
+                (
+                    wavelengthId_
+                ).boundaryField()[patchI];
 
-            scalar InOut = -n[faceI] & si;
 
-            if (InOut < 0.) // qin into the wall
+            if ((-n[faceI] & d) < 0.0) // qin into the wall
             {
-                const vector& di = dom.RadIntRay(i).Di();
-                Ir = Ir + Iface[faceI]*mag(n[faceI] & di);
+                const vector& dAve = dom.IRay(rayI).dAve();
+                Ir = Ir + Iface[faceI]*mag(n[faceI] & dAve);
             }
         }
 
-        const vector& mySi = dom.RadIntRay(myRayIndex_).Si();
-
-        scalar InOut = -n[faceI] & mySi;
+        const vector& d = dom.IRay(rayId_).d();
 
-        if (InOut > 0.) //direction out of the wall
+        if ((-n[faceI] & d) > 0.0) //direction out of the wall
         {
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
             refValue()[faceI] =
                 (
                     Ir*(1.0 - emissivity_)
-                  + emissivity_* Eb[faceI]
+                  + emissivity_*Eb[faceI]
                 )
                /mathematicalConstant::pi;
         }
-        else if (InOut < 0.) //direction into the wall
+        else  //direction into the wall
         {
             valueFraction()[faceI] = 0.0;
             refGrad()[faceI] = 0.0;
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
index 80f985b3b02941646a11edbbb02b579e91a9d1a5..17b98c5cf57f29c14588f46d7dbcc826df655995 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
@@ -60,14 +60,11 @@ class wideBandDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
-        //- Direction index
-        label myRayIndex_;
+        //- Ray index
+        label rayId_;
 
-        //- Direction index
-        label myWaveLengthIndex_;
-
-        //- Init ray flag
-        label myRayIsInit_;
+        //- Wavelength index
+        label wavelengthId_;
 
         //- Radiative heat flux on walls.
         scalarField qr_;
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
index cbea255a14c9e800c8b6df34225bb79b1251d893..d49d2e596201400f52f860a5cef713e8544ea37b 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
@@ -31,29 +31,28 @@ License
 
 Foam::radiation::blackBodyEmission::blackBodyEmission
 (
-    const fileName& fn,
+    const fileName& name,
     const word& instance,
-    label lambdaj,
+    label nLambda,
     const volScalarField& T
 )
 :
-    blackBodyEmissiveTable_(fn, instance, T.mesh()),
+    blackBodyEmissiveTable_(name, 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),
+    bj_(nLambda),
     T_(T)
 {
-    bj_.setSize(lambdaj);
-    for (label i=0; i < lambdaj; i++)
+    forAll(bj_, lambdaI)
     {
         bj_.set
         (
-            i,
+            lambdaI,
             new volScalarField
             (
                 IOobject
                 (
-                    "bj_" + Foam::name(i) ,
+                    "bj_" + Foam::name(lambdaI) ,
                     T.mesh().time().timeName(),
                     T.mesh(),
                     IOobject::NO_READ,
@@ -80,7 +79,7 @@ Foam::scalar Foam::radiation::blackBodyEmission::flambdaT
     const scalar lambdaT
 ) const
 {
-    return  blackBodyEmissiveTable_.LookUp(lambdaT*1.0e6)[1];
+    return  blackBodyEmissiveTable_.lookUp(lambdaT*1.0e6)[1];
 }
 
 
@@ -114,7 +113,7 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT
     }
     else
     {
-        forAll(T,i)
+        forAll(T, i)
         {
             scalar T1 = flambdaT(band[1]*T[i]);
             scalar T2 = flambdaT(band[0]*T[i]);
@@ -133,11 +132,11 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT
 
 void Foam::radiation::blackBodyEmission::correct
 (
-    label j,
+    const label lambdaI,
     const Vector2D<scalar>& band
 )
 {
-    bj_[j] = EbDeltaLambdaT(T_, band);
+    bj_[lambdaI] = EbDeltaLambdaT(T_, band);
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
index 985dd9a202e8a8ffbb1bc6fb07a54d73bac5c428..a91fed11437054b04a3694297397c10fc2ce1e4f 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
@@ -26,12 +26,11 @@ Class
     Foam::radiation::blackBodyEmission
 
 Description
-    Class Black Body Emission Intensity.
+    Class black body emission
 
 SourceFiles
     blackBodyEmission.C
 
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef blackModyEmission_H
@@ -50,28 +49,35 @@ namespace Foam
 {
 namespace radiation
 {
+
 /*---------------------------------------------------------------------------*\
-                           Class blackBodyEmission Declaration
+                      Class blackBodyEmission Declaration
 \*---------------------------------------------------------------------------*/
+
 class blackBodyEmission
 {
     // Private data
 
         mutable interpolationLookUpTable<scalar> blackBodyEmissiveTable_;
 
-        scalar flambdaT(const scalar lambdaT) const;
-
+        //- Constant C1
         const dimensionedScalar C1_;
 
+        //- Constant C2
         const dimensionedScalar C2_;
 
-        // Ptr List for enregy black body emission
+        // Ptr List for black body emission energy field for each wavelength
         PtrList<volScalarField> bj_;
 
         // Reference to the temperature field
         const volScalarField& T_;
 
 
+    // Private member functions
+
+        scalar flambdaT(const scalar lambdaT) const;
+
+
 public:
 
     // Constructors
@@ -95,9 +101,9 @@ public:
         // Access
 
             //- Black body spectrum
-            inline const volScalarField& bj(label i) const
+            inline const volScalarField& bj(const label lambdaI) const
             {
-                return bj_[i];
+                return bj_[lambdaI];
             }
 
             //- Spectral emission for the black body at T and lambda
@@ -121,7 +127,7 @@ public:
     // Edit
 
         // Update black body emission
-        void correct(label j, const Vector2D<scalar>& band);
+        void correct(const label lambdaI, const Vector2D<scalar>& band);
 };
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
index c0286ca0881035e6e8d72e1824318223e9934a56..faf1153467d3bbbfd9a1d9825d61fb693b1373ed 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
@@ -98,7 +98,6 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
-    aj_(0),
     e_
     (
         IOobject
@@ -125,27 +124,35 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         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)
+    nTheta_(readLabel(coeffs_.lookup("nTheta"))),
+    nPhi_(readLabel(coeffs_.lookup("nPhi"))),
+    nRay_(0),
+    nLambda_(absorptionEmission_->nBands()),
+    aj_(nLambda_),
+    blackBody_
+    (
+        fileName("blackBodyEmissivePower"),
+        mesh_.time().constant(),
+        nLambda_,
+        T
+    ),
+    IRay_(0),
+    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0))
 {
-    aj_.setSize(lambdaj_);
     if (mesh_.nSolutionD() == 3)    //3D
     {
-        RadIntRay_.setSize(4.0*Nphi_*Ntheta_);
-        Ni_ = 4.0*Nphi_*Ntheta_;
-        scalar deltaPhi = mathematicalConstant::pi/(2.0*Nphi_);
-        scalar deltaTheta = mathematicalConstant::pi/Ntheta_;
+        IRay_.setSize(4*nPhi_*nTheta_);
+        nRay_ = 4.0*nPhi_*nTheta_;
+        scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
+        scalar deltaTheta = mathematicalConstant::pi/nTheta_;
         label i = 0;
-        for (label n = 1 ; n <= Ntheta_ ; n++)
+        for (label n = 1; n <= nTheta_; n++)
         {
-            for (label m = 1 ; m <= 4*Nphi_ ; m++)
+            for (label m = 1; m <= 4*nPhi_; m++)
             {
                 scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
-                RadIntRay_.set
+                IRay_.set
                 (
                     i,
                     new radiativeIntensityRay
@@ -154,7 +161,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                         thetai,
                         deltaPhi,
                         deltaTheta,
-                        lambdaj_,
+                        nLambda_,
                         mesh_,
                         absorptionEmission_,
                         blackBody_
@@ -170,14 +177,14 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         {
             scalar thetai = mathematicalConstant::pi/2.0;
             scalar deltaTheta = mathematicalConstant::pi;
-            RadIntRay_.setSize(4.0*Nphi_);
-            Ni_ = 4.0*Nphi_;
-            scalar deltaPhi = mathematicalConstant::pi/(2.0*Nphi_);
+            IRay_.setSize(4*nPhi_);
+            nRay_ = 4.0*nPhi_;
+            scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
             label i = 0;
-            for (label m = 1 ; m <= 4*Nphi_ ; m++)
+            for (label m = 1; m <= 4*nPhi_; m++)
             {
                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
-                RadIntRay_.set
+                IRay_.set
                 (
                     i,
                     new radiativeIntensityRay
@@ -186,7 +193,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                         thetai,
                         deltaPhi,
                         deltaTheta,
-                        lambdaj_,
+                        nLambda_,
                         mesh_,
                         absorptionEmission_,
                         blackBody_
@@ -199,14 +206,14 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         {
             scalar thetai = mathematicalConstant::pi/2.0;
             scalar deltaTheta = mathematicalConstant::pi;
-            RadIntRay_.setSize(2);
-            Ni_ = 2.0;
+            IRay_.setSize(2);
+            nRay_ = 2.0;
             scalar deltaPhi = mathematicalConstant::pi;
             label i = 0;
-            for (label m = 1 ; m <= 2 ; m++)
+            for (label m = 1; m <= 2; m++)
             {
                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
-                RadIntRay_.set
+                IRay_.set
                 (
                     i,
                     new radiativeIntensityRay
@@ -215,7 +222,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                         thetai,
                         deltaPhi,
                         deltaTheta,
-                        lambdaj_,
+                        nLambda_,
                         mesh_,
                         absorptionEmission_,
                         blackBody_
@@ -228,25 +235,26 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     }
 
 
-    // Construct absorption for wave length
-    for(label i=0; i < lambdaj_; i++)
+    // Construct absorption field for each wavelength
+    forAll(aj_, lambdaI)
     {
-        volScalarField* volPtr= new volScalarField
+        aj_.set
         (
-            IOobject
+            lambdaI,
+            new volScalarField
             (
-                "aj_" + Foam::name(i) ,
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            a_
+                IOobject
+                (
+                    "aj_" + Foam::name(lambdaI) ,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                a_
+            )
         );
-
-        aj_.set(i, volPtr);
     }
-
 }
 
 
@@ -264,6 +272,9 @@ bool Foam::radiation::fvDOM::read()
     {
         // nothing to read
 
+//        coeffs_.lookup("nTheta") >> nTheta_;
+//        coeffs_.lookup("nPhi") >> nPhi_;
+
         return true;
     }
     else
@@ -286,22 +297,20 @@ void Foam::radiation::fvDOM::correct()
     updateBlackBodyEmission();
 
     scalar maxResidual = 0;
-    scalar convergenceCriterion = 0;
-    radiationModelCoeffs_.readIfPresent("convergence", convergenceCriterion);
     label radIter = 0;
     do
     {
         radIter ++;
-        for (label i = 0; i < Ni_; i++)
+        forAll(IRay_, rayI)
         {
-            maxResidual = 0;
-            scalar maxBandResidual = RadIntRay_[i].correct(this);
+            maxResidual = 0.0;
+            scalar maxBandResidual = IRay_[rayI].correct(this);
             maxResidual = max(maxBandResidual, maxResidual);
         }
 
         Info << "Radiation solver Iter: " <<  radIter << endl;
 
-    } while(maxResidual > convergenceCriterion);
+    } while(maxResidual > convergence_);
 
     updateG();
 
@@ -311,7 +320,6 @@ void Foam::radiation::fvDOM::correct()
 
 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
 {
-
     return tmp<volScalarField>
     (
         new volScalarField
@@ -348,24 +356,23 @@ Foam::radiation::fvDOM::Ru() const
 
 void Foam::radiation::fvDOM::updateBlackBodyEmission()
 {
-    for (label j=0; j < lambdaj_; j++)
+    for (label j=0; j < nLambda_; j++)
     {
-          blackBody_.correct(j, absorptionEmission_->bands(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++)
+    forAll(IRay_, rayI)
     {
-        RadIntRay_[i].addIntensity();
-        G_ +=  RadIntRay_[i].I()* RadIntRay_[i].omegai();
-        Qr_ += RadIntRay_[i].Qri();
+        IRay_[rayI].addIntensity();
+        G_ += IRay_[rayI].I()*IRay_[rayI].omega();
+        Qr_ += IRay_[rayI].Qr();
     }
 }
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
index 1f3266f15d753af51210e58824b2c1b95dbcdd54..b6139fb4d88787630fc71a0b45a42d23c96c657b 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
@@ -82,38 +82,41 @@ class fvDOM
         //- Incident radiation  [W/m2]
         volScalarField G_;
 
-        //- Total Radiative heat flux [W/m2]
+        //- Total radiative heat flux [W/m2]
         volScalarField Qr_;
 
-        //- Total Absorption coefficient  [1/m]
+        //- Total absorption coefficient [1/m]
         volScalarField a_;
 
-        //- wavelength total Absorption coefficient [1/m]
-        PtrList<volScalarField> aj_;
-
-        //- Total Emission coefficient [1/m]
+        //- Total emission coefficient [1/m]
         volScalarField e_;
 
         //- Emission contribution [Kg/m/s^3]
         volScalarField E_;
 
         //- Number of solid angles in theta
-        label Ntheta_;
+        label nTheta_;
 
         //- Number of solid angles in phi
-        label Nphi_ ;
+        label nPhi_ ;
 
-        //- Total number of directions
-        label Ni_;
+        //- Total number of rays (1 per direction)
+        label nRay_;
 
         //- Number of wavelength bands
-        label lambdaj_;
+        label nLambda_;
 
-        //- Black Body
+        //- Wavelength total absorption coefficient [1/m]
+        PtrList<volScalarField> aj_;
+
+        //- Black body
         blackBodyEmission blackBody_;
 
-        //- List of Pointers to RadiativeIntensityRay
-        PtrList<radiativeIntensityRay> RadIntRay_;
+        //- List of pointers to radiative intensity rays
+        PtrList<radiativeIntensityRay> IRay_;
+
+        //- Convergence criterion
+        scalar convergence_;
 
 
     // Private member functions
@@ -127,8 +130,8 @@ class fvDOM
         //- Update Absorption Coefficients
 //        void updateAbsorptionCoeffs(void);
 
-        //- Update Black Body Emissiom
-        void updateBlackBodyEmission(void);
+        //- Update nlack body emission
+        void updateBlackBodyEmission();
 
 
 public:
@@ -148,7 +151,7 @@ public:
         fvDOM(const volScalarField& T);
 
 
-    // Destructor
+    //- Destructor
     virtual ~fvDOM();
 
 
@@ -174,42 +177,42 @@ public:
 
         // Access
 
-            //- Ray intensity in i direction
-            inline const radiativeIntensityRay& RadIntRay(label i) const;
+            //- Ray intensity for rayI
+            inline const radiativeIntensityRay& IRay(const label rayI) const;
 
-            //- Ray intensity in i direction and j band-width
-            inline const volScalarField& RadIntRayiLambdaj
+            //- Ray intensity for rayI and lambda bandwidth
+            inline const volScalarField& IRayWave
             (
-                const label i,
-                const label j
+                const label rayI,
+                const label lambdaI
             ) const;
 
             //- Number of angles in theta
-            inline label Ntheta() const;
+            inline label nTheta() const;
 
             //- Number of angles in phi
-            inline label Nphi() const;
+            inline label nPhi() const;
 
-            //- Number of directions
-            inline label Ni() const;
+            //- Number of rays
+            inline label nRay() const;
 
             //- Number of wavelengths
-            inline label lambdaj() const;
+            inline label nLambda() const;
 
-            // Const access to a
+            // Const access to total absorption coefficient
             inline const volScalarField& a() const;
 
-            // Const access to aj
-            inline const volScalarField& aj(label i) const;
+            // Const access to wavelength total absorption coefficient
+            inline const volScalarField& aj(const label lambdaI) const;
 
-            // Const access to G
+            // Const access to incident radiation field
             inline const volScalarField& G() const;
 
-            //  Const access to Qr
+            //  Const access to total radiative heat flux field
             inline const volScalarField& Qr() const;
 
-            //  Const access to blackBody
-            virtual const blackBodyEmission& blackBody() const;
+            //  Const access to black body
+            inline const blackBodyEmission& blackBody() const;
 };
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H
index f5fcde4e1326ca7fd5e200b6ff3d3b6b98646178..74e0b2bf345a29f4299e2b1d0137b0f0d51d31a5 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H
@@ -25,43 +25,44 @@ License
 \*---------------------------------------------------------------------------*/
 
 inline const Foam::radiation::radiativeIntensityRay&
-Foam::radiation::fvDOM::RadIntRay(label i) const
+Foam::radiation::fvDOM::IRay(const label rayI) const
 {
-    return  RadIntRay_[i];
+    return  IRay_[rayI];
 }
 
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::RadIntRayiLambdaj
+inline const Foam::volScalarField&
+Foam::radiation::fvDOM::IRayWave
 (
-    const label i,
-    const label j
+    const label rayI,
+    const label lambdaI
 ) const
 {
-    return  RadIntRay_[i].Ilambdaj(j);
+    return IRay_[rayI].IWave(lambdaI);
 }
 
 
-inline Foam::label Foam::radiation::fvDOM::Ntheta() const
+inline Foam::label Foam::radiation::fvDOM::nTheta() const
 {
-    return Ntheta_;
+    return nTheta_;
 }
 
 
-inline Foam::label Foam::radiation::fvDOM::Nphi() const
+inline Foam::label Foam::radiation::fvDOM::nPhi() const
 {
-    return Nphi_;
+    return nPhi_;
 }
 
 
-inline Foam::label Foam::radiation::fvDOM::Ni() const
+inline Foam::label Foam::radiation::fvDOM::nRay() const
 {
-    return Ni_;
+    return nRay_;
 }
 
 
-inline Foam::label Foam::radiation::fvDOM::lambdaj() const
+inline Foam::label Foam::radiation::fvDOM::nLambda() const
 {
-    return lambdaj_;
+    return nLambda_;
 }
 
 
@@ -73,10 +74,10 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::a() const
 
 inline const Foam::volScalarField& Foam::radiation::fvDOM::aj
 (
-    const label i
+    const label lambdaI
 ) const
 {
-    return aj_[i];
+    return aj_[lambdaI];
 }
 
 
@@ -92,7 +93,7 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::Qr() const
 }
 
 
-virtual const Foam::radiation::blackBodyEmission&
+inline const Foam::radiation::blackBodyEmission&
 Foam::radiation::fvDOM::blackBody() const
 {
     return blackBody_;
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
index 20aa9eec885806bb2f5ee3184ebe4df443860763..e9e28658c476f0e6b9e95e5c17d98f186859dcd6 100755
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
@@ -79,18 +79,18 @@ Foam::label Foam::interpolationLookUpTable <Type>::index
 ) const
 {
     label i = 0;
-    label totalindex =
-    min
-    (
-        max
+    label totalIndex =
+        Foam::min
         (
-            label((indice - min_[i])/delta_[i]),
-            0
-        ),
-        dim_[i]
-    );
+            Foam::max
+            (
+                label((indice - min_[i])/delta_[i]),
+                0
+            ),
+            dim_[i]
+        );
 
-    return totalindex;
+    return totalIndex;
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
index 1d108a78135734d449eda7511bd04e04737787c8..d2518d4e02e75d787518e76f95a2f906f5d43128 100755
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
@@ -198,7 +198,7 @@ public:
         inline const List<scalar>& max() const;
 
         //- Return const access to the table name
-        inline const word tableName() const;
+        inline word tableName() const;
 
 
      // Member Operators
@@ -212,6 +212,10 @@ public:
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interpolationLookUpTableI.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTableI.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTableI.H
index 0c2089311e02da392acbedc5738112b71d3803c1..5829db3148e534ecdbcaaeb797427426db55acd0 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTableI.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTableI.H
@@ -84,7 +84,7 @@ Foam::interpolationLookUpTable<Type>::max() const
 
 
 template<class Type>
-inline const Foam::word Foam::interpolationLookUpTable<Type>::tableName() const
+inline Foam::word Foam::interpolationLookUpTable<Type>::tableName() const
 {
     return fileName_.name();
 }
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index 766ff0da4f8911ad07cb465dbb967518454025e2..b9c405126eb4c34720f9a3bcb964a7bbeadd88ab 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -44,11 +44,11 @@ Foam::label Foam::radiation::radiativeIntensityRay::rayId = 0;
 
 Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
 (
-    scalar& phii,
-    scalar& thetai,
-    scalar& deltaPhi,
-    scalar& deltaTheta,
-    label& lambdaj,
+    const scalar phi,
+    const scalar theta,
+    const scalar deltaPhi,
+    const scalar deltaTheta,
+    const label nLambda,
     const fvMesh& mesh,
     const absorptionEmissionModel& absEmmModel,
     const blackBodyEmission& blackBody
@@ -70,7 +70,7 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         mesh_,
         dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
     ),
-    Qri_
+    Qr_
     (
         IOobject
         (
@@ -82,90 +82,62 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         ),
         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
-(
-    const scalar phii,
-    const scalar thetai,
-    const scalar deltaPhi,
-    const scalar deltaTheta,
-    const scalar lambdaj
-)
+    ),
+    d_(vector::zero),
+    dAve_(vector::zero),
+    theta_(theta),
+    phi_(phi),
+    omega_(0.0),
+    nLambda_(nLambda),
+    IWave_(nLambda)
 {
-    phii_ = phii;
-    thetai_ = thetai;
-    nLambdaj_ = lambdaj;
-
-    scalar sinTheta = Foam::sin(thetai);
-    scalar cosTheta = Foam::cos(thetai);
-    scalar sinPhi = Foam::sin(phii);
-    scalar cosPhi = Foam::cos(phii);
-    Si_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
-    omegai_ = 2.0*Foam::sin(thetai)*Foam::sin(deltaTheta/2.0)*deltaPhi;
-    Ilambdaj_.setSize(nLambdaj_);
-    Di_ = vector
+    scalar sinTheta = Foam::sin(theta);
+    scalar cosTheta = Foam::cos(theta);
+    scalar sinPhi = Foam::sin(phi);
+    scalar cosPhi = Foam::cos(phi);
+
+    omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
+    d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
+    dAve_ = vector
     (
         sinPhi
        *Foam::sin(0.5*deltaPhi)
-       *(deltaTheta - Foam::cos(2.0*thetai)
+       *(deltaTheta - Foam::cos(2.0*theta)
        *Foam::sin(deltaTheta)),
         cosPhi
        *Foam::sin(0.5*deltaPhi)
-       *(deltaTheta - Foam::cos(2.0*thetai)
+       *(deltaTheta - Foam::cos(2.0*theta)
        *Foam::sin(deltaTheta)),
-        0.5*deltaPhi*Foam::sin(2.0*thetai)*Foam::sin(deltaTheta)
+        0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
     );
 
-    forAll(Ilambdaj_, i)
+    forAll(IWave_, i)
     {
-        IOobject header
+        IOobject IHeader
         (
             "Ilambda_" + name(rayId) + "_" + name(i),
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
         );
 
         // check if field exists and can be read
-        if (header.headerOk())
+        if (IHeader.headerOk())
         {
-            Ilambdaj_.set
+            IWave_.set
             (
                 i,
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        "Ilambda_" + name(rayId) + "_" + name(i),
-                        mesh_.time().timeName(),
-                        mesh_,
-                        IOobject::MUST_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    mesh_
-                )
+                new volScalarField(IHeader, mesh_)
             );
         }
         else
         {
-            volScalarField Idefault
+            volScalarField IDefault
             (
                 IOobject
                 (
-                    "Idefault",
+                    "IDefault",
                     mesh_.time().timeName(),
                     mesh_,
                     IOobject::MUST_READ,
@@ -174,21 +146,10 @@ void Foam::radiation::radiativeIntensityRay::init
                 mesh_
             );
 
-            Ilambdaj_.set
+            IWave_.set
             (
                 i,
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        "Ilambda_" + name(rayId) + "_"+ name(i),
-                        mesh_.time().timeName(),
-                        mesh_,
-                        IOobject::NO_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    Idefault
-                )
+                new volScalarField(IHeader, IDefault)
             );
         }
     }
@@ -196,31 +157,42 @@ void Foam::radiation::radiativeIntensityRay::init
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 Foam::scalar Foam::radiation::radiativeIntensityRay::correct
 (
     fvDOM* DomPtr
 )
 {
-    Qri_ =  dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    // reset boundary heat flux to zero
+    Qr_ =  dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 
     scalar maxResidual = 0.0;
 
-    for(label i = 0; i < nLambdaj_; i++)
+    forAll(IWave_, lambdaI)
     {
-        volScalarField k = DomPtr->aj(i);
+        volScalarField k = DomPtr->aj(lambdaI);
 
-        volScalarField E = absEmmModel_.ECont(i)/Foam::mathematicalConstant::pi;
+        volScalarField E =
+            absEmmModel_.ECont(lambdaI)/Foam::mathematicalConstant::pi;
 
-        surfaceScalarField Ji = Di_ & mesh_.Sf();
+        surfaceScalarField Ji = dAve_ & mesh_.Sf();
 
-        volScalarField Ib = blackBody_.bj(i)/Foam::mathematicalConstant::pi;
+        volScalarField Ib =
+            blackBody_.bj(lambdaI)/Foam::mathematicalConstant::pi;
 
         fvScalarMatrix IiEq
         (
-            fvm::div(Ji, Ilambdaj_[i], " div(Ji,Ii_h)")
-          + fvm::Sp(k*omegai_, Ilambdaj_[i])
+            fvm::div(Ji, IWave_[lambdaI], " div(Ji,Ii_h)")
+          + fvm::Sp(k*omega_, IWave_[lambdaI])
          ==
-            k*omegai_*Ib + E
+            k*omega_*Ib + E
         );
 
         IiEq.relax();
@@ -234,6 +206,7 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct
         maxResidual = max(eqnResidual, maxResidual);
 
     }
+
     return maxResidual;
 }
 
@@ -242,21 +215,11 @@ void Foam::radiation::radiativeIntensityRay::addIntensity()
 {
     I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 
-    for (label i = 0; i < nLambdaj_; i++)
+    forAll(IWave_, lambdaI)
     {
-        I_ += absEmmModel_.addRadInt(i, Ilambdaj_[i]);
+        I_ += absEmmModel_.addRadInt(lambdaI, IWave_[lambdaI]);
     }
 }
 
 
-void Foam::radiation::radiativeIntensityRay::add
-(
-    const scalarField& qr,
-    const label patchI
-) const
-{
-    Qri_.boundaryField()[patchI] += qr;
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
index c63a8d63d3c46f1d5e8894c58a7e653fc56956ae..d959f45f8fd5eb26de7fc59dc4f8d75866314901 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
@@ -26,6 +26,7 @@ Class
     Foam::radiation::radiativeIntensityRay
 
 Description
+    Radiation intensity for a ray in a given direction
 
 SourceFiles
     radiativeIntensityRay.C
@@ -65,32 +66,32 @@ class radiativeIntensityRay
         //- Black body
         const blackBodyEmission&  blackBody_;
 
-        //- Total radiative intensity in i direction  / [W/m2]
+        //- Total radiative intensity / [W/m2]
         volScalarField I_;
 
-        //- Total radiative heat flux on boundary in i direction
-        mutable volScalarField Qri_;
+        //- Total radiative heat flux on boundary
+        volScalarField Qr_;
 
-        //- Direction i
-        vector Si_;
+        //- Direction
+        vector d_;
 
-        //- Theta angle of direction i
-        scalar thetai_;
+        //- Average direction vector inside the solid angle
+        vector dAve_;
 
-        //- Phi angle of direction i
-        scalar phii_;
+        //- Theta angle
+        scalar theta_;
 
-        //- Solid angle of direction i
-        scalar omegai_;
+        //- Phi angle
+        scalar phi_;
 
-        //- Number of bands on i direction
-        label nLambdaj_;
+        //- Solid angle
+        scalar omega_;
 
-        //- Average vector inside the solid angle
-        vector Di_;
+        //- Number of wavelengths/bands
+        label nLambda_;
 
-        //- List of pointers to radiative intensity wave-length
-        PtrList<volScalarField> Ilambdaj_;
+        //- List of pointers to radiative intensity fields for given wavelengths
+        PtrList<volScalarField> IWave_;
 
 
     // Private member functions
@@ -111,14 +112,14 @@ public:
 
     // Constructors
 
-        //- Null constructor
+        //- Construct form components
         radiativeIntensityRay
         (
-            scalar& phii,
-            scalar& thetai,
-            scalar& deltaPhi,
-            scalar& deltaTheta,
-            label& lambdaj,
+            const scalar phi,
+            const scalar theta,
+            const scalar deltaPhi,
+            const scalar deltaTheta,
+            const label lambda,
             const fvMesh& mesh_,
             const absorptionEmissionModel& absEmmModel_,
             const blackBodyEmission& blackBody
@@ -139,75 +140,48 @@ public:
             //- Initialise the ray in i direction
             void init
             (
-                const scalar phii,
-                const scalar thetai,
+                const scalar phi,
+                const scalar theta,
                 const scalar deltaPhi,
                 const scalar deltaTheta,
                 const scalar lambda
             );
 
-            //- Add radiative heat flux on walls from the boundary patch
-            void add(const scalarField&, const label) const;
-
-            //- Add Radiative intensities from all the bands
+            //- Add radiative intensities from all the bands
             void addIntensity();
 
 
         // Access
 
-            //- Return intensity in i direction
-            inline const volScalarField& I() const
-            {
-                return I_;
-            }
+            //- Return intensity
+            inline const volScalarField& I() const;
+
+            //- Return const access to the boundary heat flux
+            inline const volScalarField& Qr() const;
 
-            //- Return heat flux on boundary in i direction
-            inline const volScalarField& Qri() const
-            {
-                return Qri_;
-            }
+            //- Return non-const access to the boundary heat flux
+            inline volScalarField& Qr();
 
-            //- Return direction i
-            inline const vector Si() const
-            {
-                return Si_;
-            }
+            //- Return direction
+            inline const vector& d() const;
 
             //- Return the average vector inside the solid angle
-            inline const vector Di() const
-            {
-                return Di_;
-            }
-
-            //- Return the number of bands on i direction
-            scalar nLambdaj() const
-            {
-                return  nlambdaj_;
-            }
-
-            //- Return the phi angle of direction i
-            scalar phii() const
-            {
-                return  phii_;
-            }
-
-            //- Return the theta angle of direction i
-            scalar thetai() const
-            {
-                return  thetai_;
-            }
-
-            //- Return the solid angle of direction i
-            scalar omegai() const
-            {
-                return  omegai_;
-            }
-
-            //- Return the list of pointers to radiative intensity wave-length
-            const volScalarField& Ilambdaj(const label i) const
-            {
-                return  Ilambdaj_[i];
-            }
+            inline const vector& dAve() const;
+
+            //- Return the number of bands
+            inline scalar nLambda() const;
+
+            //- Return the phi angle
+            inline scalar phi() const;
+
+            //- Return the theta angle
+            inline scalar theta() const;
+
+            //- Return the solid angle
+            inline scalar omega() const;
+
+            //- Return the radiative intensity for a given wavelength
+            inline const volScalarField& IWave(const label lambdaI) const;
 };
 
 
@@ -218,6 +192,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "radiativeIntensityRayI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C b/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C
index a30d7946f4c94b03eb76aeb056a22d56f3cdbf1a..c7840d9433599e97c0435390dc694b1fb134676f 100644
--- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C
+++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C
@@ -53,7 +53,7 @@ autoPtr<radiationModel> radiationModel::New
             (
                 "radiationProperties",
                 T.mesh().time().constant(),
-                T.mesh().db(),
+                T.mesh().objectRegistry::db(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
@@ -75,7 +75,7 @@ autoPtr<radiationModel> radiationModel::New
             "radiationModel::New(const volScalarField&)"
         )   << "Unknown radiationModel type " << radiationModelTypeName
             << nl << nl
-            << "Valid radiationModel types are :" << nl
+            << "Valid radiationModel types are:" << nl
             << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }
diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
index faf26873f150655f5616408d5e384e11fbb210f2..28e10cd84eca9133bdc7e7c1b41c39fd4e62eebd 100644
--- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
+++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
@@ -54,8 +54,8 @@ Foam::radiation::radiationModel::radiationModel
         IOobject
         (
             "radiationProperties",
-            T.mesh().time().constant(),
-            T.mesh().db(),
+            T.time().constant(),
+            T.mesh().objectRegistry::db(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
@@ -63,7 +63,7 @@ Foam::radiation::radiationModel::radiationModel
     T_(T),
     mesh_(T.mesh()),
     radiation_(lookup("radiation")),
-    radiationModelCoeffs_(subDict(type + "Coeffs")),
+    coeffs_(subDict(type + "Coeffs")),
     nFlowIterPerRadIter_(readLabel(lookup("nFlowIterPerRadIter"))),
     absorptionEmission_(absorptionEmissionModel::New(*this, mesh_)),
     scatter_(scatterModel::New(*this, mesh_))
@@ -83,7 +83,7 @@ bool Foam::radiation::radiationModel::read()
     if (regIOobject::read())
     {
         lookup("radiation") >> radiation_;
-        radiationModelCoeffs_ = subDict(type() + "Coeffs");
+        coeffs_ = subDict(type() + "Coeffs");
 
         return true;
     }
diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
index 4326bc7a4c52aaf40471157289144d2ac4d99812..31e5dd25198345a74f7839cde4248a96b7cc4d28 100644
--- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
+++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
@@ -62,14 +62,13 @@ class absorptionEmissionModel;
 class scatterModel;
 
 /*---------------------------------------------------------------------------*\
-                           Class radiationModel Declaration
+                       Class radiationModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class radiationModel
 :
     public IOdictionary
 {
-
 protected:
 
     // Protected data
@@ -83,13 +82,14 @@ protected:
         //- Model specific dictionary input parameters
         Switch radiation_;
 
-
         //- Radiation model dictionary
-        dictionary radiationModelCoeffs_;
+        dictionary coeffs_;
 
-        //- Number of iteration in the Flow solver per Radiative solver
+        //- Number of iterations in the flow solver per radiation solver
+        //  iteration
         label nFlowIterPerRadIter_;
 
+
         // References to the radiation sub-models
 
             //- Absorption/emission model
@@ -99,7 +99,6 @@ protected:
             autoPtr<scatterModel> scatter_;
 
 
-
 private:
 
     // Private Member Functions
@@ -111,7 +110,6 @@ private:
         void operator=(const radiationModel&);
 
 
-
 public:
 
     //- Runtime type information
@@ -135,20 +133,13 @@ public:
     // Constructors
 
         //- Construct from components
-        radiationModel
-        (
-            const word& type,
-            const volScalarField& T
-        );
+        radiationModel(const word& type, const volScalarField& T);
 
 
     // Selectors
 
          //- Return a reference to the selected radiation model
-         static autoPtr<radiationModel> New
-         (
-             const volScalarField& T
-         );
+         static autoPtr<radiationModel> New(const volScalarField& T);
 
 
     // Destructor
@@ -175,12 +166,7 @@ public:
             virtual tmp<DimensionedField<scalar, volMesh> > Ru() const = 0;
 
             //- Enthalpy source term
-            virtual tmp<fvScalarMatrix> Sh
-            (
-                basicThermo& thermo
-            ) const;
-
-
+            virtual tmp<fvScalarMatrix> Sh(basicThermo& thermo) const;
 };
 
 
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C
index 1ab41afbed222d0b854edaf03693b0bc7e994840..e787a2a8500ff3cc71fe8085149b7e0656bfcf21 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C
+++ b/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(label i) const
+Foam::radiation::absorptionEmissionModel::a(const label bandI) const
 {
-    return aDisp(i) + aCont(i);
+    return aDisp(bandI) + aCont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aCont(label i) const
+Foam::radiation::absorptionEmissionModel::aCont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -89,7 +89,7 @@ Foam::radiation::absorptionEmissionModel::aCont(label i) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::aDisp(label i) const
+Foam::radiation::absorptionEmissionModel::aDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -112,14 +112,14 @@ Foam::radiation::absorptionEmissionModel::aDisp(label i) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::e(label i) const
+Foam::radiation::absorptionEmissionModel::e(const label bandI) const
 {
-    return eDisp(i) + eCont(i);
+    return eDisp(bandI) + eCont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eCont(label i) const
+Foam::radiation::absorptionEmissionModel::eCont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -142,7 +142,7 @@ Foam::radiation::absorptionEmissionModel::eCont(label i) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::eDisp(label i) const
+Foam::radiation::absorptionEmissionModel::eDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -165,14 +165,14 @@ Foam::radiation::absorptionEmissionModel::eDisp(label i) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::E(label i) const
+Foam::radiation::absorptionEmissionModel::E(const label bandI) const
 {
-    return EDisp(i) + ECont(i);
+    return EDisp(bandI) + ECont(bandI);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::ECont(label i) const
+Foam::radiation::absorptionEmissionModel::ECont(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -195,7 +195,7 @@ Foam::radiation::absorptionEmissionModel::ECont(label i) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::absorptionEmissionModel::EDisp(label i) const
+Foam::radiation::absorptionEmissionModel::EDisp(const label bandI) const
 {
     return tmp<volScalarField>
     (
@@ -224,9 +224,8 @@ Foam::label Foam::radiation::absorptionEmissionModel::nBands() const
 
 
 const Foam::Vector2D<Foam::scalar>&
-Foam::radiation::absorptionEmissionModel::bands(label n) const
+Foam::radiation::absorptionEmissionModel::bands(const label n) const
 {
-
     return Vector2D<scalar>::one;
 }
 
@@ -240,11 +239,11 @@ bool Foam::radiation::absorptionEmissionModel::isGrey() const
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::absorptionEmissionModel::addRadInt
 (
-    const label i,
-    const volScalarField& Ilambdaj
+    const label rayI,
+    const volScalarField& IWave
 ) const
 {
-    return Ilambdaj;
+    return IWave;
 }
 
 
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H
index 8fc075c15fa59e2b1b9e4ceccc665af7db188117..5b0c3dbbd7ebfb055001bd8f75335b591d6e63c3 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H
@@ -96,12 +96,11 @@ public:
 
 
     //- Selector
-
-        static autoPtr<absorptionEmissionModel> New
-        (
-            const dictionary& dict,
-            const fvMesh& mesh
-        );
+    static autoPtr<absorptionEmissionModel> New
+    (
+        const dictionary& dict,
+        const fvMesh& mesh
+    );
 
 
     //- Destructor
@@ -124,6 +123,7 @@ public:
                 return dict_;
             }
 
+
             // Absorption coefficient
 
                 //- Absorption coefficient (net)
@@ -166,16 +166,16 @@ public:
 
             //- Const access to the bands - defaults to Vector2D::one for grey
             //  absorption/emission
-            virtual const Vector2D<scalar>& bands(label n) const;
+            virtual const Vector2D<scalar>& bands(const label n) const;
 
             //- Flag for whether the absorption/emission is for a grey gas
-            virtual bool isGrey(void) const;
+            virtual bool isGrey() const;
 
-            //- Add radiative intensity fir ray i
+            //- Add radiative intensity for ray i
             virtual tmp<volScalarField> addRadInt
             (
-                const label i,
-                const volScalarField& Ilambdaj
+                const label rayI,
+                const volScalarField& IWave
             ) const;
 
             //- Correct absorption coefficients
@@ -183,7 +183,7 @@ public:
             (
                 volScalarField& a,
                 PtrList<volScalarField>& aj
-            ) const ;
+            ) const;
 };
 
 
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 898edefaab8fb1dac2d17c75086ac3c30fe9374d..09435f352e891c0d32da78e430e7013a7d7d5c70 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -63,7 +63,7 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
         mesh.time().constant(),
         mesh
     ),
-    thermo_(mesh.db().lookupObject<basicThermo>("thermophysicalProperties")),
+    thermo_(mesh.lookupObject<basicThermo>("thermophysicalProperties")),
     EhrrCoeff_(readScalar(coeffsDict_.lookup("EhrrCoeff"))),
     Yj_(0)
 {
@@ -91,7 +91,7 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
     label j = 0;
     forAllConstIter(HashTable<label>, speciesNames_, iter)
     {
-        if (mesh.db().foundObject<volScalarField>("ft"))
+        if (mesh.foundObject<volScalarField>("ft"))
         {
             if (lookUpTable_.found(iter.key()))
             {
@@ -102,15 +102,15 @@ Foam::radiation::greyMeanAbsorptionEmission::greyMeanAbsorptionEmission
 
                 specieIndex_[iter()] = index;
             }
-            else if (mesh.db().foundObject<volScalarField>(iter.key()))
+            else if (mesh.foundObject<volScalarField>(iter.key()))
             {
                 volScalarField& Y =
                     const_cast<volScalarField&>
                     (
-                        mesh.db().lookupObject<volScalarField>(iter.key())
+                        mesh.lookupObject<volScalarField>(iter.key())
                     );
                 Yj_.set(j, &Y);
-                specieIndex_[iter()] = 0.0;
+                specieIndex_[iter()] = 0;
                 j++;
                 Info << "specie: " << iter.key() << " is being solved "
                      << endl;
@@ -154,8 +154,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
 {
     const volScalarField& T = thermo_.T();
     const volScalarField& p = thermo_.p();
-    const volScalarField& ft =
-        this->mesh().db().lookupObject<volScalarField>("ft");
+    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
 
     label nSpecies = speciesNames_.size();
 
@@ -234,7 +233,7 @@ Foam::radiation::greyMeanAbsorptionEmission::eCont(const label bandI) const
                 IOobject::NO_WRITE
             ),
             mesh(),
-            dimensionedScalar("e",dimless/dimLength, 0.0)
+            dimensionedScalar("e", dimless/dimLength, 0.0)
         )
     );
 
@@ -252,21 +251,20 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
             IOobject
             (
                 "E",
-                mesh().time().timeName(),
-                mesh(),
+                mesh_.time().timeName(),
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            mesh(),
+            mesh_,
             dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
         )
     );
 
-    if (mesh().db().foundObject<volScalarField>("hrr"))
+    if (mesh_.foundObject<volScalarField>("hrr"))
     {
-        const volScalarField& hrr =
-            mesh().db().lookupObject<volScalarField>("hrr");
-        E().internalField() =  EhrrCoeff_*hrr.internalField();
+        const volScalarField& hrr = mesh_.lookupObject<volScalarField>("hrr");
+        E().internalField() = EhrrCoeff_*hrr.internalField();
     }
 
     return E;
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index fa959e210c0e805b73759e881d7f8a41729fa510..d9f4223b2b5c648d2913fbc16e25ba29253860d7 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
@@ -63,7 +63,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
         mesh.time().constant(),
         mesh
     ),
-    thermo_(mesh.db().lookupObject<basicThermo>("thermophysicalProperties")),
+    thermo_(mesh.lookupObject<basicThermo>("thermophysicalProperties")),
     Yj_(0),
     totalWaveLength_(0)
 {
@@ -85,7 +85,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
 
         label nSpec = 0;
         const dictionary& specDicts = dict.subDict("species");
-        forAllConstIter(dictionary, SpecDicts, iter)
+        forAllConstIter(dictionary, specDicts, iter)
         {
             const word& key = iter().keyword();
             if (nBand == 0)
@@ -124,10 +124,10 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
                 << " with index: " << index << endl;
             specieIndex_[iter()] = index;
         }
-        else if (mesh.db().foundObject<volScalarField>(iter.key()))
+        else if (mesh.foundObject<volScalarField>(iter.key()))
         {
             volScalarField& Y = const_cast<volScalarField&>
-                (mesh.db().lookupObject<volScalarField>(iter.key()));
+                (mesh.lookupObject<volScalarField>(iter.key()));
 
             Yj_.set(j, &Y);
 
@@ -164,8 +164,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
 {
     const volScalarField& T = thermo_.T();
     const volScalarField& p = thermo_.p();
-    const volScalarField& ft =
-        this->mesh().db().lookupObject<volScalarField>("ft");
+    const volScalarField& ft = mesh_.lookupObject<volScalarField>("ft");
 
     label nSpecies = speciesNames_.size();
 
@@ -275,10 +274,10 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
         )
     );
 
-    if (mesh().db().foundObject<volScalarField>("hrr"))
+    if (mesh().foundObject<volScalarField>("hrr"))
     {
         const volScalarField& hrr =
-            mesh().db().lookupObject<volScalarField>("hrr");
+            mesh().lookupObject<volScalarField>("hrr");
         E().internalField() =
             iEhrrCoeffs_[bandI]
            *hrr.internalField()