#ifndef ZERNIKE_C
#define ZERNIKE_C
/**
 @author Richard J. Mathar
 @since 2007-03-23
 <a href="http://www.strw.leidenuniv.nl/~mathar">Richard J. Mathar homepage</a>.
*/
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <string.h>

#include "Zernike.h"

using namespace std ;

/** Ctor.
* @param[in] nodaln the nodal quantum number
* @since 2007-03-23
*/
Zernike::Zernike(const int nodaln) : n(nodaln)
{
	/* coefficients are only non-zero for even n-m. So we loop only
	* over these.
	*/
	for(int m= n %2 ; m <=n ; m += 2)
	{
		const int nplusmhalf = (n+m)/2 ;
		const int nminusmhalf = (n-m)/2 ;

		/* start with the coefficient at s=(n-m)/2, corresponding to r^m and x^m.
		* This is the coefficient in front of \f$r^{n-2s}\f$ in eq (2) of Noll's paper.
		* (n-s)! = [(n+m)/2]!. s!=[(n-m)/2]!. [(n+m)/2-s]!=m!. [(n-m)/2-s]!=1
		*/
		double c = ( nminusmhalf % 2) ? -1. : 1. ;
		c *= binom(nplusmhalf,m) ;

		/* so far, c is the polynomial coefficient for \f$r\f$. Move over to the
		* coefficient for \f$x=r/2\f$.
		*/
		c *= pow(2.,m) ;

		/* We are interested not in Noll's \f$R_n^m\f$ but in the normalized \f$2\surd()R_n^m\f$
		*/
		c *= 2.*sqrt((double)(2*n+2)) ;

		vector<double> Rm ;

		/* we arrange for the lowest power x^m coming first, the power x^n coming last.
		*/
		for(int s= nminusmhalf ; s >= 0 ; s--)
		{
			Rm.push_back(c) ;
			/* Update: (-)^s  brings a minus sign. Transition (n-s)!->(n-s+1)! implied by s->s-1 brings multiplication
			* by (n-s+1). Transition 1/s!->1/(s-1)! brings multiplication by s. Transition
			* 1/[(n+m)/2-s]! -> 1/[..-s+1]! brings division through (n+m)/2-s+1. Transition
			* 1/[(n-m)/2-s]! -> 1/[..-s+1]! brings division through (n-m)/2-s+1. The modification 
			* to the x powers means (2x)^(n-2s)-> (2x)^(n-2s+2) brings a factor of 4.
			*/
			c *= -4*(n-s+1)*s/(nplusmhalf-s+1.)/(nminusmhalf-s+1.) ;
		}

		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " Z(n,m) " << nodaln << " " << m << endl ;
		*/

		coef.push_back(Rm); 
	}
}

/** The value at some distance from the origin.
* @param[in] x the distance from the origin in the range 0 to 1/2.
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @return the value at that distance in the unit normalization.
* @since 2007-03-23
*/
double Zernike::at(const double x, const int m) const
{

	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " Z(n,m) query " << n << " " << m << endl ;
	*/

	/* request outside the valid distance from the origin or with invalid azimuthal number
	* are silently answered with a value of zero.
	*/
	if ( x>=0. && x <=0.5 && m>= 0 && m <=n && (n-m) %2 ==0)
	{
		double resul =0. ;

		/* sum the coefficients times x^t, t=m..n in steps of 2, coefficients in coef[m/2][(t-m)/2].
		* Use some Horner scheme to do the multiplications.
		* resul = coef[][0]*x^m+coef[][1]*x^(m+2)+...+coef[][(n-m)/2]*x^n.
		*       = x^m { [][0]+coef[][1]*x^2+...+coef[][(n-m)/2]*x^(n-m) }
		*/
		const int mhalf = m/2 ;
		const double x2 = x*x ;
		for(int nminmhalf=(n-m)/2; nminmhalf>=0 ; nminmhalf--)
		{
			resul *= x2 ;
			resul += coef[mhalf][nminmhalf] ;
		}

		/* the residual x^m is accumulated also
		*/
		return resul* pow(x,m) ;
	}
	else
		return 0. ;
}

