diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 86172659b86d597637d074394668e0b5cd6b81ab..f2178b30be2a0521866fc8a2666b99bb58325e67 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -10,40 +10,97 @@ SUBROUTINE moments_eq_rhs
   use prec_const
   IMPLICIT NONE
 
-  INTEGER:: ip,ij, ipj, ikr,ikz
-  REAL(dp) :: ip_dp, ij_dp
-  REAL(dp) :: tmp, moml, sqrtT
-  
+  INTEGER     :: ip,ij, ipj, ipp2j, ipm2j, ipjp1, ipjm1, ikr,ikz
+  INTEGER     :: kronp0, kronp2
+  INTEGER, DIMENSION(2) :: tmppj
+  REAL(dp)    :: ip_dp, ij_dp
+  REAL(dp)    :: kz, taueqeetaBi, tauiqietaBi, tauaqaetaBi
+  COMPLEX(dp) :: Napj, Napp2j, Napm2j, Napjp1, Napjm1
+   
+  taueqeetaBi = tau_e/real(q_e)*eta_B * imagu ! Factor tau_s/qu_s*eta_B*i
+  tauiqietaBi = tau_i/real(q_i)*eta_B * imagu
+
+  Napp2j = 0; Napm2j = 0;  Napjp1 = 0; Napjm1 = 0; ! higher mixing moments terms
+  kronp0 = 0; kronp2 = 0; ! Useful kronecker for phi terms
+
+  ipp2j   = 1; ipm2j   = 1; ipjp1   = 1; ! Mixing moment indices are set to one initially
+  ipjm1   = 1;                           ! in order to not overpass the moment array
+
   do ipj = ipjs, ipje
 
     if (ipj .le. Nmome + 1) then ! electrons moments
-      CALL rabe(ipj, ip, ij) ! Compute p,j electrons moments degree
+
+      tmppj   = rabe(ipj) ! Compute p,j electrons moments degree
+      ip = tmppj(1); ij = tmppj(2);
+      if (ip+2 .le. pmaxe+1) then
+        ipp2j   = bare(ip+2,ij)
+        Napp2j = moments(ipp2j,ikr,ikz,updatetlevel)
+      endif
+      if (ip-2 .ge. 1) then
+        ipm2j   = bare(ip-2,ij)
+        Napm2j = moments(ipm2j,ikr,ikz,updatetlevel)
+      endif
+      if (ij+1 .le. jmaxe+1) then
+        ipjp1   = bare(ip,ij+1)
+        Napjp1 = moments(ipjp1,ikr,ikz,updatetlevel)
+      endif
+      if (ij-1 .ge. 1) then
+        ipjm1   = bare(ip,ij-1)
+        Napjm1 = moments(ipjm1,ikr,ikz,updatetlevel)
+      endif
+      tauaqaetaBi = taueqeetaBi !species dependant factor
+
     else ! ions moments
-      CALL rabi(ipj, ip, ij) ! Compute p,j ions moments degree
+
+      tmppj = rabi(ipj) ! Compute p,j ions moments degree
+      ip = tmppj(1); ij = tmppj(2);
+      if (ip+2 .le. pmaxi+1) then
+        ipp2j   = bari(ip+2,ij)
+        Napp2j = moments(ipp2j,ikr,ikz,updatetlevel)
+      endif
+      if (ip-2 .ge. 1) then
+        ipm2j   = bari(ip-2,ij)
+        Napm2j = moments(ipm2j,ikr,ikz,updatetlevel)
+      endif
+      if (ij+1 .le. jmaxi+1) then
+        ipjp1   = bari(ip,ij+1)
+        Napjp1 = moments(ipjp1,ikr,ikz,updatetlevel)
+      endif
+      if (ij-1 .ge. 1) then
+        ipjm1   = bari(ip,ij-1)
+        Napjm1 = moments(ipjm1,ikr,ikz,updatetlevel)
+      endif
+      tauaqaetaBi = tauiqietaBi
+    endif
+
+    if (ip .eq. 0) then
+      kronp0 = 1
     endif
+    if (ip .eq. 2) then
+      kronp2 = 1
+    endif 
+
+    Napj   = moments(ipj,ikr,ikz,updatetlevel)
 
     ip_dp = real(ip,dp) 
     ij_dp = real(ij,dp) 
 
     do ikr = ikrs,ikre 
       do ikz = ikzs,ikze
+        kz = kzarray(ikz)
+      ! N_a^{pj} term
+        Napj   = 2*(ip_dp + ij_dp + 1._dp) * Napj
+      ! N_a^{p+2,j} term
+        Napp2j = SQRT(ip_dp + 1._dp) * SQRT(ip_dp + 2._dp) * Napp2j
+      ! N_a^{p-2,j} term
+        Napm2j = SQRT(ip_dp) * SQRT(ip_dp - 1._dp) * Napm2j
+      ! N_a^{p,j+1} term
+        Napjp1 = -(ij_dp + 1._dp) * Napjp1
+      ! N_a^{p,j-1} term
+        Napjm1 = -ij_dp * Napjm1
 
-      ! N_e^{p} term
-
-      ! N_e^{p+1} term
-
-
-      ! N_e^{p-1} term
-
-
-      ! N_e^{p-2} term
-
-    end do
-  end do
-  
-    do ikr = ikrs,ikre
-      do ikz = ikzs,ikze
-        moments_rhs(ipj,ikr,ikz,updatetlevel) = 0._dp ! debug
+        moments_rhs(ipj,ikr,ikz,updatetlevel) = &
+         tauaqaetaBi * kz * (Napj + Napp2j + Napm2j + Napjp1 + Napjm1)
       end do
     end do
   end do