diff --git a/diffusivityModel/Ion/Ion.H b/diffusivityModel/Ion/Ion.H
index 4db7a0d..e047bd3 100644
--- a/diffusivityModel/Ion/Ion.H
+++ b/diffusivityModel/Ion/Ion.H
@@ -131,6 +131,8 @@ public:
 
         inline scalar Zrot(scalar T) const;
 
+        inline bool hasRct(const word &target) const;
+
         inline scalar Arct(const word &target) const;
 
         inline scalar Brct(const word &target) const;
diff --git a/diffusivityModel/Ion/IonI.H b/diffusivityModel/Ion/IonI.H
index 0dddb69..30d95f3 100644
--- a/diffusivityModel/Ion/IonI.H
+++ b/diffusivityModel/Ion/IonI.H
@@ -73,6 +73,13 @@ inline Foam::scalar Foam::Ion::Zrot(const scalar T)
 }
 
 
+inline bool Foam::Ion::hasRct(const word &target)
+    const
+{
+    return rctTarget_.contains(target);
+}
+
+
 inline Foam::scalar Foam::Ion::Arct(const word &target)
     const
 {
diff --git a/diffusivityModel/N64/N64.C b/diffusivityModel/N64/N64.C
index b52b1a6..2f37a3a 100644
--- a/diffusivityModel/N64/N64.C
+++ b/diffusivityModel/N64/N64.C
@@ -26,18 +26,20 @@ License
 #include "N64.H"
 
 #include "error.H"
+#include "scalarMatrices.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::scalar Foam::N64::a11[(M+1)*(M+2)/2] = {
-        0., -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
+        1.21327441805, -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
         -0.02385399, -0.14337633, -0.18493828, 0.05466979, -0.0010004, 0.01936938,
         -0.0130214, 0.12179535, -0.05237791, -0.00172218, -0.00447904, 0.02211965,
         0.00964791, -0.07566612, 0.02527412
     };
 
+
 const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
-        0.00000000e+00, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
+        1.3678302569, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
         3.59026813e-01, -3.25082644e-01, -8.27621218e-03, -1.12825149e-01,
         -2.07675933e-01, 1.02433682e+00, -4.59772292e-04, 1.03264160e-02,
         -1.35152522e-02, 1.69188798e-01, -1.21839345e+00, -1.88957374e-03,
@@ -45,12 +47,19 @@ const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
         5.07723424e-01
     };
 
-Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
+//Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
+Foam::scalarField Foam::N64::vec((M+1), 0.0);
 
 const Foam::scalar Foam::N64::digamma3 = 0.922784335098467139393487909917597568957840664060;
 
 const Foam::scalar Foam::N64::trigamma3 = sqr(Ion::pi) / 6. - 5./4.;
 
+/*
+const Foam::scalar Foam::N64::digamma3 = 0.922784335098;
+
+const Foam::scalar Foam::N64::trigamma3 = 0.394934066848;
+*/
+
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 
@@ -71,15 +80,89 @@ Foam::N64::N64(const Ion& a, const Neutral& b)
     eij_(5.2 * b_.alpha() * sqr(a_.z()) * (1.0 + ksiij_) / pow4(sigmaij_) / (a_.k/a_.e)),
     C6ij_(0.0),
     gammaij_(0.0),
-    Wij_(1.0/(1.0/a_.W() + 1.0/b_.W()))
+    Wij_(1.0/(1.0/a_.W() + 1.0/b_.W())),
+    b11_((M+1), 0.0),
+    b22_((M+1), 0.0)
 {
     scalar aij = a_.alpha() / b_.alpha();
     scalar aji = 1.0 / aij;
     C6ij_ = 2.0 * a_.C6() * b_.C6() / (aji * a_.C6() + aij * b_.C6());
     gammaij_ = ((2.0/sqr(a_.z())) * C6ij_ + b_.alphaQ()) / (b_.alpha() * sqr(sigmaij_));
 
-    Info << "N64, gamma = " << gammaij_ << endl;
+    Info << a_.name() + " = " + b_.name() + " N64, gamma = " << gammaij_ << endl;
 
+    scalarField temp((M+1)*(M+2)/2, 0.0); // temporary feature vector with fixed gammaij value
+
+    transform(1.0, gammaij_, temp.data());
+
+    Info << "fixed gammaij feature" << endl << temp << endl;
+
+    scalarField t11 (temp * UList<scalar>(const_cast<scalar*>(a11), temp.size()));
+
+    scalarField t22 (temp * UList<scalar>(const_cast<scalar*>(a22), temp.size()));
+
+    scalarSquareMatrix mat11(M+1, 0.0);
+
+    scalarSquareMatrix mat22(M+1, 0.0);
+
+    label index = 0;
+
+    for (label s = 0; s < M+1; s++)
+    {
+        for (label i = 0; i <= s; i++)
+        {
+           mat11(i, s-i) = t11[index];
+           mat22(i, s-i) = t22[index];
+           index++;
+        }
+    }
+
+
+    for (label s = 0; s < M+1; s++)
+    {
+        b11_ += UList<scalar>(mat11[s], M+1);
+        b22_ += UList<scalar>(mat22[s], M+1);
+    }
+
+    if (false) //  (a.hasRct(b.name()))
+    {
+        Info << "RCT" << endl;
+
+        scalar A = a.Arct(b.name());
+
+        scalar B = a.Brct(b.name());
+
+        scalar C = A/B - digamma3 - loge(eij_);
+
+        scalar D = sqr(B) / (a_.pi * sqr(sigmaij_)); // Formula in [A^2] (Angstrom)
+
+        scalar a0 = D * (trigamma3 + sqr(C));
+
+        scalar a1 = D * (-2.0 * C);
+
+        scalar a2 = D;
+
+        b11_[0] += a0;
+
+        b11_[1] += a1;
+
+        b11_[2] += a2;
+    }
+
+    /*
+
+    for (label s = 0; s < M+1; s++)
+    {
+        Info << UList<scalar>(mat11[s], M+1) << endl;
+    }
+    Info << b11_ << endl;
+    for (label s = 0; s < M+1; s++)
+    {
+        Info << UList<scalar>(mat22[s], M+1) << endl;
+    }
+    Info << b22_ << endl;
+
+    */
 }
 
 
