# HG changeset patch
# User Ramom Flores <jose.ramom.flores.das.seixas@gmail.com>
# Date 1564942240 -7200
#      Sun Aug 04 20:10:40 2019 +0200
# Node ID 437f7f1fae96cb71384639cb179aeabe796505aa
# Parent  f97c278a5ac0b480bd44a41bd6b57ff5a846bfeb
Added two new functions and documentation

user: Ramom Flores <jose.ramom.flores.das.seixas@gmail.com>
branch 'default'
added doc/zernikes_and_derivatives_cartesian_OSA.pdf
added doc/zernikes_and_derivatives_cartesian_OSA.tex
added inst/zernike_osa_ansi_to_mn.m
added inst/zernikes_and_derivatives_cartesian_OSA.m

diff -r f97c278a5ac0 -r 437f7f1fae96 doc/zernikes_and_derivatives_cartesian_OSA.pdf
Binary file doc/zernikes_and_derivatives_cartesian_OSA.pdf has changed
diff -r f97c278a5ac0 -r 437f7f1fae96 doc/zernikes_and_derivatives_cartesian_OSA.tex
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/zernikes_and_derivatives_cartesian_OSA.tex	Sun Aug 04 20:10:40 2019 +0200
@@ -0,0 +1,356 @@
+%% LyX 2.3.3 created this file.  For more info, see http://www.lyx.org/.
+%% Do not edit unless you really know what you are doing.
+\documentclass[english]{article}
+\usepackage[T1]{fontenc}
+\usepackage[latin9]{inputenc}
+\usepackage[a4paper]{geometry}
+\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
+\usepackage{color}
+\usepackage{babel}
+\usepackage{array}
+\usepackage{amsmath}
+\usepackage[unicode=true,pdfusetitle,
+ bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
+ breaklinks=false,pdfborder={0 0 0},pdfborderstyle={},backref=false,colorlinks=true]
+ {hyperref}
+
+\makeatletter
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
+%% Because html converters don't know tabularnewline
+\providecommand{\tabularnewline}{\\}
+
+\makeatother
+
+\begin{document}
+\title{Calculation of the Zernike polynomials and their partial derivatives
+in Cartesian coordinates using the standard indexes OSA/ANSI.}
+\author{Jos� Ramom Flores das Seixas}
+\date{August 4, 2019}
+
+\maketitle
+\tableofcontents{}
+
+\section{Introduction}
+
+The \texttt{zernikes\_and\_derivatives\_cartesian\_OSA} function calculates
+the Zernike polynomials and their partial derivatives in Cartesian
+coordinates. For this purpose, it uses the algorithm described in
+reference \cite{Andersen_2018}. However, it should be noted that
+in the cited article the author use unit-normalized Zernike polynomials
+arranged according to the azimuthal scheme set forth by Rimmer and
+Wyant \cite{Rimmer=000026Wyant_1975}, while to implement the function
+has been used the OSA/ANSI standard notation described in references
+\cite{Thibos_et_al_2002,ANSI-2017}. This means the index scheme is
+different and, moreover, the polynomials are not unit-normalized but
+normalized to $\pi$.
+
+\section{Notation}
+
+Zernike's polynomials are usually ordered by a double-index, $Z_{n}^{m}$,
+being $n$ the radial order and $m$ the angular frequency, both integers.
+Although in programming it is usually utilized a single index, $Z_{j}$.
+
+The scheme used in reference \cite{Andersen_2018}, consider a double-index
+with $n\geq0$ and $0\leq m\leq n$. However, the standard OSA/ANSI\footnote{In fact, this double index schema is the most usual, and is employed
+by other notations than the standard OSA/ANSI.} uses a double-index with $n\geq0$ and $-n\leq m\leq n$, as shown
+in the following table:
+\noindent \begin{center}
+\begin{tabular}{|c|c|c|c|>{\centering}p{3cm}|}
+\multicolumn{2}{|c|}{Rimmer\&Wyant schema} &  & \multicolumn{2}{c|}{~~~~OSA standard~~~~}\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+$\tilde{n}$ & $\tilde{m}$ &  & $n$ & $m$\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+0 & 0 &  & 0 & 0\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+1 & 0, 1 &  & 1 & -1, 1\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+2 & 0, 1, 2 &  & 2 & -2, 0, 2\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+3 & 0, 1, 2, 3 &  & 3 & -3, -1, 1, 3\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+4 & 0, 1, 2, 3, 4 &  & 4 & -4, -2, 0, 2, 4\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+$\vdots$ & $\vdots$ &  & $\vdots$ & $\vdots$\tabularnewline
+\cline{1-2} \cline{2-2} \cline{4-5} \cline{5-5} 
+\end{tabular}
+\par\end{center}
+
+Therefore, by the way of example, the following polynomials --using
+the two schemes seen-- are equivalent:
+\noindent \begin{center}
+\begin{tabular}{ccc}
+Rimmer\&Wyant schema &  & ~~~~OSA standard~~~~\tabularnewline
+\cline{1-1} \cline{3-3} 
+$Z_{1}^{0}\left(x,y\right)$ & ~~~~ & $Z_{1}^{-1}\left(x,y\right)$\tabularnewline
+$Z_{3}^{0}\left(x,y\right)$ &  & $Z_{3}^{-3}\left(x,y\right)$\tabularnewline
+$Z_{3}^{1}\left(x,y\right)$ &  & $Z_{3}^{-1}\left(x,y\right)$\tabularnewline
+$Z_{4}^{4}\left(x,y\right)$ &  & $Z_{4}^{4}\left(x,y\right)$\tabularnewline
+$Z_{8}^{6}\left(x,y\right)$ &  & $Z_{8}^{4}\left(x,y\right)$\tabularnewline
+\end{tabular}
+\par\end{center}
+
+Using $\left(\tilde{n},\tilde{m}\right)$ to denote the R\&W indices,
+and $\left(n,m\right)$ those of the OSA scheme, both pairs of indices
+are related as follow:
+\begin{equation}
+n=\tilde{n}\qquad m=2\cdot\tilde{m}-\tilde{n}\label{eq:R=000026W_to_OSA}
+\end{equation}
+
+
+\subsection{Single index}
+
+As said before, it is more convenient to use a simple index to program.
+In the OSA standard, the conversion of indices is \href{https://en.wikipedia.org/wiki/Zernike_polynomials\#OSA/ANSI_standard_indices}{as follows}:
+\[
+j=\frac{n\left(n+2\right)+m}{2}\qquad\begin{array}{l}
+n=roundup\left[\dfrac{-3+\sqrt{9+8\,j}}{2}\right]\\
+\\
+m=2\,j-n\left(n+2\right)
+\end{array}
+\]
+
+
+\subsection{U-polynomials}
+
+Employing the notation used in reference \cite{Andersen_2018}, the
+unit-normalized Zernike's polynomials are called $U_{nm}$, and are
+related to the $\pi$-normalized ones as follows:
+\[
+Z_{n}^{m}\left(x,y\right)=N_{nm}\cdot U_{nm}\left(x,y\right)
+\]
+ with
+\[
+N_{nm}=\sqrt{\frac{2\left(n+1\right)}{1+\delta_{m0}}}\qquad\textrm{being}\quad\delta_{m0}=\left\{ \begin{array}{ll}
+1 & m=0\\
+0 & m\neq0
+\end{array}\right.
+\]
+
+
+\section{Algorithm}
+
+It is well known that for high-order Zernike's polynomials, computing
+the values using explicit expressions is not the best strategy, since
+it is inefficient and suffers form large cancellation errors. So,
+several schemes using recurrence relations have been devised.
+
+Reference \cite{Andersen_2018} presents one recurrence relation with
+coefficients that do not depend on radial or azimuthal orders and
+which contains no singularities. In addition, it also presents a recurrence
+relation to compute the partial derivatives in Cartesian coordinates.
+
+\subsection{Recurrence relations in Cartesian coordinates for Zernike U-polynomials}
+
+The general recurrence relation is
+\[
+U_{\tilde{n},\tilde{m}}=xU_{\tilde{n}-1,\tilde{m}}+yU_{\tilde{n}-1,\tilde{n}-1-\tilde{m}}+xU_{\tilde{n}-1,\tilde{m}-1}-yU_{\tilde{n}-1,\tilde{n}-\tilde{m}}-U_{\tilde{n}-2,\tilde{m}-1}
+\]
+But there are several exceptions:
+\begin{itemize}
+\item for $\tilde{m}=0$$\qquad U_{\tilde{n},0}=xU_{\tilde{n}-1,0}+yU_{\tilde{n}-1,\tilde{n}-1}$
+\item for $\tilde{m}=\tilde{n}$$\qquad U_{\tilde{n},\tilde{n}}=xU_{\tilde{n}-1,\tilde{n}-1}-yU_{\tilde{n}-1,0}$
+\item for $\tilde{n}$ odd and $\tilde{m}=\dfrac{\tilde{n}-1}{2}$
+\[
+U_{\tilde{n},\tilde{m}}=yU_{\tilde{n}-1,\tilde{n}-1-\tilde{m}}+xU_{\tilde{n}-1,\tilde{m}-1}-yU_{\tilde{n}-1,\tilde{n}-\tilde{m}}-U_{\tilde{n}-2,\tilde{m}-1}
+\]
+\item for $\tilde{n}$ odd and $\tilde{m}=\dfrac{\tilde{n}-1}{2}+1$
+\[
+U_{\tilde{n},\tilde{m}}=xU_{\tilde{n}-1,\tilde{m}}+yU_{\tilde{n}-1,\tilde{n}-1-\tilde{m}}+yU_{\tilde{n}-1,\tilde{m}-1}-U_{\tilde{n}-2,\tilde{m}-1}
+\]
+\item for $\tilde{n}$ even and $\tilde{m}=\dfrac{\tilde{n}}{2}$
+\[
+U_{\tilde{n},\tilde{m}}=2xU_{\tilde{n}-1,\tilde{m}}+2yU_{\tilde{n}-1,\tilde{m}-1}-U_{\tilde{n}-2,\tilde{m}-1}
+\]
+\end{itemize}
+The starting polynomials being
+\[
+U_{0,0}=1,\quad U_{1,0}=y,\quad U_{1,1}=x\quad.
+\]
+
+
+\subsection{Recurrence relations for the OSA/ANSI scheme}
+
+If we now ``translate'' the previous equations into the OSA scheme,
+we obtain that the general recurrence relation is
+\[
+U_{n.m}=xU_{n-1,m+1}+yU_{n-1,-\left(m+1\right)}+xU_{n-1,m-1}-yU_{n-1,1-m}-U_{n-2,m}
+\]
+the exceptions are given by
+\begin{itemize}
+\item for $m=-n$$\qquad U_{n,-n}=xU_{n-1,1-n}+yU_{n-1,n-1}$
+\item for $m=n$$\qquad U_{n,n}=xU_{n-1,n-1}-yU_{n-1,1-n}$
+\item for $m=-1$$\qquad U_{n,-1}=yU_{n-1,0}+xU_{n-1,-2}-yU_{n-1,2}-U_{n-2,-1}$
+\item for $m=1$$\qquad U_{n,1}=xU_{n-1,2}+yU_{n-1,-2}+xU_{n-1,0}-U_{n-2,1}$
+\item for $m=0$$\qquad U_{n,0}=2xU_{n-1,1}+2yU_{n-1,-1}-U_{n-2,0}$
+\end{itemize}
+and the starting polynomials are
+\[
+U_{0,0}=1,\quad U_{1,-1}=y\quad U_{1,1}=x\quad.
+\]
+
+Note that using the OSA scheme the conditions of the exceptions have
+been simplified, which makes it easier to compute.
+
+\subsubsection{Recurrence relations with a single index}
+
+As mentioned earlier, when programming it is more convenient to use
+a simple index. In this case, recurrence relations are given as follows.
+Let us start now with the exceptions:
+\begin{itemize}
+\item for $m=-n\quad\Rightarrow\quad j=\dfrac{\left(n+1\right)n}{n}$
+\[
+U_{j}=xU_{j-n}+yU_{j-1}
+\]
+\item for $m=n\quad\Rightarrow\quad j=\dfrac{n\left(n+3\right)}{2}$
+\[
+U_{j}=xU_{j-\left(n+1\right)}-yU_{j-2n}
+\]
+\item for $m=-1\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)-1}{2}$
+\[
+U_{j}=xU_{j-\left(n+1\right)}+yU_{j-n}-yU_{j-\left(n-1\right)}-U_{j-2n}
+\]
+\item for $m=1\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)+1}{2}$
+\[
+U_{j}=xU_{j-n}+xU_{j-\left(n+1\right)}+yU_{j-\left(n+2\right)}-U_{j-2n}
+\]
+\item for $m=0\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)}{2}$
+\[
+U_{j}=2xU_{j-n}+2yU_{j-\left(n+1\right)}-U_{j-2n}\,.
+\]
+\end{itemize}
+The general case is
+\[
+U_{j}=xU_{j-n}+yU_{j-\left(n+m+1\right)}+xU_{j-\left(n+1\right)}-yU_{j-\left(m+m\right)}-U_{j-2n}
+\]
+and the starting polynomials are
+\[
+U_{0}=1,\quad U_{1}=y,\quad U_{2}=x\quad.
+\]
+
+
+\subsection{Recurrence relations for Cartesian derivatives for Zernike U-polynomials}
+
+The general recursive relations for partial derivatives are:
+\[
+\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial x}=\tilde{n}U_{\tilde{n}-1,\tilde{m}}+\tilde{n}U_{\tilde{n}-1,\tilde{m}-1}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial x}
+\]
+and
+\[
+\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial y}=\tilde{n}U_{\tilde{n}-1}-\tilde{n}U_{\tilde{n}-1,\tilde{n}-\tilde{m}}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial y}
+\]
+But, as before, there are several exceptions:
+\begin{itemize}
+\item for $\tilde{m}=0$
+\[
+\frac{\partial U_{\tilde{n},0}}{\partial x}=\tilde{n}U_{\tilde{n}-1,0}\qquad\frac{\partial U_{\tilde{n},0}}{\partial y}=\tilde{n}U_{\tilde{n}-1,\tilde{n}-1}
+\]
+\item for $\tilde{m}=\tilde{n}$
+\[
+\frac{\partial U_{\tilde{n},\tilde{n}}}{\partial x}=-\tilde{n}U_{\tilde{n}-1,\tilde{n}-1}\qquad\frac{\partial U_{\tilde{n},\tilde{n}}}{\partial y}=-\tilde{n}U_{\tilde{n}-1,0}
+\]
+\item for $\tilde{n}$ odd and $\tilde{m}=\dfrac{\tilde{n}-1}{2}$
+\[
+\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial x}=\tilde{n}U_{\tilde{n}-1,\tilde{m}-1}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial x}\qquad\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial y}=\tilde{n}U_{\tilde{n}-1,\tilde{n}-\tilde{m}-1}-\tilde{n}U_{\tilde{n}-1,\tilde{n}-\tilde{m}}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial y}
+\]
+\item for $\tilde{n}$ odd and $\tilde{m}=\dfrac{\tilde{n}-1}{2}+1$
+\[
+\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial x}=\tilde{n}U_{\tilde{n}-1,\tilde{m}}+\tilde{n}U_{\tilde{n}-1,\tilde{m}-1}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial x}\qquad\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial y}=\tilde{n}U_{\tilde{n}-1,\tilde{n}-\tilde{m}-1}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial y}
+\]
+\item for $\tilde{n}$ even and $\tilde{m}=\dfrac{\tilde{n}}{2}$
+\[
+\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial x}=2\tilde{n}U_{\tilde{n}-1,\tilde{m}}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial x}\qquad\frac{\partial U_{\tilde{n},\tilde{m}}}{\partial y}=2\tilde{n}U_{\tilde{n}-1,\tilde{n}-\tilde{m}-1}+\frac{\partial U_{\tilde{n}-2,\tilde{m}-1}}{\partial y}
+\]
+\end{itemize}
+The starting expressions for the Cartesian derivatives being:
+\[
+\frac{\partial U_{0,0}}{\partial x}=\frac{\partial U_{0,0}}{\partial y}=0,\quad\frac{\partial U_{1,0}}{\partial x}=\frac{\partial U_{1,1}}{\partial y}=0,\quad\frac{\partial U_{1,1}}{\partial x}=\frac{\partial U_{1,0}}{\partial y}=1
+\]
+
+
+\subsubsection{Recurrence relations for the OSA/ANSI scheme}
+
+If we now ``translate'' the previous equations into the OSA scheme,
+we obtain that the general recurrence relations for partial derivativess
+are
+\[
+\frac{\partial U_{n,m}}{\partial x}=nU_{n-1,m+1}+nU_{n-1,m-1}+\frac{\partial U_{n-2,m}}{\partial x}
+\]
+and
+\[
+\frac{\partial U_{n,m}}{\partial y}=nU_{n-1,-\left(m+1\right)}-nU_{n-1,1-m}+\frac{\partial U_{n-2,m}}{\partial y}
+\]
+the exceptions are give by
+\begin{itemize}
+\item for $m=-n$$\qquad\dfrac{\partial U_{n,-n}}{\partial x}=nU_{n-1,1-n}\qquad\dfrac{\partial U_{n,-n}}{\partial y}=nU_{n-1,n-1}$
+\item for $m=n$$\qquad\dfrac{\partial U_{n,-n}}{\partial x}=nU_{n-1,n-1}\qquad\dfrac{\partial U_{n,-n}}{\partial y}=-nU_{n-1,1-n}$
+\item for $m=-1$
+\[
+\dfrac{\partial U_{n,-1}}{\partial x}=nU_{n-1,-2}+\dfrac{\partial U_{n-2,-1}}{\partial x}\qquad\dfrac{\partial U_{n,-1}}{\partial y}=nU_{n-1,0}-nU_{n-1,2}+\dfrac{\partial U_{n-2,-1}}{\partial y}
+\]
+\item for $m=1$
+\[
+\dfrac{\partial U_{n,1}}{\partial x}=nU_{n-1,2}+nU_{n-1,0}+\dfrac{\partial U_{n-2,1}}{\partial x}\qquad\dfrac{\partial U_{n,1}}{\partial y}=nU_{n-1,-2}+\dfrac{\partial U_{n-2,1}}{\partial y}
+\]
+\item for $m=0$
+\[
+\dfrac{\partial U_{n,0}}{\partial x}=2nU_{n-1,1}+\dfrac{\partial U_{n-2,0}}{\partial x}\qquad\dfrac{\partial U_{n,0}}{\partial y}=2nU_{n-1,-1}+\dfrac{\partial U_{n-2,0}}{\partial y}
+\]
+\end{itemize}
+and the starting derivatives are
+\[
+\frac{\partial U_{0,0}}{\partial x}=\frac{\partial U_{0,0}}{\partial y}=0,\quad\frac{\partial U_{1,-1}}{\partial x}=\frac{\partial U_{1,1}}{\partial y}=0,\quad\frac{\partial U_{1,1}}{\partial x}=\frac{\partial U_{1,-1}}{\partial y}=1
+\]
+
+
+\subsubsection{Recurrence relations with a single index}
+
+Let us start now with the exceptions:
+\begin{itemize}
+\item for $m=-n\quad\Rightarrow\quad j=\dfrac{\left(n+1\right)n}{n}$
+\[
+\frac{\partial U_{j}}{\partial x}=nU_{j-n}\qquad\frac{\partial U_{j}}{\partial y}=nU_{j-1}
+\]
+\item for $m=n\quad\Rightarrow\quad j=\dfrac{n\left(n+3\right)}{2}$
+\[
+\frac{\partial U_{j}}{\partial x}=nU_{j-\left(n+1\right)}\qquad\frac{\partial U_{j}}{\partial y}=-nU_{j-2n}
+\]
+\item for $m=-1\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)-1}{2}$
+\[
+\frac{\partial U_{j}}{\partial x}=nU_{j-\left(n+1\right)}+\frac{\partial U_{j-2n}}{\partial x}\qquad\frac{\partial U_{j}}{\partial y}=nU_{j-n}-nU_{j-\left(n-1\right)}+\frac{\partial U_{j-2n}}{\partial y}
+\]
+\item for $m=1\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)+1}{2}$
+\[
+\frac{\partial U_{j}}{\partial x}=nU_{j-n}+nU_{j-\left(n+1\right)}+\frac{\partial U_{j-2n}}{\partial x}\qquad\frac{\partial U_{j}}{\partial y}=nU_{j-\left(n+2\right)}+\frac{\partial U_{j-2n}}{\partial y}
+\]
+\item for $m=0\quad\Rightarrow\quad j=\dfrac{n\left(n+2\right)}{2}$
+\[
+\frac{\partial U_{j}}{\partial x}=2nU_{j-n}+\frac{\partial U_{j-2n}}{\partial x}\qquad\frac{\partial U_{j}}{\partial y}=2nU_{j-\left(n+1\right)}+\frac{\partial U_{j-2n}}{\partial y}
+\]
+\end{itemize}
+The general case is
+\[
+\frac{\partial U_{j}}{\partial x}=nU_{j-n}+nU_{j-\left(n+1\right)}+\frac{\partial U_{j-2n}}{\partial x}\qquad\frac{\partial U_{j}}{\partial y}=nU_{j-\left(n+m+1\right)}-nU_{j-\left(n+m\right)}+\frac{\partial U_{j-2n}}{\partial y}
+\]
+and the starting derivatives are
+
+\[
+\frac{\partial U_{0}}{\partial x}=\frac{\partial U_{0}}{\partial y}=0,\quad\frac{\partial U_{1}}{\partial x}=\frac{\partial U_{2}}{\partial y}=0,\quad\frac{\partial U_{2}}{\partial x}=\frac{\partial U_{1}}{\partial y}=1
+\]
+
+\begin{thebibliography}{1}
+\bibitem{Andersen_2018}Andersen T.B., \href{https://doi.org/10.1364/OE.26.018878}{Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates}.
+\emph{Optic Express 26}(15), 18878-18896 (2018).
+
+\bibitem{Rimmer=000026Wyant_1975}Rimmer M.P. \& Wyant J.C., \href{https://wp.optics.arizona.edu/jcwyant/wp-content/uploads/sites/13/2016/08/LSI-Analysis.pdf}{Evaluation of large aberrations using a lateral-shear interferometer having variable shear}.
+\emph{Appl. Opt. 14}(1), 142--150 (1975).
+
+\bibitem{Thibos_et_al_2002}Thibos, L.N, Applegate, R.A., Schwiegerling,
+J.T. \& Webb, R., \href{https://pdfs.semanticscholar.org/6e71/c63afa9960dc6e452f3ae2c8b20e2185ba88.pdf}{Standards for reporting the optical aberrations of eyes}.
+\emph{Journal of refractive surgery}, \emph{18}(5), S652-S660 (2002).
+
+\bibitem{ANSI-2017}ANSI Z80.28--2017. American National Standards
+of Ophthalmics:methods for reporting optical aberrations of eyes.
+\end{thebibliography}
+
+\end{document}
diff -r f97c278a5ac0 -r 437f7f1fae96 inst/zernike_osa_ansi_to_mn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/zernike_osa_ansi_to_mn.m	Sun Aug 04 20:10:40 2019 +0200
@@ -0,0 +1,135 @@
+## Copyright (C) 2019 José Ramom Flores das Seixas 
+##                                      <jose.ramom.flores.das.seixas@gmail.com>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {[@var{m}, @var{n}] =} zernike_osa_ansi_to_mn (@var{j})
+## Convert OSA/ANSI single-index @var{j} to double index @var{m} (Azimuthal 
+## degree) and @var{n} (Radial degree).
+##
+## Example
+## @example
+## [m,n] = zernike_osa_ansi_to_mn(4)
+##     @result{} [0, 2]
+## @end example
+##
+## References:
+##
+## @enumerate
+## @item Thibos, L.N, Applegate, R.A., Schwiegerling, J.T. & Webb, R., 
+##       Standards for reporting the optical aberrations of eyes. Journal of
+##       refractive surgery , 18 (5), S652-S660 (2002).
+## @item @url{https://en.wikipedia.org/wiki/Zernike_polynomials#OSA/ANSI_standard_indices,
+## OSA/ANSI standard indices of Zernike polynomials}, last retrieved on
+##       July 2109.
+## @end enumerate
+##
+## @seealso{zernike_noll_to_mn, zernikes_and_derivatives_cartesian_OSA, zernike_cartesian, zernike_name, zernike_polar, zernike_R_poly}
+##
+## @end deftypefn
+
+function [m, n] = zernike_osa_ansi_to_mn (j)
+  if (nargin != 1 || nargout != 2)
+    print_usage ();
+  elseif (any (j(:) < 0 | j(:) != fix (j(:))) || ! isvector (j))
+    error ("zernike_osa_ansi_to_mn: j has to be a vector with integers >=0");
+  endif
+  n  = ceil((-3+sqrt(9+8*j))/2);
+  m  = 2*j - n .* (n+2);
+endfunction
+
+%!test
+%! [m,n]=zernike_osa_ansi_to_mn(0);
+%! assert([m n],[0 0])
+%! [m,n]=zernike_osa_ansi_to_mn(1);
+%! assert([m n],[-1 1])
+%! [m,n]=zernike_osa_ansi_to_mn(2);
+%! assert([m n],[1 1])
+%! [m,n]=zernike_osa_ansi_to_mn(3);
+%! assert([m n],[-2 2])
+%! [m,n]=zernike_osa_ansi_to_mn(4);
+%! assert([m n],[0 2])
+%! [m,n]=zernike_osa_ansi_to_mn(5);
+%! assert([m n],[2 2])
+%! [m,n]=zernike_osa_ansi_to_mn(6);
+%! assert([m n],[-3 3])
+%! [m,n]=zernike_osa_ansi_to_mn(7);
+%! assert([m n],[-1 3])
+%! [m,n]=zernike_osa_ansi_to_mn(8);
+%! assert([m n],[1 3])
+%! [m,n]=zernike_osa_ansi_to_mn(9);
+%! assert([m n],[3 3])
+%! [m,n]=zernike_osa_ansi_to_mn(10);
+%! assert([m n],[-4 4])
+%! [m,n]=zernike_osa_ansi_to_mn(11);
+%! assert([m n],[-2 4])
+%! [m,n]=zernike_osa_ansi_to_mn(12);
+%! assert([m n],[0 4])
+%! [m,n]=zernike_osa_ansi_to_mn(13);
+%! assert([m n],[2 4])
+%! [m,n]=zernike_osa_ansi_to_mn(14);
+%! assert([m n],[4 4])
+%! [m,n]=zernike_osa_ansi_to_mn(15);
+%! assert([m n],[-5 5])
+%! [m,n]=zernike_osa_ansi_to_mn(16);
+%! assert([m n],[-3 5])
+%! [m,n]=zernike_osa_ansi_to_mn(17);
+%! assert([m n],[-1 5])
+%! [m,n]=zernike_osa_ansi_to_mn(18);
+%! assert([m n],[1 5])
+%! [m,n]=zernike_osa_ansi_to_mn(19);
+%! assert([m n],[3 5])
+%! [m,n]=zernike_osa_ansi_to_mn(20);
+%! assert([m n],[5 5])
+%! [m,n]=zernike_osa_ansi_to_mn(21);
+%! assert([m n],[-6 6])
+%! [m,n]=zernike_osa_ansi_to_mn(22);
+%! assert([m n],[-4 6])
+%! [m,n]=zernike_osa_ansi_to_mn(23);
+%! assert([m n],[-2 6])
+%! [m,n]=zernike_osa_ansi_to_mn(24);
+%! assert([m n],[0 6])
+%! [m,n]=zernike_osa_ansi_to_mn(25);
+%! assert([m n],[2 6])
+%! [m,n]=zernike_osa_ansi_to_mn(26);
+%! assert([m n],[4 6])
+%! [m,n]=zernike_osa_ansi_to_mn(27);
+%! assert([m n],[6 6])
+%! [m,n]=zernike_osa_ansi_to_mn(28);
+%! assert([m n],[-7 7])
+%! [m,n]=zernike_osa_ansi_to_mn(29);
+%! assert([m n],[-5 7])
+%! [m,n]=zernike_osa_ansi_to_mn(30);
+%! assert([m n],[-3 7])
+%! [m,n]=zernike_osa_ansi_to_mn(31);
+%! assert([m n],[-1 7])
+%! [m,n]=zernike_osa_ansi_to_mn(32);
+%! assert([m n],[1 7])
+%! [m,n]=zernike_osa_ansi_to_mn(33);
+%! assert([m n],[3 7])
+%! [m,n]=zernike_osa_ansi_to_mn(34);
+%! assert([m n],[5 7])
+%! [m,n]=zernike_osa_ansi_to_mn(35);
+%! assert([m n],[7 7])
+
+## vector test
+%!test
+%! [m,n]=zernike_osa_ansi_to_mn([2 7 15 35]);
+%! assert(m,[1 -1 -5 7])
+%! assert(n,[1 3 5 7])
+
+
+
+
diff -r f97c278a5ac0 -r 437f7f1fae96 inst/zernikes_and_derivatives_cartesian_OSA.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/zernikes_and_derivatives_cartesian_OSA.m	Sun Aug 04 20:10:40 2019 +0200
@@ -0,0 +1,215 @@
+## Copyright (C) 2019 José Ramom Flores das Seixas 
+##                                      <jose.ramom.flores.das.seixas@gmail.com>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {[@var{Z}, @var{dZx}, @var{dZy}] =} zernikes_and_derivarives_cartesian_OSA (@var{x}, @var{y}, @var{n})
+## @deftypefnx {} {[@var{Z}, @var{dZx}, @var{dZy}] =} zernikes_and_derivarives_cartesian_OSA (@var{x}, @var{y}, @var{n}, @var{nan_zero})
+## Return the cartesian Zernike's pollynomials and its partial 
+## derivatives up to radial degree @var{n}, i.e. until Z[@var{n},@var{n}]
+## 
+## @var{x} is a matrix with the X coordinates of the points where the Zernike's
+## polynomials and its derivatives are computed.
+## @var{y} is a matrix with the Y coordinates of the same points.
+## @var{n} is an integer with the maximum radial degree desired.
+## @var{nan_zero} is a string that determines the values of polynomial 
+## and derived values outside the radio unit circle.
+##
+## Strictly, the polynoms are only defined for 0 <= X²+Y² <= 1.
+## If variable @var{nan_zero} = 'nan', the values of the polynomials for which
+## it is verified that (X²+Y²)>1 are set = NaN.
+## If variable @var{nan_zero} = 'zero', the values of the polynomials for which
+## it is verified that (X²+Y²)>1 are set = 0.
+##
+## @var{Z} is a 3D matrix. Each page contains a 2D matrix with the values of a 
+## Zernike's polynomial at the points defined by @var{x} and @var{y}.  
+##
+## @var{dZx} is a 3D matrix. Each page contains the values of the 
+## partial derivative in x.
+##
+## @var{dZy} is a 3D matrix. Each page contains the values of the 
+## partial derivative in y.
+## 
+## It should be noted that in standard OSA/ANSI the simple-index j starts 
+## at 0, but in octave the indices of the vectors and matrices start at 1.
+## So that page 1 of the 3D Z, dZx and dZy matrices corresponds to the 
+## single-index j = 0, and therefore to the double-index m = 0 and n = 0.
+## Page 2 corresponds to j = 1, page 3 --> j = 2, etc.
+## 
+## Example
+## @example
+## x = linspace(-1,1,101);
+## [X,Y] = meshgrid(x,x);
+## [Z,dZx,dZy] = zernikes_and_derivatives_cartesian_OSA (X,Y,7,'zero');
+## Z_00 = Z(:,:,1);
+##    # Z_00 is a 2D matrix with the values of Zernike's polynomial
+##    # with simple-index j = 0, and double-index m = 0 & n = 0. 
+## dZx_-24 = dZx(:,:,11);
+##    # Z_-44 is a 2D matrix with the values of the partial 
+##    # derivative in x of Zernike's polynomial with 
+##    # simple-index j = 10, and double-index m = -4 & n = 4.
+## @end example
+##
+## Run the demo to see a more complete example.
+##
+
+##
+## Size of @var{x} must be equal size of @var{y}.
+##
+## References:
+##
+## @enumerate
+## @item Andersen T.B., @url{https://doi.org/10.1364/OE.26.018878, "Efficient 
+##       and robust recurrence relations for the Zernike circle polynomials and
+##       their derivatives in Cartesian coordinates"}. Optic Express 26(15), 
+##       18878-18896 (2018).
+## @item Thibos, L.N, Applegate, R.A., Schwiegerling, J.T. & Webb, R., 
+##       Standards for reporting the optical aberrations of eyes. Journal of 
+##       refractive surgery, 18(5), S652-S660 (2002).
+## @end enumerate
+##
+## @seealso{zernike_osa_ansi_to_nm, zernike_cartesian, zernike_name,
+##          zernike_polar, zernike_R_poly}
+## @end deftypefn
+
+function [Z,dZx,dZy] = zernikes_and_derivatives_cartesian_OSA (x, y, N, nan_zero = 'NaN');
+  if (nargin < 3 || nargin > 4)
+    print_usage ();
+  elseif (! isscalar (N) || N < 1 || N != fix (N))
+    error ("zernike_and_derivatives_cartesian_OA: Maximum radial degree must be a integer >=1");
+  elseif (any (size (x) != size (y)))
+    error ("zernike_and_derivatives_cartesian_OA: X and Y must have the same size");
+  endif
+  r2 = x.^2 + y.^2;  out_p = find(r2>1);
+  MaxJ = N*(N+3)/2;
+  U1 = ones(size(x));
+  switch(nan_zero)
+    case "zero"
+      x(out_p) = 0;
+      y(out_p) = 0;
+      U1(out_p) = 0;
+    case "NaN"
+      x(out_p) = NaN;
+      y(out_p) = NaN;
+      U1(out_p) = NaN;
+  endswitch
+  U = zeros ([size(x) (MaxJ+1)]);    dUx = U;   dUy = U;
+  Z = U;		     dZx = U;		     dZy = U;
+  U(:,:,1) = U1;	dUx(:,:,1) = U1*0;	dUy(:,:,1) = U1*0;	
+  Z(:,:,1) = U1;	dZx(:,:,1) = U1*0;	dZy(:,:,1) = U1*0;
+  U(:,:,2) = y;	dUx(:,:,2) = U1*0;	dUy(:,:,2) = U1;
+  Z(:,:,2) = 2*y;	dZx(:,:,2) = U1*0;	dZy(:,:,2) = 2*U1;
+  U(:,:,3) = x;	dUx(:,:,3) = U1;	dUy(:,:,3) = U1*0;
+  Z(:,:,3) = 2*x;	dZx(:,:,3) = 2*U1;	dZy(:,:,3) = U1*0;
+  for k = 4:(MaxJ+1)
+    [m, n] = zernike_osa_ansi_to_mn (k-1);
+    switch(m)
+      case -n   
+	   U(:,:,k) = x.*U(:,:,(k-n))+y.*U(:,:,(k-1));
+	   dUx(:,:,k) = n*U(:,:,(k-n));
+	   dUy(:,:,k) = n*U(:,:,(k-1));
+      case  n   
+	   U(:,:,k) = x.*U(:,:,(k-(n+1)))-y.*U(:,:,(k-2*n));
+	   dUx(:,:,k) = n*U(:,:,(k-(n+1)));
+	   dUy(:,:,k) = -n*U(:,:,(k-2*n));
+      case -1   
+	   U(:,:,k) = x.*U(:,:,(k-(n+1)))+y.*U(:,:,(k-n))-y.*U(:,:,(k-(n-1)))-...
+		         U(:,:,(k-2*n));
+	   dUx(:,:,k) = n*U(:,:,(k-(n+1)))+dUx(:,:,(k-2*n));
+	   dUy(:,:,k) = n*U(:,:,(k-n))-n*U(:,:,(k-(n-1)))+dUy(:,:,(k-2*n));
+      case  1   
+	   U(:,:,k) = x.*U(:,:,(k-n))+x.*U(:,:,(k-(n+1)))+y.*U(:,:,(k-(n+2)))-...
+		         U(:,:,(k-2*n));
+	   dUx(:,:,k) = n*U(:,:,(k-n))+n*U(:,:,(k-(n+1)))+dUx(:,:,(k-2*n));
+	   dUy(:,:,k) = n*U(:,:,(k-(n+2)))+dUy(:,:,(k-2*n));
+      case  0   
+	   U(:,:,k) = 2*x.*U(:,:,(k-n))+2*y.*U(:,:,(k-(n+1)))-U(:,:,(k-2*n));
+	   dUx(:,:,k) = 2*n*U(:,:,(k-n))+dUx(:,:,(k-2*n));
+	   dUy(:,:,k) = 2*n*U(:,:,(k-(n+1)))+dUy(:,:,(k-2*n));
+      otherwise 
+	   U(:,:,k) = x.*U(:,:,(k-n))+y.*U(:,:,(k-(n+m+1)))+x.*U(:,:,(k-(n+1)))-...
+		         y.*U(:,:,(k-(n+m)))-U(:,:,(k-2*n));
+	   dUx(:,:,k) = n*U(:,:,(k-n))+n*U(:,:,(k-(n+1)))+dUx(:,:,(k-2*n));
+	   dUy(:,:,k) = n*U(:,:,(k-(n+m+1)))-n*U(:,:,(k-(n+m)))+dUy(:,:,(k-2*n));
+    endswitch
+    Nnm = sqrt(2*(n+1));
+    if m == 0 
+       Nnm = Nnm/sqrt(2);
+    endif
+    Z(:,:,k) = Nnm * U(:,:,k);
+    dZx(:,:,k) = Nnm * dUx(:,:,k);
+    dZy(:,:,k) = Nnm * dUy(:,:,k);
+  endfor
+endfunction
+
+%!demo
+%! Nmgr = 7;                        # Maximum radial degree number, default = 7
+%! Nmodos = ((Nmgr+3)*Nmgr/2)+1;    # Mode number, default = 36
+%!
+%! pr = 20;                         # Pupil radius, mm
+%! # Let us suppose that we obtain the heights of a biconical surface
+%! # for a series of points in the xy plane
+%! Rx = 0.002;			          # mm⁻¹
+%! px = 0.8;
+%! Ry = 0.001;			          # mm⁻¹
+%! py = 0.5;
+%! x = pr*(2*rand(200,1) - 1);	     # mm
+%! y = pr*(2*rand(200,1) - 1);	     # mm
+%! surface = (Rx*x.^2 + Ry*y.^2)/(1+sqrt(1-(px*(Rx*x).^2+py*(Ry*y).^2)));
+%! sag = max(max(surface));
+%! surface = sag - surface;         # surface height, mm
+%!
+%! # Coordinates normalization
+%! xnorm=x/pr;
+%! ynorm=y/pr;
+%! usefuldata = find(sqrt(x.^2+y.^2)<pr);     # data inside the pupil circle
+%! ixnorm = xnorm(usefuldata);
+%! iynorm = ynorm(usefuldata);
+%! isurface = surface(usefuldata);
+%!
+%! # Fitting the surface			Note that ixnomr and iynorm are vectors
+%! [Z,dZx,dZy] = zernikes_and_derivatives_cartesian_OSA (ixnorm, iynorm, Nmgr);
+%! Mrec = inv(Z'*Z)*Z';             # LSQ fitting matrix
+%! coefZ = Mrec*isurface;           # Fitting Zernike's coefficients
+%!
+%! # Reconstructed surface and partial derivatives
+%! nx = linspace(-1,1,129);
+%! [Xn,Yn] = meshgrid(nx,nx);
+%! [Z,dZx,dZy] = zernikes_and_derivatives_cartesian_OSA (Xn, Yn, Nmgr);
+%! S = zeros(size(Xn));		# Note that Xn and Yn are square matrices
+%! dSx = zeros(size(Xn));
+%! dSy = zeros(size(Xn));
+%! for m = 1:Nmodos
+%!    S = S + Z(:,:,m)*coefZ(m);		# Reconstructed surface
+%!    dSx = dSx + dZx(:,:,m) * coefZ(m);	# Derivative of S with respect to x
+%!    dSy = dSy + dZy(:,:,m) * coefZ(m);	# Derivative of S with respect to y
+%! endfor
+%!
+%! # Plotting the reconstructed surface and their derivatives
+%! nx = linspace(-pr,pr,129);
+%! [Xn,Yn] = meshgrid(nx,nx);
+%! colormap ("default");
+%! subplot (3,2,1:4)
+%! mesh(Xn,Yn,S);
+%! title("Reconstructed surface");
+%! xlabel("x (mm)");      ylabel("y (mm)");        zlabel("z (mm)");
+%! subplot(3,2,5)
+%! imagesc(nx,nx,dSx);
+%! axis square;      axis xy;       colorbar;
+%! title("Derivative of Surface with respect to x");
+%! subplot(3,2,6)
+%! imagesc(nx,nx,dSy);
+%! axis square;      axis xy;       colorbar;
+%! title("Derivative of Surface with respect to y");