#if 0

/** The value at some distance from the origin.
* @param[in] x the distance from the origin in the range 0 to 1/2.
* @param[in] n the nodal quantum number
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @return the value at that distance in the unit normalization.
* The functionality is the same as with a Zernike::Zernike(n) followed by Zernike::at(x,m) .
* @since 2007-03-23
*/
double Zernike::at(const double x, const int n, const int m)
{
	/* request outside the valid distance from the origin or with invalid azimuthal number
	* are silently answered with a value of zero.
	*/
	if ( x>=0. && x <=0.5 && m>= 0 && m <=n && (n-m) %2 ==0)
	{
		const double r = 2.*x ;
		const double r2 = r*r ;
		const int nplusmhalf = (n+m)/2 ;
		const int nminusmhalf = (n-m)/2 ;

		/* start with the coefficient at s=0, corresponding to r^n.
		* This is the coefficient in front of \f$r^{n-2s}\f$ in eq (2) of Noll's paper.
		* (n-s)! = n!. s!=1. [(n+m)/2-s]!=[(n+m)/2]!. [(n-m)/2-s]!=[(n-m)/2]!
		*/
		double c = binom(n,nminusmhalf) ;

		double resul =0. ;
		/* sum the coefficients times r^s, s=0..(n-m)/2 in steps of 2.
		* Use some Horner scheme to do the multiplications.
		*/
		for(int s=0 ; s <= nminusmhalf ; s++)
		{
			resul *= r2 ;
			resul += c ;

			/* Update: (-)^s  brings a minus sign. Transition (n-s)!->(n-s-1)! implied by s->s+1 brings division
			* through (n-s). Transition 1/s!->1/(s+1)! brings division through s. Transition
			* 1/[(n+m)/2-s]! -> 1/[..-s-1]! brings multiplication by (n+m)/2-s. Transition
			* 1/[(n-m)/2-s]! -> 1/[..-s-1]! brings multiplication by (n-m)/2-s.
			*/
			c *= -(nplusmhalf-s)*(nminusmhalf-s)/((double)(n-s)*s) ;
		}

		/* the residual r^m is accumulated also.
		* We are interested not in Noll's \f$R_n^m\f$ but in the normalized \f$2\surd()R_n^m\f$
		*/
		return resul*2.*M_SQRT2*sqrt((double)(n+1))*pow(r,m) ;
	}
	else
		return 0. ;
}

#endif /* 0 */

/** The coefficient of a specified power of x.
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @param[in] i the power of \f$x\f$
* @return the coefficient in front of \f$x^i\f$.
* @since 2007-03-23
*/
double Zernike::powcoef(const int m, const int i) const
{

	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " Z(n,m) query " << n << " " << m << " x^" << i << endl ;
	*/

	/* request outside the valid distance from the origin or with invalid azimuthal number
	* are silently answered with a value of zero.
	*/
	if ( m>= 0 && m <=n && (n-m) %2 ==0 && i>=m && i<=n && (i-m) %2 == 0)
	{
		/* the coefficients in front of x^i, i=m..n in steps of 2, are in coef[m/2][(i-m)/2].
		*/
		return coef[m/2][(i-m)/2] ;
	}
	else
		return 0. ;
}

/** The integral from xlo to xhi
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @param[in] xlo the lower limit
* @param[in] xhi the upper limit. Actually this will also return the correct
*   value of the upper limit is smaller than the lower one.
* @return the \f$\int_x Z dx\f$. Note that this does not involve any additional
*  power of \f$x\f$.
* @warn there is no check built in here that the two limits are in the interval
*   between 0 and 1/2.
* @since 2007-03-28
*/
double Zernike::intgrl(const int m, const double xlo, const double xhi) const
{
	/* request outside the valid distance from the origin or with invalid azimuthal number
	* are silently answered with a value of zero.
	*/
	if ( m>= 0 && m <=n && (n-m) %2 ==0 )
	{
		double resul = 0. ;
		/* the coefficients in front of x^i, i=m..n in steps of 2, are in coef[m/2][(i-m)/2].
		*/
		for(int i=m; i <= n ; i+= 2)
			resul += powcoef(m,i)*( pow(xhi,i+1)-pow(xlo,i+1))/(i+1) ;
		return resul ;
	}
	else
		return 0. ;
}

