From 718b267ee981e6426683ec392bfc0caa9414f4c7 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 10 Aug 2018 16:39:35 -0500
Subject: [PATCH 01/18] Add math functions to calculate normalization arrays

---
 openmc/capi/math.py | 93 +++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 89 insertions(+), 4 deletions(-)

diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 390ee466ff..697ed558d7 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -13,7 +13,10 @@ _dll.calc_pn_c.restype = None
 _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
 _dll.evaluate_legendre_c.restype = c_double
-_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
+_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double, c_int]
+
+_dll.calc_norm_pn_c.restype = None
+_dll.calc_norm_pn_c.argtypes =[c_int, c_double, c_double, ndpointer(c_double)]
 
 _dll.calc_rn_c.restype = None
 _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
@@ -24,6 +27,12 @@ _dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
 _dll.calc_zn_rad_c.restype = None
 _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
+_dll.evaluate_zernike_rad_c.restype = c_double
+_dll.evaluate_zernike_rad_c.argtypes = [c_int, POINTER(c_double), c_double]
+
+_dll.calc_norm_zn_rad_c.restype = None
+_dll.calc_norm_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
+
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
                                 POINTER(c_double)]
@@ -81,7 +90,7 @@ def calc_pn(n, x):
     return pnx
 
 
-def evaluate_legendre(data, x):
+def evaluate_legendre(data, x, flag=1):
     """ Finds the value of f(x) given a set of Legendre coefficients
     and the value of x.
 
@@ -91,6 +100,8 @@ def evaluate_legendre(data, x):
         Legendre coefficients
     x : float
         Independent variable to evaluate the Legendre at
+    flag : int
+        Flags to indicate use default normalization or not
 
     Returns
     -------
@@ -102,7 +113,31 @@ def evaluate_legendre(data, x):
     data_arr = np.array(data, dtype=np.float64)
     return _dll.evaluate_legendre_c(len(data),
                                     data_arr.ctypes.data_as(POINTER(c_double)),
-                                    x)
+                                    x, flag)
+
+
+def calc_norm_pn(l, zmin, zmax):
+    """ Calculate normalization arrays for n orders in Legendre polynomial.
+
+    Parameters
+    ----------
+    l : int
+        Maximum order
+    zmin : float
+        The lowest boundary of the domain
+    zmax : float
+        The highest boundary of the domain
+
+    Returns
+    -------
+    numpy.ndarray
+        Corresponding resulting list of normalization coefficients
+
+    """
+
+    norm_pn = np.empty(l + 1, dtype=np.float64)
+    _dll.calc_norm_pn_c(l, zmin, zmax, norm_pn)
+    return norm_pn
 
 
 def calc_rn(n, uvw):
@@ -182,7 +217,57 @@ def calc_zn_rad(n, rho):
     zn_rad = np.zeros(num_bins, dtype=np.float64)
     _dll.calc_zn_rad_c(n, rho, zn_rad)
     return zn_rad
-    
+
+
+def evaluate_zernike_rad(data, r):
+    """ Finds the value of f(x) given a set of even order Zernike coefficients
+    and the value of r.
+
+    Parameters
+    ----------
+    data : iterable of float
+        Legendre coefficients
+    r : float
+        Independent variable to evaluate the radial part of Zernike at
+
+    Returns
+    -------
+    float
+        Corresponding Zernike expansion result
+
+    """
+
+    data_arr = np.array(data, dtype=np.float64)
+    n = int(2 * (len(data) - 1))
+    return _dll.evaluate_zernike_rad_c(n,
+                                        data_arr.ctypes.data_as(POINTER(
+                                        c_double)),
+                                        r)
+
+
+def calc_norm_zn_rad(n, rho):
+    """ Calculate normalization arrays for the even orders in modified Zernike
+    polynomial moment with no azimuthal dependency (m=0).
+
+    Parameters
+    ----------
+    n : int
+        Maximum order
+    rho : float
+        Radial location in the unit disk
+
+    Returns
+    -------
+    numpy.ndarray
+        Corresponding resulting list of normalization coefficients
+
+    """
+
+    num_bins = n // 2 + 1
+    norm_zn_rad = np.zeros(num_bins, dtype=np.float64)
+    _dll.calc_norm_zn_rad_c(n, rho, norm_zn_rad)
+    return norm_zn_rad
+
 
 def rotate_angle(uvw0, mu, phi=None):
     """ Rotates direction cosines through a polar angle whose cosine is

From a52f31c95b65f4e97b4ee56ea1b1e5ea51991860 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 10 Aug 2018 16:40:17 -0500
Subject: [PATCH 02/18] Add cpp math functions to calc normalization arrays Add
 math funcs to reconstruct FET

---
 src/math_functions.cpp | 52 +++++++++++++++++++++++++++++++++++++++---
 src/math_functions.h   | 46 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 94 insertions(+), 4 deletions(-)

diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 8a36f7b922..b5c718d5d0 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -103,17 +103,36 @@ void calc_pn_c(int n, double x, double pnx[]) {
 }
 
 
