diff --git a/src/forward.jl b/src/forward.jl
index 77d15eb..76fa3fa 100644
--- a/src/forward.jl
+++ b/src/forward.jl
@@ -393,7 +393,7 @@ function alifForward(   zit,
                         # compute epsilonRec
                         epsilonRec[i1,i2,i3,i4] = (alpha[i1,i2,i3,i4] * epsilonRec[i1,i2,i3,i4])
 
-                        # compute epsilonRecA
+                        # compute epsilonRecA   use eq.26
                         epsilonRecA[i1,i2,i3,i4] = (rho[i1,i2,i3,i4] * 
                                                         (phi[i1,i2,i3,i4] * epsilonRec[i1,i2,i3,i4]))
 
@@ -424,14 +424,13 @@ function alifForward(   zit,
                         phi[i1,i2,i3,i4] = (gammaPd[i1,i2,i3,i4] / vth[i1,i2,i3,i4]) * 
                                 max(0, 1 - ((vt[i1,i2,i3,i4] - vth[i1,i2,i3,i4]) / vth[i1,i2,i3,i4]))
 
+                        # compute epsilonRec
+                        epsilonRec[i1,i2,i3,i4] = (alpha[i1,i2,i3,i4] * epsilonRec[i1,i2,i3,i4]) + 
+                                                        (zit[i1,i2,i3,i4] * subscription[i1,i2,i3,i4])
                         # compute epsilonRecA   use eq.26
                         epsilonRecA[i1,i2,i3,i4] = (rho[i1,i2,i3,i4] * 
                                                         (phi[i1,i2,i3,i4] * epsilonRec[i1,i2,i3,i4])) + 
                                                         (zit[i1,i2,i3,i4] * subscription[i1,i2,i3,i4])
-
-                        # compute epsilonRec
-                        epsilonRec[i1,i2,i3,i4] = (alpha[i1,i2,i3,i4] * epsilonRec[i1,i2,i3,i4]) + 
-                                                        (zit[i1,i2,i3,i4] * subscription[i1,i2,i3,i4])
                 end
         end
         return nothing