#if 0

/** The coefficient of a specified power of x.
* @param[in] n the nodal quantum number
* @param[in] m the azimuthal quantum number between 0 and n
* @param[in] i the power of \f$x\f$
* @return the coefficient in front of \f$x^i\f$.
* @since 2007-03-23
* The functionality is the same as with a Zernike::Zernike(n) followed by Zernike::powcoef(m,i) .
*/
double Zernike::powcoef(const int n, const int m, const int i)
{
	/* request outside the valid distance from the origin or with invalid azimuthal number
	* are silently answered with a value of zero.
	*/
	if ( m>= 0 && m <=n && (n-m) %2 ==0 && i>=m && i<=n && (i-m) %2 == 0)
	{
		const int s = (n-i)/2 ;

		const int nplusmhalf = (n+m)/2 ;
		const int nminusmhalf = (n-m)/2 ;

		/* Compute first coefficient for R_n^m(r=2x).
		* Noll's equation (2). Start with factor (-)^s
		*/
		double resul = ( s % 2) ? -1. : 1. ;

		/* rewrite (n-s)!/s! as binomial(n-s,s)*(n-2s)!. Then rewrite
		* (n-2s)!/ (nplusmhalf-s)! / (nminusmhalf-s)! as binomial(n-2s,nminusmhalf-s).
		* Note that this could also be nplusmhalf-s, but since there is a loop over the second
		* parameter in KarhLoeve::binom(), we try to minimize it. Finally use that n-2s=i.
		*/
		resul *= binom(n-s,s) * binom(i,nminusmhalf-s) ;

		/* The transition from r=2x to x induces another 2^i
		*/
		resul *= pow(2.,i) ;

		/*
		* We are interested not in Noll's \f$R_n^m\f$ but in the normalized \f$2\surd()R_n^m\f$
		*/
		return resul*2.*M_SQRT2*sqrt((double)(n+1)) ;
	}
	else
		return 0. ;
}

#endif /* 0 */

/** A overlap integral with another Zernike radial function.
* @param[in] m the azimuthal quantum number between 0 and n
* @param[in] oth the other radial function
* @param[in] othm the azimuthal quantum number of the other radial function
* @return the integral \f$\int_0^{1/2} Z(x)* othZ(x) dx\f$.
* @since 2007-03-23
* @warn this is not the ordinary overlap which would have a factor \f$x\f$ extra
*   from the Jacobian of the polar coordinates and which would be loosely related
*   to the usual orthogonalization properties.
*/
double Zernike::Overl(const int m, const Zernike & oth, const int othm) const
{
	double resul =0. ;
	for (int i= m; i<= n ; i += 2)
	for (int j= othm;j <= oth.n; j += 2)
		/* we use that the integral of x^i*x^j from 0 to 1/2 is 
		* x^(i+j+1)/(i+j+1).
		*/
		resul += powcoef(m,i)* oth.powcoef(othm,j)/ ((i+j+1)*pow(2.,i+j+1)) ;
	return resul ;
}

/** The overlap integral of the derivative with another radial function.
* @param[in] m the azimuthal quantum number between 0 and n
* @param[in] oth the other radial function in the overlap integral
* @param[in] othm the azimuthal quantum number of the other radial function
* @return the value of \f$\int_0^{1/2}dx x [\partial_x this (x)] other(x)\f$.
* @author Richard J. Mathar
* @since 2007-03-23
* @warning unlike Overl(), this function here is not commutative in the two functions.
*/
double Zernike::gradOverl(const int m, const Zernike & oth, const int othm) const
{
	double resul =0. ;
	for (int i= m; i<= n ; i += 2)
	for (int j= othm;j <= oth.n; j += 2)
		/* we use that the integral of x d_x(x^i)*x^j from 0 to 1/2 is 
		* i*x^(i+j+1)/(i+j+1). The only difference to Overl() is the
		* additional factor of i from the derivative of the first factor.
		*/
		resul += i*powcoef(m,i)* oth.powcoef(othm,j)/ ((i+j+1)*pow(2.,i+j+1)) ;
	return resul ;
}