-double evaluate_legendre_c(int n, const double data[], double x) {
+double evaluate_legendre_c(int n, const double data[], double x, int flag) {
   double pnx[n + 1];
   double val = 0.0;
   calc_pn_c(n, x, pnx);
-  for (int l = 0; l <= n; l++) {
-    val += (l + 0.5) * data[l] * pnx[l];
+
+  if (flag == 1) {
+    for (int l = 0; l <= n; l++) {
+      val += (l + 0.5) * data[l] * pnx[l];
+    }
+  } else if (flag == 2) {
+    for (int l = 0; l <= n; l++) {
+      val += data[l] * pnx[l];
+    }
   }
+
   return val;
 }
 
 
+void calc_norm_pn_c(int l, double zmin, double zmax, double norm_vec[]) {
+  // ===========================================================================
+  // Set up the normalization vector for FET Pn reconstruction
+  double d = zmax - zmin;
+
+  for (int i = 0; i <= l; i++) {
+    norm_vec[i] = (2 * i + 1)/d;
+  }
+}
+
+
 void calc_rn_c(int n, const double uvw[3], double rn[]){
   // rn[] is assumed to have already been allocated to the correct size
 
@@ -615,6 +634,33 @@ void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
 }
 
 
+void calc_norm_zn_rad_c(int n, double r, double norm_vec[]) {
+  // ===========================================================================
+  // Set up the normalization vector for FET Zn Radial only reconstruction
+  int index;
+
+  for (int i = 0; i <= n; i+=2) {
+    index = int(i/2);
+    norm_vec[index] = (i + 1)/(PI * r * r);
+  }
+}
+
+
+double evaluate_zernike_rad_c(int n, const double data[], double r) {
+  double rnr[n / 2 + 1];
+  double val = 0.0;
+  calc_zn_rad_c(n, r, rnr);
+  int index;
+
+  for (int i = 0; i <= n; i+ = 2) {
+    index = int(i/2);
+    val += data[index] * rnr[index];
+  }
+
+  return val;
+}
+
+
 void rotate_angle_c(double uvw[3], double mu, double* phi) {
   // Copy original directional cosines
   double u0 = uvw[0]; // original cosine in x direction
diff --git a/src/math_functions.h b/src/math_functions.h
index 7ef39c4461..30e210e500 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -53,11 +53,27 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
 //! @param data The polynomial expansion coefficient data; without the (2l+1)/2
 //!   factor.
 //! @param x    The value to evaluate at; x is expected to be within [-1,1]
+//! @param flag A flag of which method to calculate Legendre polynomials. 1 for
+//!   default (2l+1)/2; 2 for custom normalized coefficients
 //! @return The requested Legendre polynomials of order 0 to n (inclusive)
 //!   evaluated at x
 //==============================================================================
 
-extern "C" double evaluate_legendre_c(int n, const double data[], double x);
+extern "C" double evaluate_legendre_c(int n, const double data[], double x, int flag = 1);
+
+//==============================================================================
+//! Evaluate the normalization vector when reconstructing functional expansion
+//! from Legendre polynomials. 
+//!
+//! The normalization is (2 * l + 1)/(zmax - zmin)
+//!
+//! @param l         The maximum order requested
+//! @param zmin      The minimum location of 1-d domain
+//! @param zmax      The maximum location of 1-d domain
+//! @param norm_vec  The requested normalization of order 0 to n (inclusive)
+//==============================================================================
+
+extern "C" void calc_norm_pn_c(int l, double zmin, double zmax, double norm_vec[]);
 
 //==============================================================================
 //! Calculate the n-th order real spherical harmonics for a given angle (in
@@ -69,6 +85,7 @@ extern "C" double evaluate_legendre_c(int n, const double data[], double x);
 //!   evaluated at uvw.
 //==============================================================================
 
+
 extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
 
 //==============================================================================
@@ -109,6 +126,33 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 
 extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
 
+//==============================================================================
+//! Evaluate the normalization vector when reconstructing functional expansion
+//! from radial only Zernike polynomials. 
+//!
+//! The normalization is (2 * l + 1)/(zmax - zmin)
+//!
+//! @param n         The maximum order requested
+//! @param r      	 The radius of the disk
+//! @param norm_vec  The requested normalization of even orders from 0 to n
+//==============================================================================
+
+extern "C" void calc_norm_zn_rad_c(int n, double r, double norm_vec[]);
+
+//==============================================================================
+//! Find the value of f(x) given a set of even order Zernike coefficients and 
+//! the value of x if there is only radial dependency.
+//!
+//! @param n    The maximum order of the expansion
+//! @param data The overall polynomial expansion coefficient data;
+//!   normalization should be considered by users
+//! @param r    The value to evaluate at; r is expected to be within [0,1]
+//! @return The requested Zernike polynomials of order 0 to n (inclusive)
+//!   evaluated at r
+//==============================================================================
+
+extern "C" double evaluate_zernike_rad_c(int n, const double data[], double r);
+
 //==============================================================================
 //! Rotate the direction cosines through a polar angle whose cosine is mu and
 //! through an azimuthal angle sampled uniformly.

From 61431be559fbceb8a5b3fa809056a52702f90eb0 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 10 Aug 2018 16:41:24 -0500
Subject: [PATCH 03/18] Add unit test. Still need to add another tests for
 eva_zn

---
 tests/unit_tests/test_math.py | 49 ++++++++++++++++++++++++++++++++++-
 1 file changed, 48 insertions(+), 1 deletion(-)

diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 6825f6b930..a63561854d 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -57,9 +57,30 @@ def test_evaluate_legendre():
     test_vals = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x)
                           for x in test_xs])
 
+    # When turn on the flag 2, evaluate_legendre should be the same as
+    # np.polynomial.legendre.legval().
+    test_coeffs = [(2. * l + 1.) / 10 for l in range(max_order + 1)]
+    ref_vals_1 = np.polynomial.legendre.legval(test_xs, test_coeffs)
+
+    test_vals_1 = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x,
+                            2)for x in test_xs])
+
+    assert np.allclose(ref_vals, test_vals)
+    assert np.allclose(ref_vals_1, test_vals_1)
+
+
+def test_calc_norm_pn():
+    max_order = 10
+    zmin = -10
+    zmax = 10
+
+    ref_vals = np.divide(2 * np.arange(max_order + 1) + 1, (zmax-zmin))
+
+    test_vals = openmc.capi.math.calc_norm_pn(max_order, zmin, zmax)
+
     assert np.allclose(ref_vals, test_vals)
 
-
+    
 def test_calc_rn():
     max_order = 10
     test_ns = np.array([i for i in range(0, max_order + 1)])
@@ -137,6 +158,32 @@ def test_calc_zn():
     assert np.allclose(ref_vals, test_vals)
 
 
+def test_calc_zn_rad():
+    n = 10
+    rho = 0.5
+    phi = 0.5
+
+    # Reference solution from running the Fortran implementation
+    ref_vals = np.array([
+        1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
+        4.37500000e-01, -2.89062500e-01,-8.98437500e-02])
+
+    test_vals = openmc.capi.math.calc_zn_rad(n, rho)
+
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_calc_norm_zn_rad():
+    max_order = 10
+    rho = 0.5
+
+    ref_vals = np.divide(np.arange(0, max_order + 1, 2) + 1, np.pi * rho * rho)
+
+    test_vals = openmc.capi.math.calc_norm_zn_rad(max_order, rho)
+
+    assert np.allclose(ref_vals, test_vals)
+
+
 def test_rotate_angle():
     uvw0 = np.array([1., 0., 0.])
     phi = 0.

From 7e54eb0e976bd8a1d1e1f0c75638f3dbde555ffc Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Mon, 13 Aug 2018 10:21:22 -0500
Subject: [PATCH 04/18] Add some comments

