diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index c7f938fb4d1ab36cf84d8660c3bfe007135bb363..39e309dd5ab213f68124b402383f0a1ac5ecb35c 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2019OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -193,13 +193,23 @@ updateCoeffs()
 
     const vector& myRayId = dom.IRay(rayId).d();
 
-    // Use updated Ir while iterating over rays
-    // avoids to used lagged qin
-    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
-
-    for (label rayI=1; rayI < dom.nRay(); rayI++)
+    scalarField Ir(patch().size(), Zero);
+    forAll(Iw, facei)
     {
-        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
+        for (label rayi=0; rayi < dom.nRay(); rayi++)
+        {
+            const vector& d = dom.IRay(rayi).d();
+
+            if ((-n[facei] & d) < 0.0)
+            {
+                // q into the wall
+                const scalarField& IFace =
+                    dom.IRay(rayi).ILambda(lambdaId).boundaryField()[patchi];
+
+                const vector& rayDave = dom.IRay(rayi).dAve();
+                Ir[facei] += IFace[facei]*(n[facei] & rayDave);
+            }
+        }
     }
 
     if (dom.useSolarLoad())
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
index 148c582870e50ddd8802fd51e99662c0acf9cb1e..b9f061fde10d7c24b6570a0813703bd27b027939 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -188,17 +188,6 @@ updateCoeffs()
     scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
     scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
-    // Use updated Ir while iterating over rays
-    // avoids to used lagged qin
-    /*
-    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
-
-    for (label rayI=1; rayI < dom.nRay(); rayI++)
-    {
-        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
-    }
-    */
-
     // Calculate Ir into the wall on the same lambdaId
     scalarField Ir(patch().size(), Zero);
     forAll(Iw, facei)
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index 89c6942cc7275a1eefaa7330de13d9087091e8bf..12bd0c2f69d0cacfc08dfb231e15c80e3eaea23f 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2019 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -270,6 +270,8 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 {
     // Reset boundary heat flux to zero
     qr_.boundaryFieldRef() = 0.0;
+    qem_.boundaryFieldRef() = 0.0;
+    qin_.boundaryFieldRef() = 0.0;
 
     scalar maxResidual = -GREAT;