/** Ctor.
* There are no parameters and we grow the internal set of basis functions as needed,
* which means as requested through one of the query functions.
* @since 2007-03-23
*/
ZernikeBase::ZernikeBase()
{
}

vector<Zernike> ZernikeBase::Z ;

/** Enlarge the basis set on request.
* @param[in] n the minimum nodal quantum number
* @since 2007-03-23
*/
void ZernikeBase::growto(const int n)
{
	/* if n is larger than the components created so far, enlarge the basis set
	*/
	while( n >= Z.size() )
	{
		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " Zernike::growto " << n << endl ;
		*/
		Zernike Zn(Z.size()) ;
		Z.push_back(Zn) ;
	}
}

/** The value at some distance from the origin.
* @param[in] x the distance from the origin in the range 0 to 1/2.
* @param[in] n the nodal quantum number
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @return the value at that distance in the unit normalization.
* The functionality is the same as with a Zernike::Zernike(n) followed by Zernike::at(x,m) .
* @since 2007-03-23
*/
double ZernikeBase::at(const double x, const int n, const int m)
{
	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " ZBas(n,m) query " << n << " " << m << endl ;
	*/

	/* if n is larger than the components created so far, enlarge the basis set
	*/
	growto(n) ;

	/* the static version which doesn't use the pre-computed values.
	* return Z[n].at(x,n,m) ;
	*/

	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " ZBas(n,m) query " << x << " " << n << " " << m << " " << Z[n].at(x,m) << endl ;
	*/

	return Z[n].at(x,m) ;
}

/** The value at some distance from the origin.
* @param[in] x the distance from the origin in the range 0 to 1/2.
* @param[in] c the expansion coefficients in the condensed format.
*   The condensed format is the one with all coefficiens of zero value omitted.
* @param[in] m the azimuthal quantum number between 0 and the nodal number.
* @return the value at that distance in the unit normalization.
* @since 2007-03-23
*/
double ZernikeBase::at(const double x, const vector<double> & coe, const int m)
{
	double resul=0. ;
	for(int c=0 ; c < coe.size() ; c++)
	{
		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " " << c << " " << coe[c] << " " << at(x,m+2*c,m) << endl ;
		*/

		resul += coe[c]* at(x,m+2*c,m) ;
	}
	return resul ;
}

/** The coefficient of a specified power of x.
* @param[in] n the nodal quantum number
* @param[in] m the azimuthal quantum number between 0 and n
* @param[in] i the power of \f$x\f$
* @return the coefficient in front of \f$x^i\f$.
* @since 2007-03-23
* The functionality is the same as with a Zernike::Zernike(n) followed by Zernike::powcoef(m,i) .
*/
double ZernikeBase::powcoef(const int n, const int m, const int i)
{
	/* if n is larger than the components created so far, enlarge the basis set
	*/
	growto(n) ;

	/* the static version which doesn't use the pre-computed values.
	* return Z[n].powcoef(n,m,i) ;
	*/

	return Z[n].powcoef(m,i) ;
}

/** The coefficient of a specified power of x.
* @param[in] n the nodal quantum number
* @param[in] m the azimuthal quantum number between 0 and n
* @param[in] xlo the lower limit of the integral
* @param[in] xhi the upper limit of the integral
* @return the \f$\int Z(n,m)\f$ between the two limits.
* @since 2007-03-28
* The functionality is the same as with a Zernike::Zernike(n) followed by Zernike::intgrl(m,xlo,xhi) .
*/
double ZernikeBase::intgrl(const int n, const int m, const double xlo, const double xhi)
{
	/* if n is larger than the components created so far, enlarge the basis set
	*/
	growto(n) ;

	return Z[n].intgrl(m,xlo,xhi) ;
}

