diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
index 1b50a08cf..81547766e 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/absorptionCoeffs/absorptionCoeffs.C
@@ -40,7 +40,7 @@ void Foam::radiation::absorptionCoeffs::checkT(const scalar T) const
     if (T < Tlow_ || T > Thigh_)
     {
         WarningInFunction
-            << "usinf absCoeff out of temperature range:" << nl
+            << "using absorptionCoeffs out of temperature range:" << nl
             << "    " << Tlow_ << " -> " << Thigh_ << ";  T = " << T
             << nl << endl;
     }
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
index b2fdb85d1..30b3467a1 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
@@ -197,6 +197,51 @@ Foam::scalar Foam::radiation::blackBodyEmission::fLambdaT
 }
 
 
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::blackBodyEmission::DeltaLambdaT
+(
+    const volScalarField& T,
+    const Vector2D<scalar>& band
+) const
+{
+    tmp<volScalarField> DeltaLambdaT
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Eb",
+                T.mesh().time().timeName(),
+                T.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+	    T.mesh(),
+	    dimensionedScalar("", dimless, 1.0)
+        )
+    );
+
+
+    if (band == Vector2D<scalar>::one)
+    {
+        return DeltaLambdaT;
+    }
+    else
+    {
+        scalarField& DeltaLambdaTf = DeltaLambdaT.ref();
+
+        forAll(T, i)
+        {
+            scalar T1 = fLambdaT(band[1]*T[i]);
+            scalar T2 = fLambdaT(band[0]*T[i]);
+            DeltaLambdaTf[i] *= T1 - T2;
+        }
+
+        return DeltaLambdaT;
+    }
+}
+
+
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::blackBodyEmission::EbDeltaLambdaT
 (
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
index d7f175fc5..ae9cee868 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.H
@@ -130,6 +130,12 @@ public:
                 const Vector2D<scalar>& band
             ) const;
 
+            //- Proportion of total energy at T from lambda1 to lambda2
+            tmp<Foam::volScalarField> DeltaLambdaT
+            (
+                const volScalarField& T,
+                const Vector2D<scalar>& band
+            ) const;
 
     // Edit
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
index 94d365081..b20d70781 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
@@ -455,7 +455,8 @@ void Foam::radiation::fvDOM::calculate()
 
 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
 {
-    return tmp<volScalarField>
+    // Construct using contribution from first frequency band
+    tmp<volScalarField> tRp
     (
         new volScalarField
         (
@@ -468,28 +469,73 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
                 IOobject::NO_WRITE,
                 false
             ),
-            // Only include continuous phase emission
-            4*absorptionEmission_->aCont()*physicoChemical::sigma
+	    (
+	        4*physicoChemical::sigma
+	      * (aLambda_[0]-absorptionEmission_->aDisp(0)())
+	      * blackBody_.DeltaLambdaT(T_,absorptionEmission_->bands(0))
+	    )
         )
     );
+
+    volScalarField& Rp=tRp.ref();
+
+    // Add contributions over remaining frequency bands
+    for (label j=1; j < nLambda_; j++)
+    {
+        Rp +=
+	(
+	    4*physicoChemical::sigma
+	  * (aLambda_[j]-absorptionEmission_->aDisp(j)())
+	  * blackBody_.DeltaLambdaT(T_,absorptionEmission_->bands(j))
+	);
+    }
+    
+    return tRp;
 }
 
 
 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
 Foam::radiation::fvDOM::Ru() const
 {
+    tmp<DimensionedField<scalar, volMesh>> tRu
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "Ru",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+	    mesh_,
+	    dimensionedScalar("zero",dimensionSet(1,-1,-3,0,0), 0.0)
+        )
+    );
 
-    const volScalarField::Internal& G =
-        G_();
+    DimensionedField<scalar, volMesh>& Ru=tRu.ref();
 
-    const volScalarField::Internal E =
-        absorptionEmission_->ECont()()();
+    // Sum contributions over all frequency bands
+    for (label j=0; j < nLambda_; j++)
+    {
+	// Compute total incident radiation within frequency band
+	tmp<DimensionedField<scalar, volMesh>> Gj
+	(
+	    IRay_[0].ILambda(j)()*IRay_[0].omega()
+	);
 
-    // Only include continuous phase absorption
-    const volScalarField::Internal a =
-        absorptionEmission_->aCont()()();
+	for (label rayI=1; rayI < nRay_; rayI++)
+        {
+	    Gj.ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
+	}
+	
+	Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
+	    - absorptionEmission_->ECont(j)()();
+    }
 
-    return a*G - E;
+    return tRu;
 }