---
 src/math_functions.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index b5c718d5d0..97dd7b2859 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -112,7 +112,7 @@ double evaluate_legendre_c(int n, const double data[], double x, int flag) {
     for (int l = 0; l <= n; l++) {
       val += (l + 0.5) * data[l] * pnx[l];
     }
-  } else if (flag == 2) {
+  } else if (flag == 2) { // if flag is 2, we don't apply any normalization
     for (int l = 0; l <= n; l++) {
       val += data[l] * pnx[l];
     }
@@ -647,6 +647,7 @@ void calc_norm_zn_rad_c(int n, double r, double norm_vec[]) {
 
 
 double evaluate_zernike_rad_c(int n, const double data[], double r) {
+  // data here is already normalized outside this function
   double rnr[n / 2 + 1];
   double val = 0.0;
   calc_zn_rad_c(n, r, rnr);

From 1b127cd2e8baca88686948fa8359c74c78513d6a Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Mon, 13 Aug 2018 10:35:16 -0500
Subject: [PATCH 05/18] Add test for evaluate zn_rad

---
 src/math_functions.cpp        |  2 +-
 tests/unit_tests/test_math.py | 17 ++++++++++++++++-
 2 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 97dd7b2859..72fb85ae87 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -653,7 +653,7 @@ double evaluate_zernike_rad_c(int n, const double data[], double r) {
   calc_zn_rad_c(n, r, rnr);
   int index;
 
-  for (int i = 0; i <= n; i+ = 2) {
+  for (int i = 0; i <= n; i+= 2) {
     index = int(i/2);
     val += data[index] * rnr[index];
   }
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index a63561854d..65fcd728f5 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -161,7 +161,6 @@ def test_calc_zn():
 def test_calc_zn_rad():
     n = 10
     rho = 0.5
-    phi = 0.5
 
     # Reference solution from running the Fortran implementation
     ref_vals = np.array([
@@ -172,6 +171,22 @@ def test_calc_zn_rad():
 
     assert np.allclose(ref_vals, test_vals)
 
+def test_evaluate_legendre():
+    max_order = 10
+    rho = 0.5
+
+    raw_zn = np.array([
+        1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
+        4.37500000e-01, -2.89062500e-01,-8.98437500e-02])
+    raw_coeff = np.ones(6)
+    norm_vec = np.divide(np.arange(0, max_order + 1, 2) + 1, np.pi * rho * rho)
+    real_coeff = np.multiply(raw_coeff,norm_vec)
+
+    ref_vals = np.sum(np.multiply(real_coeff,raw_zn))
+
+    test_vals = openmc.capi.math.evaluate_zernike_rad(real_coeff, rho)
+
+    assert np.allclose(ref_vals, test_vals)
 
 def test_calc_norm_zn_rad():
     max_order = 10

From 57387bf1361048b02fc35f8b8c91e28a87744f6f Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 10:39:25 -0500
Subject: [PATCH 06/18] Add Polynomial Class for Pn and Zn rad

---
 openmc/__init__.py    |  1 +
 openmc/polynomials.py | 35 +++++++++++++++++++++++++++++++++++
 2 files changed, 36 insertions(+)
 create mode 100644 openmc/polynomials.py

diff --git a/openmc/__init__.py b/openmc/__init__.py
index 191528327c..c7db3171c1 100644
--- a/openmc/__init__.py
+++ b/openmc/__init__.py
@@ -28,6 +28,7 @@ from openmc.particle_restart import *
 from openmc.mixin import *
 from openmc.plotter import *
 from openmc.search import *
+from openmc.polynomials import *
 from . import examples
 
 # Import a few convencience functions that used to be here
diff --git a/openmc/polynomials.py b/openmc/polynomials.py
new file mode 100644
index 0000000000..b3de1b9925
--- /dev/null
+++ b/openmc/polynomials.py
@@ -0,0 +1,35 @@
+import numpy as np
+import openmc
+import openmc.capi as capi
+
+class Polynomials(object):
+    """Class for Polynomials"""
+    def __init__(self, coeff):
+        self.coeff = coeff
+        self.n_coeff = len(coeff)
+
+class Legendre(Polynomials):
+    def __init__(self, coeff, domain_min = -1, domain_max = 1):
+        super().__init__(coeff)
+        self.order = self.n_coeff - 1
+        self.domain_min = domain_min
+        self.domain_max = domain_max
+        self.norm_vec = (2 * np.arange(self.n_coeff) + 1) / (self.domain_max
+            - self.domain_min)
+        self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
+
+    def __call__(self, x):
+        pn = capi.calc_pn(self.order, x)
+        return np.sum(np.multiply(self.norm_coeff, pn))
+
+class ZernikeRadial(Polynomials):
+    def __init__(self,coeff,radius = 1):
+        super().__init__(coeff)
+        self.order = 2 * (self.n_coeff - 1)
+        self.radius = radius
+        self.norm_vec = (2 * np.arange(self.n_coeff) + 1)/(np.pi * radius ** 2)
+        self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
+        
+    def __call__(self, r):
+        zn_rad = capi.calc_zn_rad(self.order, r)
+        return np.sum(np.multiply(self.norm_coeff, zn_rad))
\ No newline at end of file

From a67aab9c4091ccca24425c79c79c1ccb9c3c2f4a Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 10:46:34 -0500
Subject: [PATCH 07/18] revert math funcs for reconstruction, since it is not
 very necessary

---
 src/math_functions.cpp | 53 +++----------------------------------
 src/math_functions.h   | 59 +++++++-----------------------------------
 2 files changed, 13 insertions(+), 99 deletions(-)

diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 72fb85ae87..8a36f7b922 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -103,36 +103,17 @@ void calc_pn_c(int n, double x, double pnx[]) {
 }
 
 
-double evaluate_legendre_c(int n, const double data[], double x, int flag) {
+double evaluate_legendre_c(int n, const double data[], double x) {
   double pnx[n + 1];
   double val = 0.0;
   calc_pn_c(n, x, pnx);
-
-  if (flag == 1) {
-    for (int l = 0; l <= n; l++) {
-      val += (l + 0.5) * data[l] * pnx[l];
-    }
-  } else if (flag == 2) { // if flag is 2, we don't apply any normalization
-    for (int l = 0; l <= n; l++) {
-      val += data[l] * pnx[l];
-    }
+  for (int l = 0; l <= n; l++) {
+    val += (l + 0.5) * data[l] * pnx[l];
   }
-
   return val;
 }
 
 
-void calc_norm_pn_c(int l, double zmin, double zmax, double norm_vec[]) {
-  // ===========================================================================
-  // Set up the normalization vector for FET Pn reconstruction
-  double d = zmax - zmin;
-
-  for (int i = 0; i <= l; i++) {
-    norm_vec[i] = (2 * i + 1)/d;
-  }
-}
-
-
 void calc_rn_c(int n, const double uvw[3], double rn[]){
   // rn[] is assumed to have already been allocated to the correct size
 
@@ -634,34 +615,6 @@ void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
 }
 
 
-void calc_norm_zn_rad_c(int n, double r, double norm_vec[]) {
-  // ===========================================================================
-  // Set up the normalization vector for FET Zn Radial only reconstruction
-  int index;
-
-  for (int i = 0; i <= n; i+=2) {
-    index = int(i/2);
-    norm_vec[index] = (i + 1)/(PI * r * r);
-  }
-}
-
-
-double evaluate_zernike_rad_c(int n, const double data[], double r) {
-  // data here is already normalized outside this function
-  double rnr[n / 2 + 1];
-  double val = 0.0;
-  calc_zn_rad_c(n, r, rnr);
-  int index;
-
-  for (int i = 0; i <= n; i+= 2) {
-    index = int(i/2);
-    val += data[index] * rnr[index];
-  }
-
-  return val;
-}
-
-
 void rotate_angle_c(double uvw[3], double mu, double* phi) {
   // Copy original directional cosines
   double u0 = uvw[0]; // original cosine in x direction
diff --git a/src/math_functions.h b/src/math_functions.h
index 30e210e500..e0ecad60ec 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -53,27 +53,11 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
 //! @param data The polynomial expansion coefficient data; without the (2l+1)/2
 //!   factor.
 //! @param x    The value to evaluate at; x is expected to be within [-1,1]
-//! @param flag A flag of which method to calculate Legendre polynomials. 1 for
-//!   default (2l+1)/2; 2 for custom normalized coefficients
 //! @return The requested Legendre polynomials of order 0 to n (inclusive)
 //!   evaluated at x
 //==============================================================================
 
-extern "C" double evaluate_legendre_c(int n, const double data[], double x, int flag = 1);
-
-//==============================================================================
-//! Evaluate the normalization vector when reconstructing functional expansion
-//! from Legendre polynomials. 
-//!
-//! The normalization is (2 * l + 1)/(zmax - zmin)
-//!
-//! @param l         The maximum order requested
-//! @param zmin      The minimum location of 1-d domain
-//! @param zmax      The maximum location of 1-d domain
-//! @param norm_vec  The requested normalization of order 0 to n (inclusive)
-//==============================================================================
-
-extern "C" void calc_norm_pn_c(int l, double zmin, double zmax, double norm_vec[]);
+extern "C" double evaluate_legendre_c(int n, const double data[], double x);
 
 //==============================================================================
 //! Calculate the n-th order real spherical harmonics for a given angle (in
@@ -85,7 +69,6 @@ extern "C" void calc_norm_pn_c(int l, double zmin, double zmax, double norm_vec[
 //!   evaluated at uvw.
 //==============================================================================
 
-
 extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
 
 //==============================================================================
@@ -109,13 +92,18 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
 extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 
 //==============================================================================
-//! Calculate only the even order components of n-th order modified Zernike 
-//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho)
-//! location on the unit disk.
+//! Calculate only the even radial components of n-th order modified Zernike 
+//! polynomial moment with azimuthal dependency m = 0 for a given angle 
+//! (rho, theta) location on the unit disk.
 //!
 //! Since m = 0, n could only be even orders. Z_q0 = R_q0
 //!
-//! See calc_zn_c for methodology.
+//! This procedure uses the modified Kintner's method for calculating Zernike
+//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
+//! R. (2003). A comparative analysis of algorithms for fast computation of
+//! Zernike moments. Pattern Recognition, 36(3), 731-742.
+//! The normalization of the polynomials is such that the integral of Z_pq^2
+//! over the unit disk is exactly pi.
 //!
 //! @param n       The maximum order requested
 //! @param rho     The radial parameter to specify location on the unit disk
@@ -126,33 +114,6 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 
 extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
 
-//==============================================================================
-//! Evaluate the normalization vector when reconstructing functional expansion
-//! from radial only Zernike polynomials. 
-//!
-//! The normalization is (2 * l + 1)/(zmax - zmin)
-//!
-//! @param n         The maximum order requested
-//! @param r      	 The radius of the disk
-//! @param norm_vec  The requested normalization of even orders from 0 to n
-//==============================================================================
-
-extern "C" void calc_norm_zn_rad_c(int n, double r, double norm_vec[]);
-
-//==============================================================================
-//! Find the value of f(x) given a set of even order Zernike coefficients and 
-//! the value of x if there is only radial dependency.
-//!
-//! @param n    The maximum order of the expansion
-//! @param data The overall polynomial expansion coefficient data;
-//!   normalization should be considered by users
-//! @param r    The value to evaluate at; r is expected to be within [0,1]
-//! @return The requested Zernike polynomials of order 0 to n (inclusive)
-//!   evaluated at r
-//==============================================================================
-
-extern "C" double evaluate_zernike_rad_c(int n, const double data[], double r);
-
 //==============================================================================
 //! Rotate the direction cosines through a polar angle whose cosine is mu and
 //! through an azimuthal angle sampled uniformly.

From 1e7b5bb3274cb6e2863a677b2a3277adac708fc4 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 10:57:32 -0500
Subject: [PATCH 08/18] revert capi bc we will have new class

---
 openmc/capi/math.py | 93 ++-------------------------------------------
 1 file changed, 4 insertions(+), 89 deletions(-)

diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 697ed558d7..390ee466ff 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -13,10 +13,7 @@ _dll.calc_pn_c.restype = None
 _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
 _dll.evaluate_legendre_c.restype = c_double
-_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double, c_int]
-
-_dll.calc_norm_pn_c.restype = None
-_dll.calc_norm_pn_c.argtypes =[c_int, c_double, c_double, ndpointer(c_double)]
+_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
 
 _dll.calc_rn_c.restype = None
 _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
@@ -27,12 +24,6 @@ _dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
 _dll.calc_zn_rad_c.restype = None
 _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
-_dll.evaluate_zernike_rad_c.restype = c_double
-_dll.evaluate_zernike_rad_c.argtypes = [c_int, POINTER(c_double), c_double]
-
-_dll.calc_norm_zn_rad_c.restype = None
-_dll.calc_norm_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
-
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
                                 POINTER(c_double)]
@@ -90,7 +81,7 @@ def calc_pn(n, x):
     return pnx
 
 
-def evaluate_legendre(data, x, flag=1):
+def evaluate_legendre(data, x):
     """ Finds the value of f(x) given a set of Legendre coefficients
     and the value of x.
 
@@ -100,8 +91,6 @@ def evaluate_legendre(data, x, flag=1):
         Legendre coefficients
     x : float
         Independent variable to evaluate the Legendre at
-    flag : int
-        Flags to indicate use default normalization or not
 
     Returns
     -------
@@ -113,31 +102,7 @@ def evaluate_legendre(data, x, flag=1):
     data_arr = np.array(data, dtype=np.float64)
     return _dll.evaluate_legendre_c(len(data),
                                     data_arr.ctypes.data_as(POINTER(c_double)),
-                                    x, flag)
-
-
-def calc_norm_pn(l, zmin, zmax):
-    """ Calculate normalization arrays for n orders in Legendre polynomial.
-
-    Parameters
-    ----------
-    l : int
-        Maximum order
-    zmin : float
-        The lowest boundary of the domain
-    zmax : float
-        The highest boundary of the domain
-
-    Returns
-    -------
-    numpy.ndarray
-        Corresponding resulting list of normalization coefficients
-
-    """
-
-    norm_pn = np.empty(l + 1, dtype=np.float64)
-    _dll.calc_norm_pn_c(l, zmin, zmax, norm_pn)
-    return norm_pn
+                                    x)
 
 
 def calc_rn(n, uvw):
@@ -217,57 +182,7 @@ def calc_zn_rad(n, rho):
     zn_rad = np.zeros(num_bins, dtype=np.float64)
     _dll.calc_zn_rad_c(n, rho, zn_rad)
     return zn_rad
-
-
-def evaluate_zernike_rad(data, r):
-    """ Finds the value of f(x) given a set of even order Zernike coefficients
-    and the value of r.
-
-    Parameters
-    ----------
-    data : iterable of float
-        Legendre coefficients
-    r : float
-        Independent variable to evaluate the radial part of Zernike at
-
-    Returns
-    -------
-    float
-        Corresponding Zernike expansion result
-
-    """
-
-    data_arr = np.array(data, dtype=np.float64)
-    n = int(2 * (len(data) - 1))
-    return _dll.evaluate_zernike_rad_c(n,
-                                        data_arr.ctypes.data_as(POINTER(
-                                        c_double)),
-                                        r)
-
-
-def calc_norm_zn_rad(n, rho):
-    """ Calculate normalization arrays for the even orders in modified Zernike
-    polynomial moment with no azimuthal dependency (m=0).
-
-    Parameters
-    ----------
-    n : int
-        Maximum order
-    rho : float
-        Radial location in the unit disk
-
-    Returns
-    -------
-    numpy.ndarray
-        Corresponding resulting list of normalization coefficients
-
-    """
-
-    num_bins = n // 2 + 1
-    norm_zn_rad = np.zeros(num_bins, dtype=np.float64)
-    _dll.calc_norm_zn_rad_c(n, rho, norm_zn_rad)
-    return norm_zn_rad
-
+    
 
 def rotate_angle(uvw0, mu, phi=None):
     """ Rotates direction cosines through a polar angle whose cosine is

From 847e2cc28075e367c04c2beb8e54e717cdf0a13a Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 14:13:03 -0500
Subject: [PATCH 09/18] Revert test_math.py

---
 tests/unit_tests/test_math.py | 50 +----------------------------------
 1 file changed, 1 insertion(+), 49 deletions(-)

diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 65fcd728f5..7ef749dece 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -57,30 +57,9 @@ def test_evaluate_legendre():
     test_vals = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x)
                           for x in test_xs])
 
-    # When turn on the flag 2, evaluate_legendre should be the same as
-    # np.polynomial.legendre.legval().
-    test_coeffs = [(2. * l + 1.) / 10 for l in range(max_order + 1)]
-    ref_vals_1 = np.polynomial.legendre.legval(test_xs, test_coeffs)
-
-    test_vals_1 = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x,
-                            2)for x in test_xs])
-
-    assert np.allclose(ref_vals, test_vals)
-    assert np.allclose(ref_vals_1, test_vals_1)
-
-
-def test_calc_norm_pn():
-    max_order = 10
-    zmin = -10
-    zmax = 10
-
-    ref_vals = np.divide(2 * np.arange(max_order + 1) + 1, (zmax-zmin))
-
-    test_vals = openmc.capi.math.calc_norm_pn(max_order, zmin, zmax)
-
     assert np.allclose(ref_vals, test_vals)
 
-    
+
 def test_calc_rn():
     max_order = 10
     test_ns = np.array([i for i in range(0, max_order + 1)])
@@ -171,33 +150,6 @@ def test_calc_zn_rad():
 
     assert np.allclose(ref_vals, test_vals)
 
-def test_evaluate_legendre():
-    max_order = 10
-    rho = 0.5
-
-    raw_zn = np.array([
-        1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
-        4.37500000e-01, -2.89062500e-01,-8.98437500e-02])
-    raw_coeff = np.ones(6)
-    norm_vec = np.divide(np.arange(0, max_order + 1, 2) + 1, np.pi * rho * rho)
-    real_coeff = np.multiply(raw_coeff,norm_vec)
-
-    ref_vals = np.sum(np.multiply(real_coeff,raw_zn))
-
-    test_vals = openmc.capi.math.evaluate_zernike_rad(real_coeff, rho)
-
-    assert np.allclose(ref_vals, test_vals)
-
-def test_calc_norm_zn_rad():
-    max_order = 10
-    rho = 0.5
-
-    ref_vals = np.divide(np.arange(0, max_order + 1, 2) + 1, np.pi * rho * rho)
-
-    test_vals = openmc.capi.math.calc_norm_zn_rad(max_order, rho)
-
-    assert np.allclose(ref_vals, test_vals)
-
 
 def test_rotate_angle():
     uvw0 = np.array([1., 0., 0.])

From f7d098f0758231587b9d087c6d12ca0b6677dabd Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 14:14:05 -0500
Subject: [PATCH 10/18] Add some comments

---
 openmc/polynomials.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/openmc/polynomials.py b/openmc/polynomials.py
index b3de1b9925..6a9ba23c49 100644
--- a/openmc/polynomials.py
+++ b/openmc/polynomials.py
@@ -19,6 +19,7 @@ class Legendre(Polynomials):
         self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
 
     def __call__(self, x):
+        # x is the normalized position on [-1,1]
         pn = capi.calc_pn(self.order, x)
         return np.sum(np.multiply(self.norm_coeff, pn))
 
@@ -27,9 +28,11 @@ class ZernikeRadial(Polynomials):
         super().__init__(coeff)
         self.order = 2 * (self.n_coeff - 1)
         self.radius = radius
-        self.norm_vec = (2 * np.arange(self.n_coeff) + 1)/(np.pi * radius ** 2)
+        self.norm_vec = (2 * np.arange(self.n_coeff) + 1) / (np.pi * 
+            radius ** 2)
         self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
         
     def __call__(self, r):
+        # r is the normalized position on radius [0,1]
         zn_rad = capi.calc_zn_rad(self.order, r)
         return np.sum(np.multiply(self.norm_coeff, zn_rad))
\ No newline at end of file

From dc3064a79549c4f68d02e7c2b48cfaafcae54da2 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Thu, 16 Aug 2018 14:14:42 -0500
Subject: [PATCH 11/18] Test for the new Polynomial class

---
 tests/unit_tests/test_polynomials.py | 67 ++++++++++++++++++++++++++++
 1 file changed, 67 insertions(+)
 create mode 100644 tests/unit_tests/test_polynomials.py

diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py
new file mode 100644
index 0000000000..6dd220df36
--- /dev/null
+++ b/tests/unit_tests/test_polynomials.py
@@ -0,0 +1,67 @@
+import numpy as np
+import scipy as sp
+
+import openmc
+import openmc.capi
+
+def test_legendre():
+    coeff = np.asarray([2.5, 1.2, -0.2, 0.1])
+    pn = openmc.Legendre(coeff)
+    assert pn.order == 3
+    assert pn.domain_min == -1
+    assert pn.domain_max == 1
+
+    norm_vec = (2 * np.arange(4) + 1) / 2
+    assert np.allclose(pn.norm_vec,norm_vec)
+    norm_coeff = np.multiply(norm_vec,coeff)
+    assert np.allclose(pn.norm_coeff, norm_coeff)
+
+    coeff = np.asarray([2.5, 1.2, -0.2])
+    pn = openmc.Legendre(coeff, -10, 10)
+    assert pn.order == 2
+    assert pn.domain_min == -10
+    assert pn.domain_max == 10
+    norm_vec = (2 * np.arange(3) + 1) / 20
+    assert np.allclose(pn.norm_vec,norm_vec)
+    norm_coeff = np.multiply(norm_vec,coeff)
+    assert np.allclose(pn.norm_coeff, norm_coeff)
+
+    test_xs = np.linspace(-1., 1., num=5, endpoint=True)
+    ref_vals = np.polynomial.legendre.legval(test_xs, norm_coeff)
+
+    test_vals = np.array([pn(x) for x in test_xs])
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_zernike_radial():
+    coeff = np.asarray([1.3, -3.0, 9e-1, -6e-1, 0.11])
+    zn_rad = openmc.ZernikeRadial(coeff)
+    assert zn_rad.order == 8
+    assert zn_rad.radius == 1
+    norm_vec = (2 * np.arange(5) + 1) / np.pi
+    assert np.allclose(zn_rad.norm_vec,norm_vec)
+    norm_coeff = np.multiply(norm_vec,coeff)
+    assert np.allclose(zn_rad.norm_coeff,norm_coeff)
+
+    coeff = np.asarray([1.3, -3.0, 9e-1, -6e-1, 0.11, 0.222])
+    zn_rad = openmc.ZernikeRadial(coeff, 0.392)
+    assert zn_rad.order == 10
+    assert zn_rad.radius == 0.392
+    norm_vec = (2 * np.arange(6) + 1) / (np.pi * 0.392 ** 2)
+    assert np.allclose(zn_rad.norm_vec,norm_vec)
+    norm_coeff = np.multiply(norm_vec,coeff)
+    assert np.allclose(zn_rad.norm_coeff,norm_coeff)
+
+
+    rho = 0.5
+    # Reference solution from running the Fortran implementation   
+    raw_zn = np.array([
+        1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
+        4.37500000e-01, -2.89062500e-01,-8.98437500e-02])
+
+    ref_vals = np.sum(np.multiply(norm_coeff,raw_zn))
+
+    test_vals = zn_rad(rho)
+
+    assert ref_vals == test_vals
+

From 15be98de928f7a4a4a56c1a1a926e4fadb646e88 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 17 Aug 2018 14:59:37 -0500
Subject: [PATCH 12/18] Remove odd files

---
 openmc/polynomials.py | 38 --------------------------------------
 1 file changed, 38 deletions(-)
 delete mode 100644 openmc/polynomials.py

diff --git a/openmc/polynomials.py b/openmc/polynomials.py
deleted file mode 100644
index 6a9ba23c49..0000000000
--- a/openmc/polynomials.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import numpy as np
-import openmc
-import openmc.capi as capi
-
-class Polynomials(object):
-    """Class for Polynomials"""
-    def __init__(self, coeff):
-        self.coeff = coeff
-        self.n_coeff = len(coeff)
-
-class Legendre(Polynomials):
-    def __init__(self, coeff, domain_min = -1, domain_max = 1):
-        super().__init__(coeff)
-        self.order = self.n_coeff - 1
-        self.domain_min = domain_min
-        self.domain_max = domain_max
-        self.norm_vec = (2 * np.arange(self.n_coeff) + 1) / (self.domain_max
-            - self.domain_min)
-        self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
-
-    def __call__(self, x):
-        # x is the normalized position on [-1,1]
-        pn = capi.calc_pn(self.order, x)
-        return np.sum(np.multiply(self.norm_coeff, pn))
-
-class ZernikeRadial(Polynomials):
-    def __init__(self,coeff,radius = 1):
-        super().__init__(coeff)
-        self.order = 2 * (self.n_coeff - 1)
-        self.radius = radius
-        self.norm_vec = (2 * np.arange(self.n_coeff) + 1) / (np.pi * 
-            radius ** 2)
-        self.norm_coeff = np.multiply(self.norm_vec,self.coeff)
-        
-    def __call__(self, r):
-        # r is the normalized position on radius [0,1]
-        zn_rad = capi.calc_zn_rad(self.order, r)
-        return np.sum(np.multiply(self.norm_coeff, zn_rad))
\ No newline at end of file

From 98a32391be03095a2d9cc7c9cc89662275bbd8ee Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 17 Aug 2018 15:00:39 -0500
Subject: [PATCH 13/18] Change Polynomial class and fix PEP 8 style

---
 openmc/__init__.py   |  2 +-
 openmc/polynomial.py | 71 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 72 insertions(+), 1 deletion(-)
 create mode 100644 openmc/polynomial.py

diff --git a/openmc/__init__.py b/openmc/__init__.py
index c7db3171c1..4627c340b2 100644
--- a/openmc/__init__.py
+++ b/openmc/__init__.py
@@ -28,7 +28,7 @@ from openmc.particle_restart import *
 from openmc.mixin import *
 from openmc.plotter import *
 from openmc.search import *
-from openmc.polynomials import *
+from openmc.polynomial import *
 from . import examples
 
 # Import a few convencience functions that used to be here
diff --git a/openmc/polynomial.py b/openmc/polynomial.py
new file mode 100644
index 0000000000..fce068f188
--- /dev/null
+++ b/openmc/polynomial.py
@@ -0,0 +1,71 @@
+import numpy as np
+import openmc
+import openmc.capi as capi
+
+
+def legendre_from_expcoef(coef, domain):
+    """Download file from a URL
+
+    Parameters
+    ----------
+    coef : Iterable of float
+        A list of coefficients of each term in Legendre polynomials
+
+    Returns
+    -------
+    Legendre class
+        A numpy Legendre series object
+
+    """
+    n = np.arange(len(coef) + 1)
+    c = (2*n + 1) * np.asarray(coef) / (domain[1] - domain[0])
+    return np.polynomial.Legendre(c, domain)
+
+
+class Polynomial(object):
+    """Abstract Polynomial Class for creating polynomials.
+    """
+    def __init__(self, coef):
+        self.coef = np.asarray(coef)
+
+
+class ZernikeRadial(Polynomial):
+    """Create radial only Zernike polynomials given coefficients and domian.
+
+    The radial only Zernike polynomials are defined as in
+    :class:`ZernikeRadialFilter`.
+
+    Parameters
+    ----------
+    coef : Iterable of float
+        A list of coefficients of each term in radial only Zernike polynomials
+    radius : float
+        Domain of Zernike polynomials to be applied on. Default is 1.
+    r : float
+        Position to be evaluated, normalized on radius [0,1]
+
+    Attributes
+    ----------
+    order : int
+        The maximum (even) order of Zernike polynomials.
+    radius : float
+        Domain of Zernike polynomials to be applied on. Default is 1.
+    norm_coef : iterable of float
+        The list of coefficients of each term in the polynomials after
+        normailization.
+
+    """
+    def __init__(self, coef, radius=1):
+        super().__init__(coef)
+        self._order = 2 * (self.n_coef - 1)
+        self.radius = radius
+        norm_vec = (2 * np.arange(self.n_coef) + 1) / (np.pi * radius**2)
+        self._norm_coef = norm_vec * self.coef
+
+    @property
+    def order(self):
+        return self._order
+
+    def __call__(self, r):
+        zn_rad = capi.calc_zn_rad(self.order, r)
+        return np.sum(self._norm_coef * zn_rad)

From 33be934308ba3c8ad1f9e3c40143d8ed3992142e Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 17 Aug 2018 15:27:51 -0500
Subject: [PATCH 14/18] Add test of Zernike radial.

There is really nothing to test for Legendre
---
 openmc/polynomial.py                 |  8 +++--
 tests/unit_tests/test_polynomials.py | 44 +++-------------------------
 2 files changed, 9 insertions(+), 43 deletions(-)

diff --git a/openmc/polynomial.py b/openmc/polynomial.py
index fce068f188..7264524cbc 100644
--- a/openmc/polynomial.py
+++ b/openmc/polynomial.py
@@ -4,12 +4,14 @@ import openmc.capi as capi
 
 
 def legendre_from_expcoef(coef, domain):
-    """Download file from a URL
+    """Return a Legendre object
 
     Parameters
     ----------
     coef : Iterable of float
         A list of coefficients of each term in Legendre polynomials
+    domain : (2,) List of float
+        Domain of the Legendre polynomial
 
     Returns
     -------
@@ -57,9 +59,9 @@ class ZernikeRadial(Polynomial):
     """
     def __init__(self, coef, radius=1):
         super().__init__(coef)
-        self._order = 2 * (self.n_coef - 1)
+        self._order = 2 * (len(self.coef) - 1)
         self.radius = radius
-        norm_vec = (2 * np.arange(self.n_coef) + 1) / (np.pi * radius**2)
+        norm_vec = (2 * np.arange(len(self.coef)) + 1) / (np.pi * radius**2)
         self._norm_coef = norm_vec * self.coef
 
     @property
diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py
index 6dd220df36..4ce767d2ea 100644
--- a/tests/unit_tests/test_polynomials.py
+++ b/tests/unit_tests/test_polynomials.py
@@ -4,64 +4,28 @@ import scipy as sp
 import openmc
 import openmc.capi
 
-def test_legendre():
-    coeff = np.asarray([2.5, 1.2, -0.2, 0.1])
-    pn = openmc.Legendre(coeff)
-    assert pn.order == 3
-    assert pn.domain_min == -1
-    assert pn.domain_max == 1
-
-    norm_vec = (2 * np.arange(4) + 1) / 2
-    assert np.allclose(pn.norm_vec,norm_vec)
-    norm_coeff = np.multiply(norm_vec,coeff)
-    assert np.allclose(pn.norm_coeff, norm_coeff)
-
-    coeff = np.asarray([2.5, 1.2, -0.2])
-    pn = openmc.Legendre(coeff, -10, 10)
-    assert pn.order == 2
-    assert pn.domain_min == -10
-    assert pn.domain_max == 10
-    norm_vec = (2 * np.arange(3) + 1) / 20
-    assert np.allclose(pn.norm_vec,norm_vec)
-    norm_coeff = np.multiply(norm_vec,coeff)
-    assert np.allclose(pn.norm_coeff, norm_coeff)
-
-    test_xs = np.linspace(-1., 1., num=5, endpoint=True)
-    ref_vals = np.polynomial.legendre.legval(test_xs, norm_coeff)
-
-    test_vals = np.array([pn(x) for x in test_xs])
-    assert np.allclose(ref_vals, test_vals)
-
 
 def test_zernike_radial():
     coeff = np.asarray([1.3, -3.0, 9e-1, -6e-1, 0.11])
     zn_rad = openmc.ZernikeRadial(coeff)
     assert zn_rad.order == 8
     assert zn_rad.radius == 1
-    norm_vec = (2 * np.arange(5) + 1) / np.pi
-    assert np.allclose(zn_rad.norm_vec,norm_vec)
-    norm_coeff = np.multiply(norm_vec,coeff)
-    assert np.allclose(zn_rad.norm_coeff,norm_coeff)
 
     coeff = np.asarray([1.3, -3.0, 9e-1, -6e-1, 0.11, 0.222])
     zn_rad = openmc.ZernikeRadial(coeff, 0.392)
     assert zn_rad.order == 10
     assert zn_rad.radius == 0.392
     norm_vec = (2 * np.arange(6) + 1) / (np.pi * 0.392 ** 2)
-    assert np.allclose(zn_rad.norm_vec,norm_vec)
-    norm_coeff = np.multiply(norm_vec,coeff)
-    assert np.allclose(zn_rad.norm_coeff,norm_coeff)
-
+    norm_coeff = norm_vec * coeff
 
     rho = 0.5
-    # Reference solution from running the Fortran implementation   
+    # Reference solution from running the Fortran implementation
     raw_zn = np.array([
         1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
-        4.37500000e-01, -2.89062500e-01,-8.98437500e-02])
+        4.37500000e-01, -2.89062500e-01, -8.98437500e-02])
 
-    ref_vals = np.sum(np.multiply(norm_coeff,raw_zn))
+    ref_vals = np.sum(norm_coeff * raw_zn)
 
     test_vals = zn_rad(rho)
 
     assert ref_vals == test_vals
-

From cc295001ff2812b7b4b36072c98be963be03adf9 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 17 Aug 2018 16:39:56 -0500
Subject: [PATCH 15/18] Fixed a minor bug

---
 openmc/polynomial.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/openmc/polynomial.py b/openmc/polynomial.py
index 7264524cbc..d26d7edf83 100644
--- a/openmc/polynomial.py
+++ b/openmc/polynomial.py
@@ -3,8 +3,11 @@ import openmc
 import openmc.capi as capi
 
 
-def legendre_from_expcoef(coef, domain):
-    """Return a Legendre object
+def legendre_from_expcoef(coef, domain= [-1,1]):
+    """Return a Legendre object.
+
+    Given a list of coefficients from FET tally and a array of down, return
+    the numpy Legendre object.
 
     Parameters
     ----------
@@ -19,7 +22,8 @@ def legendre_from_expcoef(coef, domain):
         A numpy Legendre series object
 
     """
-    n = np.arange(len(coef) + 1)
+
+    n = np.arange(len(coef))
     c = (2*n + 1) * np.asarray(coef) / (domain[1] - domain[0])
     return np.polynomial.Legendre(c, domain)
 

From 23f71ccfad0bf1825139ed3acc160491c1744747 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 17 Aug 2018 16:40:17 -0500
Subject: [PATCH 16/18] Add Poly in docs

---
 docs/source/pythonapi/base.rst | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 752cf417e8..497b8fc3dc 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -178,6 +178,24 @@ Post-processing
    openmc.StatePoint
    openmc.Summary
 
+Polynomial
+---------------
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclass.rst 
+
+   openmc.ZernikeRadial
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myfunction.rst
+
+   openmc.legendre_from_expcoef
+
+
 Various classes may be created when performing tally slicing and/or arithmetic:
 
 .. autosummary::

From 058ceda821fff5f14ae2c95fbbaaac04ab5cd9b8 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Tue, 21 Aug 2018 10:16:06 -0500
Subject: [PATCH 17/18] Change based on suggestion

---
 openmc/polynomial.py                 | 10 +++++-----
 tests/unit_tests/test_polynomials.py |  2 --
 2 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/openmc/polynomial.py b/openmc/polynomial.py
index d26d7edf83..b120659976 100644
--- a/openmc/polynomial.py
+++ b/openmc/polynomial.py
@@ -3,8 +3,8 @@ import openmc
 import openmc.capi as capi
 
 
-def legendre_from_expcoef(coef, domain= [-1,1]):
-    """Return a Legendre object.
+def legendre_from_expcoef(coef, domain= (-1,1)):
+    """Return a Legendre series object based on expansion coefficients.
 
     Given a list of coefficients from FET tally and a array of down, return
     the numpy Legendre object.
@@ -18,8 +18,8 @@ def legendre_from_expcoef(coef, domain= [-1,1]):
 
     Returns
     -------
-    Legendre class
-        A numpy Legendre series object
+    numpy.polynomial.Legendre
+        A numpy Legendre series class
 
     """
 
@@ -36,7 +36,7 @@ class Polynomial(object):
 
 
 class ZernikeRadial(Polynomial):
-    """Create radial only Zernike polynomials given coefficients and domian.
+    """Create radial only Zernike polynomials given coefficients and domain.
 
     The radial only Zernike polynomials are defined as in
     :class:`ZernikeRadialFilter`.
diff --git a/tests/unit_tests/test_polynomials.py b/tests/unit_tests/test_polynomials.py
index 4ce767d2ea..52c30bb000 100644
--- a/tests/unit_tests/test_polynomials.py
+++ b/tests/unit_tests/test_polynomials.py
@@ -1,8 +1,6 @@
 import numpy as np
-import scipy as sp
 
 import openmc
-import openmc.capi
 
 
 def test_zernike_radial():

From dcf6cfe0a1ae399ead325d947f3be74e9e263e45 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Tue, 21 Aug 2018 10:21:01 -0500
Subject: [PATCH 18/18] Change docs

---
 docs/source/pythonapi/base.rst | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 497b8fc3dc..84bb3cc54d 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -178,8 +178,7 @@ Post-processing
    openmc.StatePoint
    openmc.Summary
 
-Polynomial
----------------
+The following classes and functions are used for functional expansion reconstruction.
 
 .. autosummary::
    :toctree: generated