From cdfa3565c9f95f793206a194563a1948a48e1572 Mon Sep 17 00:00:00 2001
From: Zhuoran Han <zhuoranh@mit.edu>
Date: Fri, 3 Aug 2018 12:42:18 -0500
Subject: [PATCH] Add math function for radial Zernike

---
 src/math.F90           | 10 ++++++++++
 src/math_functions.cpp | 29 +++++++++++++++++++++++++++++
 src/math_functions.h   | 23 +++++++++++++++++++++++
 3 files changed, 62 insertions(+)

diff --git a/src/math.F90 b/src/math.F90
index 33f1ab4d30..87c747f579 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -11,6 +11,7 @@ module math
   public :: calc_pn
   public :: calc_rn
   public :: calc_zn
+  public :: calc_zn_rad
   public :: evaluate_legendre
   public :: rotate_angle
   public :: maxwell_spectrum
@@ -70,6 +71,15 @@ module math
       real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
     end subroutine calc_zn
 
+    pure subroutine calc_zn_rad(n, rho, phi, zn_rad) bind(C, name='calc_zn_rad_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), value, intent(in) :: rho
+      real(C_DOUBLE), value, intent(in) :: phi
+      real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1)
+    end subroutine calc_zn
+
     subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
       use ISO_C_BINDING
       implicit none
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 31f86ea5f6..2a79bb4b24 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -587,6 +587,35 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) {
 
 }
 
+void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) {
+  // Calculate R_p0(rho) as Zn_p0(rho)
+  // Set up the array of the coefficients
+  int length = int(n/2) + 1;
+  double zn_rad[length];
+  double q = 0;
+
+  // R_00 is always 1
+  zn_rad[0] = 1;
+
+  // Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong)
+  for (int p = 2; p <= n; p += 2) {
+    int index = int(p/2);
+    if (p == 2) {
+    // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong)
+      R_22 = std::pow(rho, 2);
+      zn_rad[index] = 2 * R_22 - zn_rad[0]
+    }
+    else {
+      double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;
+      double k2 = 2 * p * (p - 1) * (p - 2);
+      double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);
+      double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;
+      zn_rad[index] = 
+        ((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1;
+    }
+  }
+}
+
 
 void rotate_angle_c(double uvw[3], double mu, double* phi) {
   // Copy original directional cosines
diff --git a/src/math_functions.h b/src/math_functions.h
index 2ceea98d87..98882463df 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -91,6 +91,29 @@ 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 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
+//!
+//! 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
+//! @param phi     The angle parameter to specify location on the unit disk
+//! @param zn_rad  The requested moments of order 0 to n (inclusive)
+//!   evaluated at rho and phi when m = 0.
+//==============================================================================
+
+extern "C" void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]);
+
 //==============================================================================
 //! Rotate the direction cosines through a polar angle whose cosine is mu and
 //! through an azimuthal angle sampled uniformly.