diff --git a/diffusivityModel/N64/N64.H b/diffusivityModel/N64/N64.H
index 37dd4e2..dc4a904 100644
--- a/diffusivityModel/N64/N64.H
+++ b/diffusivityModel/N64/N64.H
@@ -93,11 +93,20 @@ class N64
         //- Combined molecular weight
         scalar Wij_;
 
+        //- Omega(1,1) Fit Coefficient with fixed gamma_ij
+        scalarField b11_;
+
+        //- Omega(2,2) Fit Coefficient with fixed gamma_ij
+        scalarField b22_;
+
 
         //- Construct null
         static const scalar a11[(M+1)*(M+2)/2];
+
         static const scalar a22[(M+1)*(M+2)/2];
-        static scalar vec[(M+1)*(M+2)/2];
+
+        // static scalar vec[(M+1)*(M+2)/2];
+        static scalarField vec;
 
         //- Construct null
         static constexpr scalar K1 = 1.767;
@@ -125,9 +134,28 @@ class N64
 
         // primary template
         template <label Order>
-            class BivarPolyFeat
+            struct UnivarPolyFeat
+            {
+                    static scalar transform (scalar x1, scalar *vec)
+                    {
+                        scalar highest = UnivarPolyFeat<Order-1>::transform(x1,vec) * x1;
+
+                        *(vec+Order) = highest;
+
+                        return highest;
+                    }
+            };
+
+        // convenience function
+        inline void transform (scalar x1, scalar *vec)
+        {
+            UnivarPolyFeat<M>::transform(x1,vec);
+        }
+
+        // primary template
+        template <label Order>
+            struct BivarPolyFeat
             {
-                public:
                     static scalar *transform (scalar x1, scalar x2, scalar *vec)
                     {
                         scalar *src = BivarPolyFeat<Order-1>::transform(x1,x2,vec);
@@ -149,6 +177,8 @@ class N64
             BivarPolyFeat<M>::transform(x1,x2,vec);
         }
 
+
+        // primary template
         template <label N>
             struct ScalarProduct
             {
@@ -229,6 +259,17 @@ public:
 };
 
 // partial specialization as end criteria
+template <>
+class N64::UnivarPolyFeat<0>
+{
+    public:
+        static scalar transform (scalar x1, scalar *vec)
+        {
+            *vec = 1.0;
+            return 1.0;
+        }
+};
+
 template <>
 class N64::BivarPolyFeat<0>
 {
@@ -241,10 +282,10 @@ class N64::BivarPolyFeat<0>
 };
 
 template <>
-struct N64::ScalarProduct<0>
+struct N64::ScalarProduct<1>
 {
     static scalar product (const scalar* first, const scalar* second) {
-        return 0.0;
+        return *first * *second;
     };
 };
 
diff --git a/diffusivityModel/N64/N64I.H b/diffusivityModel/N64/N64I.H
index 99ac9fe..1c0ed05 100644
--- a/diffusivityModel/N64/N64I.H
+++ b/diffusivityModel/N64/N64I.H
@@ -59,15 +59,15 @@ inline Foam::scalar Foam::N64::bivarPolyFeat(const scalar Tstar)
 
 inline Foam::scalar Foam::N64::Omega11(const scalar Tstar)
 {
-    transform (loge(Tstar), gammaij_, vec);
-    return scalarProduct (a11, vec);
+    transform (loge(Tstar), vec.data());
+    return ScalarProduct<M+1>::product (b11_.data(), vec.data());
 }
 
 
 inline Foam::scalar Foam::N64::Omega22(const scalar Tstar)
 {
-    transform (loge(Tstar), gammaij_, vec);
-    return scalarProduct (a22, vec);
+    transform (loge(Tstar), vec.data());
+    return ScalarProduct<M+1>::product (b22_.data(), vec.data());
 }