/** A overlap integral between two Zernike radial functions.
* @param[in] n the nodal parameter of the first function
* @param[in] m the azimuthal quantum number of the first function
* @param[in] nprime the nodal parameter of the 2nd function
* @param[in] mprime the azimuthal quantum number of the 2nd function
* @return the integral \f$\int_0^{1/2} Z(x)* othZ(x) dx\f$.
* @since 2007-03-23
* @warn this is not the ordinary overlap which would have a factor \f$x\f$ extra
*   from the Jacobian of the polar coordinates and which would be loosely related
*   to the usual orthogonalization properties.
*/
double ZernikeBase::Overl(const int n, const int m, const int nprime, const int mprime)
{
	/* if n is larger than the components created so far, enlarge the basis set
	*/
	growto(n) ;
	growto(nprime) ;

	return Z[n].Overl(m, Z[nprime], mprime) ;
}

/** A overlap integral between two Zernike radial functions.
* @param[in] m the azimuthal quantum number of the first function
* @param[in] coe the expansion coefficients of the first function
* @param[in] mprime the azimuthal quantum number of the 2nd function
* @param[in] coeprime the expansion coefficients of the first function
* @return the integral \f$\int_0^{1/2} Z(x)* othZ(x) dx\f$.
* @since 2007-03-23
* @warn this is not the ordinary overlap which would have a factor \f$x\f$ extra
*   from the Jacobian of the polar coordinates and which would be loosely related
*   to the usual orthogonalization properties.
* @todo there is room for speedup since the the Overl() is a symmetric function of the
*   two essential parameter pairs. Combining the coe and coeprime instead of running
*   the loop over all independent pairs might save a factor of 2.
*   
*/
double ZernikeBase::Overl(const int m, const vector<double> & coe, const int mprime, const vector<double> &coeprime)
{
	double resul =0. ;
	for(int i = 0 ; i < coe.size(); i++)
		for(int j = 0 ; j < coeprime.size(); j++)
			resul += coe[i]*coeprime[j]*Overl(m+2*i,m,mprime+2*j,mprime) ;

	return resul ;
}

/** A overlap integral between two Zernike radial functions.
* @param[in] n the nodal parameter of the first function
* @param[in] m the azimuthal quantum number of the first function
* @param[in] nprime the nodal parameter of the 2nd function
* @param[in] mprime the azimuthal quantum number of the 2nd function
* @return the integral \f$\int_0^{1/2} x [grad_x Z(x)]* othZ(x) dx\f$.
* @since 2007-03-23
*/
double ZernikeBase::gradOverl(const int n, const int m, const int nprime, const int mprime)
{
	/* if n is larger than the components created so far, enlarge the basis set
	*/
	growto(n) ;
	growto(nprime) ;

	return Z[n].gradOverl(m, Z[nprime], mprime) ;
}

/** A overlap integral between two Zernike radial functions.
* @param[in] m the azimuthal quantum number of the first function
* @param[in] coe the expansion coefficients of the first function
* @param[in] mprime the azimuthal quantum number of the 2nd function
* @param[in] coeprime the expansion coefficients of the first function
* @return the integral \f$\int_0^{1/2} x [grad_x Z(x)]* othZ(x) dx\f$.
* @since 2007-03-23
*/
double ZernikeBase::gradOverl(const int m, const vector<double> & coe, const int mprime, const vector<double> &coeprime)
{
	double resul =0. ;
	for(int i = 0 ; i < coe.size(); i++)
		for(int j = 0 ; j < coeprime.size(); j++)
			resul += coe[i]*coeprime[j]*gradOverl(m+2*i,m,mprime+2*j,mprime) ;

	return resul ;
}

/** Binomial coefficients.
* @param[in] n the ``numerator'' of the coefficient.
* @param[in] m the ``denominator'' of the coefficient.
*    This is is limited to the range 0,1,2,...
* @return \f$x \choose l\f$
* @since 2007-03-23
*/
double Zernike::binom(const int n, int m)
{
        double resul =1. ;
	/* optimize with the formula C(n,m)=C(n,n-m) to minimize the loop count.
	*/
	if ( n < 2*m)
		m = n-m ;
        for(int c=1; c <=m ; c++)
                resul *= (n+1.-c)/c ;
        return resul ;
}


#endif /* ZERNIKE_C */