#ifndef KARHLOE_C
#define KARHLOE_C
/**
 @author Richard J. Mathar
 @since 2005-12-06
 <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 "KarhLoe.h"

extern "C" {
/** The clapack header file \c clapack.h.
* This is obtained from <a href="http://www.netlib.org/clapack">www.netlib.org/clapack</a>.
*/
#include <clapack.h>

#ifdef HAVE_GSL
#include <gsl_sf.h>
#endif
}


using namespace std ;

/** Print the Gamma coupling matrix for pupil resizes, if defined.
*/
#define KARHLOE_SHOW_RESIZE 2

/** Type 1 without lambda in the Gamma definition , Type 2 with
*/
#define KARHLOE_SHOW_TYPE 1

/** Macro to transpose matrix entries from the F77 to the C domain.
* The parameters are a row index (in the F77 program, but 0-based), a column index
* (in the F77 program, but 0-based), and the number of rows/columns of the square matrix.
* @since 2007-03-09
*/
#define KARHLOE_F77_FLAT(f77row,f77col,matdim) ((matdim)*(f77col)+(f77row))

/** \f$ 3.44=6.88/2 = (24/5*\Gamma(6/5))^{5/6}\f$.
* \latexonly
* See \cite[p 1372]{FriedJOSA56_1372}
* \endlatexonly
*/
#define FRIED_CONSTANT 3.441938591

/** The square root of Pi, \f$\sqrt\pi\f$.
*/
#define M_SQRT_PI 1.77245385090551602729816748

/** Two to the power of 11/6, \f$2^{11/6}\f$.
*/
#define M_2_POW_11_6 3.56359487256135721896090482236

/** Ctor.
* @param maxthetord
*/
KarhLoeve::KarhLoeve(const int maxthetord) : ord(maxthetord)
{
}

#ifdef KARHLOE_RADIAL_POLYN

/** Solve the KL radial equation to generate the eigenvalues and eigenfunctions.
* The symmetrized eigenvectors \f$\tilde K_p^q(x)\equiv \sqrt{x}K_p^q(x)\f$
* are used for the intermediate representation, represented by the variable \c Kpqsymm .
* The support mesh contains N basis functions \f$x^0\f$ up to \f$x^{N-1}\f$.
* The symmetrized kernel is \f$R_q(x,x')\equiv x'\tilde R_q(x,x')\f$.
* The KL equation is \f$ \int_{x'=0}^{1/2}\sqrt{xx'}\tilde R(x,x') \tilde K_p^q(x')dx'=\lambda_{pq}^2\tilde K_p^q(x)\f$.
* @param[in] q the azimuthal quantum number of the integral kernel
* @param[out] K the list of \c N eigenvectors (radial functions). 
*        \f$\int_{x=0}^{1/2}xK^2(x)dx=1=\sum_0^{N-1}x K^2(N)/(2N)\f$
* @param[out] lambda the list of \c N squared eigenvalues \f$\lambda_{pq}^2\f$
* @param[in] N the number of points in the interval \f$[0,1/2]\f$ on which the
*       integral kernel and eigenfunctions are discretized.
*/
void KarhLoeve::Kpq(const int q, double *K, double *lambda, const int N, ZernikeBase & zbase)
{
	/** matr contains the N by N matrix of the integral kernel.
	*/
	double *matr = new double[N*N] ;
	for(int i=0; i< N; i++)
	{
		for(int j=0; j<=i ; j++)
		{
			/** The values of \f$\tilde R_q(x,x')\f$ are collected by a
			* call to S(). The implicit representation is x^(1/2+q) times a polynomial of order N.
			*/
			matr[i+j*N]=matr[j+N*i] = -FRIED_CONSTANT*2*M_PI*binom(q-11./6.,q)*S(i,j,q,zbase) ;
		}
	}

	printf("Matrix q= %d N=%d\n",q,N) ;
	for(int i=0; i< N; i++)
	{
		for(int j=0; j<=i ; j++)
			printf("%.4e ",matr[i+j*N]) ;
		printf("\n") ;
	}

	/** Some intermediate variables are set up since we will be using the
	* standard LAPACK Fortran interface.
	*/
	char jobz[2] = "V" ;
	char uplo[2] = "U" ;
	integer lda = N ;
	integer lwork = 3*N ;
	integer info ;
	double *work = new double[lwork] ;

	/** The finite element version of the integral equation with symmetric matrix
	* is solved by a call to LAPACK's dsyev().
	*/
	dsyev_(jobz,uplo,&lda,matr,&lda,lambda,work,&lwork,&info) ;
	if ( info)
		cerr << __LINE__ << " dsyev info " << info << endl ;
	else
	{
		/* the result represents tilde K, but since the sqrt(x) was only
		* implicit, silently dropping this sqrt(x) does not initiate action.
		* The eigenvectors represent the prefactors of x^(q+i), i=0,1,2..N-1
		* and are normalized to a sum of squares equal unity, and we want
		* to have an integral over x*K^2=1 (x=0..1/2) = 
		* int x * sum_j,s k_j k_s x^{j+q}x^{s+q} = sum_j,s.
		*/
#if 0
		// re-normalize such that they are positive at x=1/2 : TODO
		for(int e=0; e < N ;e++)
		{
			if( matr[N-1+N*e] < 0.)	// last value of the eigenvector
				for(int j=0; j < N ; j++)
					matr[j+N*e] *= -1 ;

		}
#endif
		/* save the eigenvectors into return area
		* They are now normalized to sum_{0..N} ( x*K^2* Delta(x)=1 where
		* Delta(x)=0.5/N.
		*/
		memcpy(K,matr,N*N*sizeof(double)) ;
	}
	delete work ;
	delete matr ;
}

#else /* KARHLOE_RADIAL_POLYN */

/** Solve the KL radial equation to generate the eigenvalues and eigenfunctions.
* The symmetrized eigenvectors \f$\tilde K_p^q(x)\equiv \sqrt{x}K_p^q(x)\f$
* are used for the intermediate representation, represented by the variable \c Kpqsymm .
* The support mesh contains N points \f$x=(i+1/2)/(2N)\f$, \f$i=0\ldots N-1\f$, \f$0\le x\le 1/2\f$.
* The symmetrized kernel is \f$R_q(x,x')\equiv x'\tilde R_q(x,x')\f$.
* The KL equation is \f$ \int_{x'=0}^{1/2}\sqrt{xx'}\tilde R(x,x') \tilde K_p^q(x')dx'=\lambda_{pq}^2\tilde K_p^q(x)\f$.
* @param[in] q the azimuthal quantum number of the integral kernel
* @param[out] K the list of \c N eigenvectors (radial functions). 
*        \f$\int_{x=0}^{1/2}xK^2(x)dx=1=\sum_0^{N-1}x K^2(N)/(2N)\f$
* @param[out] lambda the list of \c N squared eigenvalues \f$\lambda_{pq}^2\f$
* @param[in] N the number of points in the interval \f$[0,1/2]\f$ on which the
*       integral kernel and eigenfunctions are discretized.
*/
void KarhLoeve::Kpq(const int q, double *K, double *lambda, const int N)
{
	/** matr contains the N by N matrix of the integral kernel.
	*/
	double *matr = new double[N*N] ;
	for(int i=0; i< N; i++)
	{
		/** The abscissa points on the radial axis are N equidistant points
		* between 0 and 1/2.
		*/
		const double x=(i+0.5)/(2.*N) ;
		for(int j=0; j<=i ; j++)
		{
			double xprime=(j+0.5)/(2.*N) ;

			/** The values of \f$\tilde R_q(x,x')\f$ are collected by a
			* call to Rsymm(). A factor \f$1/(2N)\f$ is introduced;
			* in this finite element representation it represents the stride on the
			* abscissa in the Riemann sense of integration.
			*/
			matr[i+j*N]=matr[j+N*i] = sqrt(x*xprime)*Rsymm(x,xprime,q)/(2.*N) ;
		}
	}

	/** Some intermediate variables are set up since we will be using the
	* standard LAPACK Fortran interface.
	*/
	char jobz[2] = "V" ;
	char uplo[2] = "U" ;
	integer lda = N ;
	integer lwork = 3*N ;
	integer info ;
	double *work = new double[lwork] ;

	/** The finite element version of the integral equation with symmetric matrix
	* is solved by a call to LAPACK's dsyev().
	*/
	dsyev_(jobz,uplo,&lda,matr,&lda,lambda,work,&lwork,&info) ;
	if ( info)
		cerr << __LINE__ << " dsyev info " << info << endl ;
	else
	{
#if 0
		/* The piston mode never appears here... */
		for(int i=N-1 ; i ; i--)
		{
			if( lambda[i] >= 0.)
			{
				cout << "# p=" << N-1-i << " q=" << q << " lam^2= "<< lambda[i] << endl ;
			}
			else
				cerr << "# " << __LINE__ << " lambda " << i << " " << lambda[i] <<  endl ;
		}
#endif

		// re-normalize such that they are positive at x=1/2
		for(int e=0; e < N ;e++)
		{
			if( matr[N-1+N*e] < 0.)	// last value of the eigenvector
				for(int j=0; j < N ; j++)
					matr[j+N*e] *= -1 ;

		}
		// the dsyev normalization: the sum of these squares is 1 now, still incorporating the symmetrizer,
		// equivalent to sum(j=0..N) x*Kpq(x)^2=1.
		// symmetrize such that integral(x=0..1/2)x*K^2= 1 according to eq 43), = sum(j=0..N-1) x*Kpq(x)^2/(2N)
		// since the step width is 1/(2N) in the Riemann sense for an upper limit of 1/2.
		const double nor = 2.*N ;	// this means multiply K^2 by 2N
				// and implies the angular functions are cos(q*theta)/sqrt(pi),
				// sin(q*theta)/sqrt(pi) for q>=1, or 1/sqrt(2pi) for q=0.

		/* remove the symmetrizer factor, and multiply with the the normalization to have the correct 2N or N sum
		* They are now normalized to sum_{0..N} ( x*K^2* Delta(x)=1 where
		* Delta(x)=0.5/N.
		*/
		for(int e=0; e < N ;e++)
			for(int j=0; j < N ; j++)
			{
				const double x=(j+0.5)/nor ;
				matr[j+N*e] /= sqrt(x/nor) ;
			}

#if 0
		/* debugging */
		cout << "q=" << q << " f matrix lambda^2 " << endl ;
		for(int e=0; e < N ;e++)
		{
			if ( lambda[e] > 1.e-3 )
			{
				cout << "# lambda^2 " << lambda[e] << endl ;
				for(int j=0; j < N ; j++)
				{
					cout << e << " " << (j+0.5)/nor << " " << matr[j+N*e] << endl ;
				}
				cout << endl << endl ;
			}
		}
#endif

#ifdef KARHLOE_RADIAL_PFIT
		double *matrXp = new double[N*N] ;

#if 0

		for(int e=0; e < N ;e++)
			for(int j=0; j < N ; j++)
			{
				const double x=(j+0.5)/nor ;
				//matrXp[j+N*e] = pow(x,0.5*q+2*e) ;
				//matrXp[j+N*e] = pow(x,q+2*e) ;
				/* Zernike decomposition */
				matrXp[j+N*e] = Kq::zbase.at(x,q+2*e,q) ;
			}

		/* debugging */
		cout << "q=" << q << " x power matrix " << endl ;
		for(int e=0; e < N ;e++)
		{
			for(int j=0; j < N ; j++)
			{
				cout << "x_" << j << " ^ " << q+2*e << " " << matrXp[j+N*e] << endl ;
			}
		}
#endif /* 0 */


#if 0
		integer *ipiv = new integer[N] ;

		/** The system of linear equations 
		* is solved by a call to LAPACK's dgesv().
		*/
		dgesv_(&lda,&lda,matrXp,&lda,ipiv,matr,&lda,&info) ;
		if ( info)

			cerr << __LINE__ << " dgesv info " << info << endl ;

		delete ipiv ;

		/* save the eigenvectors into return area
		*/
		memcpy(K,matr,N*N*sizeof(double)) ;
#else /* 0 */

		/* for the radial solutions enumebered by e, expand them individually
		* into the lowest N applicable Zernike polynomials.
		*/
		for(int e=0; e < N ;e++)
		{
			/* just as a check of numerical fitness, accumulate the squared sum of
			* the expansion coefficients, which ought add to unity because the
			* Zernike polynomials as well as the KL eigenfunctions habe been normoalized
			* to unite with the weight x.
			*/
			double rintg2=0. ;
			for(int zc = 0 ; zc < N ; zc++)
			{
				double rintg =0.;
				for(int j=0; j < N ; j++)
				{
					const double x=(j+0.5)/nor ;

					/* a strictly Rieman version that ignores the wiggles in the
					* Zernike polynomials and is therefore worse when we compare the sum of the
					* squared expansion coefficients below with the precise 1.
					* rintg += x* Kq::zbase.at(x,q+2*zc,q) *matr[j+N*e]/nor ;
					*/

					/* Use a step function for x*KL(x), which is the accurate representation
					* of the normalization used for the discrete solution, but within this
					* step use a more accurate integrated version of the Zernike polynomials
					*/
					rintg += x* Kq::zbase.intgrl(q+2*zc,q,j/nor,(j+1)/nor)*matr[j+N*e] ;
				}


				/* debugging
				if ( lambda[e] > 1.e-4 )
					cout << "store q=" << q << " lambda^2 " << lambda[e] << " coef " << zc << " " << rintg << endl ;
				*/
				rintg2 += rintg*rintg ;

				matrXp[zc+N*e] = rintg ;
			}

			/* delete strong wiggles and bad fits */
			if ( rintg2 > 4.5)
			{
			cout << "q=" << q << " lambda^2 " << lambda[e] << " norm " << rintg2 << endl ;
				for(int zc =0 ; zc < N ;zc++)
					matrXp[zc+N*e] = 0. ;
			}

			/* debugging
			* cout << "q=" << q << " lambda^2 " << lambda[e] << " norm " << rintg2 << endl ;
			*/
		}

#if 0
		/* debugging */
		cout << "q=" << q << " reprod " << endl ;
		for(int e=0; e < N ;e++)
		{
			if ( lambda[e] > 1.e-3 )
			{
				cout << "# lambda^2 " << lambda[e] << endl ;
				for(int j=0; j < N ; j++)
				{
					const double x=(j+0.5)/nor ;
					double fat =0. ;
					for(int i=0; i < N ;i++)
						fat += matrXp[i+N*e]*Kq::zbase.at(x,q+2*i,q) ;
					cout << x << " " << fat << endl ;
				}
				cout << endl << endl ;
			}
		}
#endif

		/* save the eigenvectors into return area
		*/
		memcpy(K,matrXp,N*N*sizeof(double)) ;
#endif /* 0 */

#if 0
		/* debugging */
		cout << "Zern coeff exp q=" << q << endl ;
		for(int e=0; e < N ;e++)
		{
			cout << "radial # " << e << " lambda^2 " << lambda[e] <<  endl ;
			double fat =0. ;	// check sumsquared
			for(int j=0; j < N ; j++)
			{
				cout << "x^" << q+2*j << " " << matr[j+N*e] << endl ;
				fat += pow(matr[j+N*e],2)  ;
			}
			cout << "radial # " << e << " coef^2 " << fat <<  endl ;

			/*
			for(int j=0; j < N ; j++)
			{
				const double x=(j+0.5)/nor ;
				double fat =0. ;
				for(int i=0; i < N ;i++)
					fat += matr[i+N*e]*Kq::zbase.at(x,q+2*i,q) ;
				cout << x << " " << fat << endl ;
			}
			*/
		}
#endif /* 0 */

		delete matrXp ;
#else /* KARHLOE_RADIAL_PFIT */

		/* save the eigenvectors into return area
		*/
		memcpy(K,matr,N*N*sizeof(double)) ;

#endif /* KARHLOE_RADIAL_PFIT */
			
	}
	delete work ;
	delete matr ;
}
#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN
/*
B_0, B_1 and B_2 of equation (40)-(42).
Digits := 26 ;
assume(x>0) ;
ikern := (x^2+x2p^2-2*x*x2p*cost)^(5/6) ;
las  := taylor(ikern,cost=0,20) ;
las := convert(las,polynom) ;
cost := cos(t) ;
las := int(las,t=0..2*Pi) ;
b := las ;
las := las*x2p ;
las := int(las,x2p=0..1/2) ;
las := 3.44*4*las/Pi ;
las := expand(las) ;
las := simplify(las,power) ;
las := simplify(las,radical) ;
evalf(las) ;

# instead from x=0..1/2 and x2p=0..1/2 only up to diagonal..
b := 2*b*x*x2p ;
b := int(b,x=0..x2p) ;
b := int(b,x2p=0..1/2) ;
evalf(b*3.44*4/Pi*8) ;

accuracy about 0.002 for ord=8
accuracy about 0.0003 for ord=20
*/
double KarhLoeve::B(const int n, const double x)
{
	switch(n)
	{
		double resul, aux ;
	case 0:
		return FRIED_CONSTANT*pow(x,5./3.) ;
		break ;
	case 1:
		/* int(x''=0..1/2) dx'' (x^2+x''^2)^(5/6) int(0..2pi) dtehta [1-2xx''/()]^(5/6)
		=int(x''=0..1/2) dx'' (x^2+x''^2)^(5/6) int(0..2pi) dtehta sum(l=0,2,4,..)(5/6 chose l)[-2xx''/()]^l Int(..)
		=sum(l=0,2,4,...) (5/6 chose l) (-2x)^l int(x''=0..1/2) x''^(l+1)/(x^2+x''^2)^(l-5/6) cosInt(l)
		where with substitution x''^2=y
		int(x''=0..1/2) x''^(l+1)/(x^2+x''^2)^(l-5/6)dx''=(1/2)int(y=0..1/4) x''^l/(x^2+y)^(l-5/6)
		= (1/2)int(y=0..1/4) y^(l/2)/(x^2+y)^(l-5/6) and with substit z=x^2+y
		= (1/2)int(z=x^2..x^2+1/4) (z-x^2)^(l/2)/z^(l-5/6) and another binomial series
		= (1/2)int(z=x^2..x^2+1/4) sum(m=0..l/2) (l/2 chose m) z^(m-l+5/6) (-x^2)^(l/2-m)
		= (1/2) sum(m=0..l/2) (l/2 chose m) (-x^2)^(l/2-m) [z^(m-l+11/6)/(m-l+11/6)]between x^2+0.25 and x^2
		*/
		resul=0. ;
		for(int l=0; l <= ord ; l += 2)
		{
			double msum=0. ;
			for(int m=0; m <= l/2; m++)
				msum += binom(0.5*l,m)*pow(-x*x,l/2-m)
					*(pow(x*x+0.25,m-l+11./6.)-pow(x*x,m-l+11./6.))/(m-l+11./6.) ;
			resul += 0.5*msum*binom(5./6.,l)*pow(-2.*x,l)*cosInt(l);
		}
		return FRIED_CONSTANT*4./M_PI*resul ;
		break ;
	case 2:
		return B2() ;
		break ;
	default:
		cerr << __LINE__ << " n " << n << " not supported\n" ;
		exit(EXIT_FAILURE) ;
	}
}

#endif /* KARHLOE_RADIAL_POLYN */

/** Binomial coefficients.
* @param[in] x the ``numerator'' of the coefficient.
*    This may be an arbitrary, positive or negative real number.
* @param[in] l the ``denominator'' of the coefficient.
*    This is is limited to the range 0,1,2,...
* @return \f$x \choose l\f$
*/
double KarhLoeve::binom(double x, int l)
{
	double resul =1. ;
	for(int c=1; c <=l ; c++)
		resul *= (x+1.-c)/c ;
	return resul ;
}

#ifndef KARHLOE_RADIAL_POLYN

/** the integral over \f$\cos^l x\f$ over one period.
* @param[in] l the integer power of the integrand
* @return \f$\int_0^{2\pi}\cos^l xdx\f$.
*/
double KarhLoeve::cosInt(int l)
{
	/** The value is zero for odd \f$l\f$.
	*/
	if( l %2 )	/* zero if l is odd=1,3,5,-1,-3,-5,.. */
		return 0 ;
	double resul ;
	switch(l)
	{
	case 0:
		/** The value is \f$2\pi\f$ for \f$l=0\f$.
		*/
		resul = 2.*M_PI ;	// int(t=0..2*pi,1)=2*pi 
		break ;
	default :
		resul=M_PI ;
		l /= 2 ;
		/** In the general case of even \f$l\f$ we use
		* \latexonly
		* \cite[2.512.2]{GR}.
		* \endlatexonly
		*/
		// (oldl-1)!!/2^newl/newl!*2*pi; no loop for oldl=2, newl=1
		for(int c=2 ; c <= l; c++)
			resul *= 0.5*(2.*c-1.)/(double)c ;
	}
	return resul ;
}

/** the integral over \f$\cos^l x\cos(qx)\f$ over one period.
* @param[in] l the integer power of the integrand
* @param[in] q the second integer of the integrand
* @return \f$\int_0^{2\pi}\cos^l x\cos (qx)dx\f$.
*/
double KarhLoeve::cosInt(const int l, int q)
{
	q = abs(q) ;
	switch(q)
	{
		double resul ;
	case 0 :
		/** the cases \f$q=0,1\f$ are forwarded to cosInt(const int l) .
		*/
		return cosInt(l) ;
		break ;
	case 1 :
		return cosInt(l+1) ;
		break ;
	default :
		/*
		int (t=0..2*pi, cos^l t *cos(qt)); q=+-1, +-2,..GR 2.512.2.
		with cos^lt cos(qt)=cos^(l-2)t cos^2t cos(qt)
		= (1/2)cos^(l-2)t [1+cos(2t)] cos(qt)
		= (1/2)cos^(l-2)t cos(qt) + (1/2) cos^(l-2)t cos(2t)cos(qt)
		= (1/2)cos^(l-2)t cos(qt) + (1/4) cos^(l-2)t {cos[(2+q)t]+cos[(2-q)t]}
		or via binomial expansion of cos^l t = [(e^i*t+e^-i*t/2]^l=sum(k=0..l)binom(l,k)..
		*/
		switch(l)
		{
		case 0 :
			return 0. ;
			break ;
		case 1:
			// cos(t)*cos(qt)=(1/2)[cos(1+q)t+cos(1-q)t] only result if 1+q=0 or 1-q=0: (1/2)*2*pi
			return ( q == 1 ) ? M_PI : 0. ;
			break ;
		default:
			resul=0. ;
			// if ( l-q >=0 && l-q <= 2*l) superfluous the l-q<=2*l since -q<=l always for q>=0
			if ( l-q >=0 )
			{
				if( (l-q) %2 == 0 )	// would meet 2k
					resul += binom( (double)l, (l-q)/2 ) ;
			}
			//if ( l+q >=0 && l+q <= 2*l) superfluous the l+q>=0 since q>0 and l>0
			if ( l+q <= 2*l)
			{
				if( (l+q) %2 == 0 )	// would meet 2k
					resul += binom( (double)l, (l+q)/2 ) ;
			}
			return resul*M_PI/pow(2.,l) ;
		}
	}
}

#endif /* KARHLOE_RADIAL_POLYN */

#if 0
/**
* @param[in] j the common integer of the powers of the integral, non-negative integer.
* @param[in] q the azimuthal quantum number of the mode, non-negative
* @param[in] o half of the effective order of the hypergeometric function, non-negative
* @return the integral of \f$\sin^{j+1+2(q+o)}\alpha/\cos^{j+14/3}\alpha\f$
*   integrated from \f$\alpha=0\f$ to \f$\alpha=\pi/4.\f$
* @since 2007-03-19
*/
double KarhLoeve::I2(const int j, const int q, const int o)
{
	if ( j % 2 )
	{
		/* case of odd j, used with recursion
		*/

		if ( o > 0 )
		{
			double resul=0. ;

			/* naive, non-recursive implementation of the binomial ..
			*/
			for(int l=0 ; l <= o ; l++)
				if ( l % 2 )
					resul -= binom(o,l)*I2(j,q+o+l,0) ;
				else
					resul += binom(o,l)*I2(j,q+o+l,0) ;

			/* debugging
			*/
			cout << __FILE__ << " " << __LINE__ << " " << j << " " << q << " " << o << " I2 " <<  resul << endl ;
			return pow(4.,o)*resul ;
		}
		else
		{
			if ( q > 0 )
			{
				return ((j+2*q)*I2(j,q-1,0)-pow(2.,11./6.-q))/(2*q-11./3.) ;
			}
			else
			{
				if ( j == 1)
					/* here qo=0, int(0..pi/4) sin^2 alpha / cos^17/3 alpha	
					*/
					return 0.492200573702716302791048473793 ;
				else
					return (M_2_POW_11_6-j*I2(j-2,0,0))/(j+11./3.) ;
			}
		}
	}
	else
	{
#ifdef HAVE_GSL
		/* 2F1(-j/2-q-o, -j/2-11/6+o;-5/6-j/2+o;x^2) at the upper limit x=1 
		* using
		* @latexonly
		* \cite[15.1.20]{AS}
		* @endlatexonly
		*/

		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " gamma " << -5./6.-0.5*j+o << " " 
		* 	<< 1.+0.5*j+q+o << " " << -5./6.+q+2*o << endl ;
		*/

		const double ulim = gsl_sf_gamma(-5./6.-0.5*j+o) *gsl_sf_gamma(1.+0.5*j+q+o) /gsl_sf_gamma(-5./6.+q+2*o);

		/* debugging
		*	cout << __FILE__ << " " << __LINE__ << " ulim " << ulim << endl ;
		* cout << __FILE__ << " " << __LINE__ << " 2f1 " << -0.5*j-q-o << " " 
		* << -11./6.-0.5*j+o << " " << -5./6.-0.5*j+o << endl ;
		*/

		/* 2F1(-j/2-q-o, -j/2-11/6+o;-5/6-j/2+o;x^2) at the lower limit x=1/sqrt(2).
		* If the first parameter is zero, the GSL implementation returns an error, which we
		* avoid by knowing that this case means 2F1()=1.
		* If the first parameter is a negative integer smaller than -10 (larger in absolute terms than 10),
		* the gsl interface will also report an error. If this is a negative integer >= -9, it will
		* use hyperg_2F1_series(). See gsl_sf_hyperg_2F1_e() which is declared static in the C source and not
		* available through the library interface.
		*/
		const int jhalf = j/2 ;

		double bin = 1.0;

		/* debuging
		* double llim = (0.5*j+q+o > 0.2) ? gsl_sf_hyperg_2F1(-0.5*j-q-o,-11./6.-0.5*j+o,-5./6.-0.5*j+o,0.5) : 1.0 ;
		*/
		double llim = 1.0 ;

		/* upper limit is jhalf+q+o but this contributes with a term of zero. So using < instead <= is sufficient
		*/
		for( int sigm = 1 ; sigm < jhalf+q+o; sigm++)
		{
			/* update the binomial coefficient, including the (-)^sigm
			*/
			bin *= (jhalf+q+o-sigm+1)/(-11./6.-jhalf+o+sigm) ;
			llim += bin ;
		}

		llim /= pow(2.,-11./6.+q+2*o) ;
	cout << __FILE__ << " " << __LINE__ << " llim " << llim << endl ;
	cout << __FILE__ << " " << __LINE__ << " j " << j << " q " << q << " " << o << " I2 " <<  pow(4.,o)*(ulim-llim)/(-j-11./3.+2*o) << endl ;
		return pow(4.,o)*(ulim-llim)/(-j-11./3.+2*o) ;
#else /* HAVE_GSL */
		/* case of even j gives a finite sum
		*/
		const int jhalf = j/2 ;

		/* binomial(j/2+q+o,sigm), initialized at sigm=0
		*/
		double bin = 1.0;

		/* 1/2^(sigm+o-j/2-11/6), initialized at sigm=0
		*/
		double pow2 = pow(2.,jhalf-o+11./6.) ;
		int sigm =0 ;

		/* the resulting value of the integral, initialized at sigm=0
		*/
		double resul= (1.-pow2)/(2*o-j-11./3.) ;
		for( sigm = 1 ; sigm <= jhalf+q+o; sigm++)
		{
			/* update the component 1/2^sigm in the power of 2
			*/
			pow2 /= 2. ;

			/* update the binomial coefficient, including the (-)^sigm
			*/
			bin *= -(double)(jhalf+q+o+1-sigm)/sigm ;
			resul += bin*(1.-pow2)/(2*(sigm+o) -j-11./3.) ;
		}

		/* debugging
		* cout << __FILE__ << " " << __LINE__ << " j " << j << " q " << q << " " << o << " I2 " <<  resul << endl ;
		*/
		return pow(4.,o)*resul ;
#endif /* HAVE_GSL */
	}
}
#endif /* 0 */

#ifdef KARHLOE_RADIAL_POLYN

/**
* @param[in] j power of the basis of the bra
* @param[in] s power of the basis of the ket
* @param[in] q azimuthal parameter of the KL mode
* @return the integral of \f$\int_0^{2\pi)d\theta \int_0^{\pi/4}d\alpha \sin^{1+q}(2\alpha)
*   [1-\sin(2\alpha)\cos\theta]^{5/6}cos(q\theta)/(1+2q+j+s+14/3)/(2\cos\alpha)^{1+j+q+2q+14/3}\f$
* @since 2007-03-26
* @warning the factor -3.44 or 2pi or (-5/6)_q/q! is not yet in here.
* @warning the overlap integral is for the symmetrized kernel.
* @warning the overlap integral is incomplete for q=0.
* @warning the index offset of the power basis is already worked into the response.
*   The parameter j might already have been reduced by q.
*/
double KarhLoeve::Ialpha(int j, int s, const int q)
{

	if ( j < s )
	{
		const int tmp=j ; j=s ; s=tmp ;
	}
		

	/* store[q][j][s] */
	static vector< vector< vector<double> > > store ;
	static int oldq = -1 ;
	static vector<double> f21store ;

	if ( oldq != q )
		f21store.clear() ;

	/* debuging
	* cout << __FILE__ << " " << __LINE__ <<  " " << j << " " << s << " " << q << endl ;
	*/

	/* supposed we need the value of q=2, stored in store[2], we fill store[0] and store[1] if needed,
	* at the end of which we'll have store.size()=q=2.
	*/
	while ( q >= store.size())
	{
		vector< vector<double> > tmp ;
		store.push_back(tmp) ;
	}

	while ( j >= store[q].size() )
	{
		vector<double> tmp ;
		store[q].push_back(tmp) ;
	}


	/* Two possibilities here, s = store[q][j].size() or s < store[q][j].size().
	* Add/compute store[q][j][s] if not yet done.
	*/
	while ( s >= store[q][j].size() )
	{
		const int sloc = store[q][j].size() ;

		double resul=0. ;

		/* alp=0 to pi/4 of
		* [sin^j(alp)/cos^(s+14/3)(alp)+(j<->s)] sin^(1+2q)(alp) 2F1[q/2-5/12,q/2+1/12;q+1;sin^2(2alp)]
		*/
		int N= 4000 ;

#if 1
		/* Simpson */
		for(int i=0 ; i <= N ; i++)
		{
			/* weights 1,4,2,4,2,...,4,1 */
			int wei = (i==0 || i==N)? 1 : ( (i %2)?4 :2 ) ;
			const double alp = i*M_PI_4/(double)N ;
			const double salp = sin(alp) ;
			const double calp = cos(alp) ;

			/* avoid domain error with sin(2*alpha)=1 at 2*alpha=pi/2 or alpha=pi/4.
			*/
			double f21  ;

			if ( oldq == q)
				f21 = f21store[i] ;
			else
			{
#if 0
				f21 = ( i==N) ? 
				0.897686900876088227830096388406*pow(2.,q)* gsl_sf_fact(q) / gsl_sf_gamma(q+11./6.)
				: gsl_sf_hyperg_2F1(0.5*q-5./12.,0.5*q+1./12.,q+1.,pow(2.*salp*calp,2.)) ;
#else
				if ( i == N )
					f21 = 0.897686900876088227830096388406*pow(2.,q)* gsl_sf_fact(q) / gsl_sf_gamma(q+11./6.) ;
				else
				{
					const double a= 0.5*q-5./12. ;
					const double b= 0.5*q+5./12. ;
					const double c=q+1. ;
					double z= 2.*salp*calp ;	/* my z, argument is z^2 */
					/* AS 15.3.20 */
					//f21 = pow(1+z,-2.*a)*gsl_sf_hyperg_2F1(2.*a,c-1./2.,2.*c-1.,2.*z/(1.+z)) ;
					/* AS 15.3.21 */
					//f21 = pow(1-z*z,-a)*gsl_sf_hyperg_2F1(2.*a,2.*(c-a)-1.,c,0.5*(sqrt(1.-z*z)-1.)/sqrt(1.-z*z)) ;
					/* AS 15.3.19 */
					z = sqrt(1.-z*z) ;
					f21 = pow(0.5+0.5*z,-2.*a)*gsl_sf_hyperg_2F1(2.*a,2.*a-c+1.,c,(1.-z)/(1.+z)) ;
				}
#endif
				f21store.push_back(f21) ;
			}

			resul += wei*( pow(salp,j)/pow(calp,j+14./3.)+pow(salp,sloc)/pow(calp,sloc+14./3.) )
				*pow(salp,1+2*q) *f21 ;
		}


		/* multiply by Simpson weight, which is a third of the step width \f$\pi(4N)\f$ */
		resul *= M_PI_4/(3*N) ;
#else
		/* Simpson */
		for(int i=0 ; i <= N ; i++)
		{
			/* weights 1,4,2,4,2,...,4,1 */
			int wei = (i==0 || i==N)? 1 : ( (i %2)?4 :2 ) ;
			const double x = M_SQRT1_2+i*(1.-M_SQRT1_2)/(double)N ;
			const double salp = sqrt(1.-x*x) ;
			const double calp = x;

			/* avoid domain error with sin(2*alpha)=1 at 2*alpha=pi/2 or alpha=pi/4.
			*/
			double f21  ;

			if ( oldq == q)
				f21 = f21store[i] ;
			else
			{
				f21 = ( i==0) ? 
				0.897686900876088227830096388406*pow(2.,q)* gsl_sf_fact(q) / gsl_sf_gamma(q+11./6.)
				: gsl_sf_hyperg_2F1(0.5*q-5./12.,0.5*q+1./12.,q+1.,pow(2.*salp*calp,2.)) ;
				f21store.push_back(f21) ;
			}

			resul += wei*( pow(salp,j)/pow(calp,j+14./3.)+pow(salp,sloc)/pow(calp,sloc+14./3.) )
				*pow(salp,2*q) *f21 ;
		}


		/* multiply by Simpson weight, which is a third of the step width \f$\pi(4N)\f$ */
		resul *= (1.-M_SQRT1_2)/(3*N) ;
#endif /* 0 */

		const double qjs17 = 2*q+j+sloc+17./3. ;

		store[q][j].push_back( resul/( qjs17* pow(2.,qjs17) ) ) ;
	}
	oldq =q ;
	return store[q][j][s] ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#if 0

/**
* @param[in] j power of the basis of the bra
* @param[in] s power of the basis of the ket
* @param[in] q azimuthal parameter of the KL mode
* @return the integral of \f$\int_0^{2\pi)d\theta \int_0^{\pi/4}d\alpha \sin^{1+q}(2\alpha)
*   [1-\sin(2\alpha)\cos\theta]^{5/6}cos(q\theta)/(1+2q+j+s+14/3)/(2\cos\alpha)^{1+j+q+2q+14/3}\f$
* @since 2007-03-19
* @warning the factor -3.44 or 2pi or (-5/6)_q/q! is not yet in here.
* @warning the overlap integral is for the symmetrized kernel.
* @warning the overlap integral is incomplete for q=0.
* @warning the index offset of the power basis is already worked into the response.
*   The parameter j might already have been reduced by q.
*/
double KarhLoeve::I1(const int j, const int s, const int q, const int o)
{
	const double qjs17 = 2*q+j+s+17./3. ;
	return ( I2(j,q,o)+I2(s,q,o) )/( qjs17* pow(2.,qjs17) ) ;
}

#endif /* 0 */

#ifdef KARHLOE_RADIAL_POLYN

/**
* @param[in] j power of the basis of the bra
* @param[in] s power of the basis of the ket
* @param[in] q azimuthal parameter of the KL mode
* @return the integral of \f$\int_0^{2\pi)d\theta \int_0^{\pi/4}d\alpha \sin^{1+q}(2\alpha)
*   [1-\sin(2\alpha)\cos\theta]^{5/6}cos(q\theta)/(1+2q+j+s+14/3)/(2\cos\alpha)^{1+j+q+2q+14/3}\f$
* @since 2007-03-19
* @warning the factor -3.44 or 2pi or (-5/6)_q/q! is not yet in here.
* @warning the overlap integral is for the symmetrized kernel.
* @warning the overlap integral is incomplete for q=0.
* @warning the index offset of the power basis is already worked into the response.
*   The parameter j might already have been reduced by q.
*/
double KarhLoeve::S(const int j, const int s, const int q)
{
	/* write the theta integral
	* sum(o=0..infinity) binomial(5/6,o) (-z)^o int_0^2pi cos^o theta cos(q theta)d theta
	*/

#if 1
	double resul= Ialpha(j,s,q) ;

	if ( q == 0)
	{
		/* contribution from G1, dropping the 2pi, the 3.44 and
		* switching the sign
		*/
		resul -= Ialpha(0,s,q)/((2+j)*pow(2.,j)) +Ialpha(j,0,q)/((2+s)*pow(2.,s)) ;
	}

#else
	double bin= 1. ;
	double resul= 0. ;
	double oresul= 0. ;

	/* todo: use member parameter ord to limit the loop
	*/
	for( int o=0; o<10 || fabs(resul-oresul) > 1.e-5*fabs(resul) ; o++)
	{
		oresul = resul ;
		resul += bin * I1(j,s,q,o) ;
		if ( q == 0)
		{
			/* contribution from G1, dropping the 2pi, the 3.44 and
			* switching the sign.
			*/
			resul -= I1(0,s,0,o)/(2+j)/pow(2.,(double)j) ;
			resul -= I1(j,0,0,o)/(2+s)/pow(2.,(double)s) ;
		}
		/* the sign is from the term (-sin(2*alpha)*cos theta)^o
		* in the binomial expansion of the theta-integral.
		* Update (q/2-5/12)_o (q/2+1/12)_0 / (q+1)_o/o! to the next o.
		* (q/2-5/12+o-1)-> (q/2-5/12+o)
		* (q/2+1/12+o-1)-> (q/2+1/12+o)
		* (q+1+o-1)-> (q+o+1)
		* o-> (o+1)
		*/
		bin *= -(0.5*q+o-5./12.)*(0.5*q+o+1./12.)/((q+o+1)*(o+1)) ;
	}
#endif

	if ( q == 0)
	{
		/* contribution from G2, dropping the 2pi, the 3.44 and
		* the negative sign.
		*/
		resul += 0.29995353261705446/((2+j)*(2+s)*pow(2.,4+j+s)) ;
	}
	
	/* optionally include the important factors left out up to here
	*  resul *= -FRIED_PARAMETER*2.0*M_PI*Pochhammer(-5/6,q)/factorial(q)
	*/
	return resul ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#ifdef KARHLOE_RADIAL_POLYN

/**
* @param[in] j index to the Zernike polynomials, 0-based
* @param[in] s index to the Zernike polynomials, 0-based
* @param[in] q azimuthal parameter of the KL mode
* @param[in] z the generic global radial Zernike functions
* @return the integral of \f$\int_0^{2\pi)d\theta \int_0^{\pi/4}d\alpha \sin^{1+q}(2\alpha)
*   [1-\sin(2\alpha)\cos\theta]^{5/6}cos(q\theta)/(1+2q+j+s+14/3)/(2\cos\alpha)^{1+j+q+2q+14/3}\f$
* @since 2007-03-19
* @warning the factor -3.44 or 2pi or (-5/6)_q/q! is not yet in here.
* @warning the overlap integral is for the symmetrized kernel.
* @warning the overlap integral is incomplete for q=0.
* @warning the index offset of the power basis is already worked into the response.
*   The parameter j might already have been reduced by q.
*/
double KarhLoeve::S(const int j, const int s, const int q, ZernikeBase & z)
{
	double resul= 0. ;

	/* The values of Noll's n parameter are starting at q and increasing to
	* whatever we need to have j (s) independent degrees of freedom
	*/
	const int nj = q+2*j ;
	const int ns = q+2*s ;

	/* the j and s terms are polynomials in z.Z[j] and z.Z[s], implicitly those
	* where the azimuthal parameter is q. This function here is called always in a context of
	* fixed q, so we may assume that the indices j and s are 0 if the nodal parameter is n=q,
	* that they are 1 if it is n=q+2 etc. The powers of the x^i are q(=m)<=i <= n, and we
	* build the distributed product over the two polynomials.
	*/
	for(int powj = q ; powj <= nj ; powj +=2)
		for(int pows = q ; pows <= ns ; pows +=2)
		{
			/* The interface to Zernike::powcoef() demands the Noll parameter n,
			* the Noll parameter m, and the nominal exponent of the radial power.
			* The interface to the reduces KarhLoeve::S() uses the powers that
			* are left over after subtracting the m parameter.
			*/
			resul += z.powcoef(nj,q,powj)* z.powcoef(ns,q,pows)* S(powj-q,pows-q,q) ;
		}
	
	/* optionally include the important factors left out up to here
	*  resul *= -FRIED_PARAMETER*2.0*M_PI*Pochhammer(-5/6,q)/factorial(q)
	*/
	return resul ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN

/** The auxiliary integral
* \latexonly
* \cite[(42)]{WangJOSA68}
* \endlatexonly
* @return \f$\int_0^{1/2}xdx\int_0^{1/2}x'dx' \int_0^{2\pi}d\theta (x^2+x'^2-2xx'\cos\theta)^{5/6}\f$.
*/
double KarhLoeve::B2()
{
	double resul =0. ;
	/*
	int(x=0..1/,x''=0..1/2) x x'' int(theta=0..2pi) [x^2+x'^2-2xx''cos(theta)]^(5/6)
	 the binomial expansion after breaking (x^2+x''^2) away from the square bracket
	=int(x,x'') x x'' (x^2+x''^2)^(5/6) sum(l=0..infity,here:ord) (5/6 choose l) [-2xx''/(x^2+x''^2)]^l int(theta)cos^l(theta)
	=int(x,x'') x x'' (x^2+x''^2)^(5/6) sum(l=0..infity,here:ord) (5/6 choose l) [-2xx''/(x^2+x''^2)]^l cosInt(l)
	=sum(l=0..infity,here:ord) (5/6 choose l) B2int(l) cosInt(l)
	*/

	/** The numerics uses a binomial expansion up to order \c ord . The number
	* of terms used is only half as large because only terms of even index contribute to the sum.
	*/
	for(int l=0; l <= ord ; l += 2)
		resul += B2int(l)*cosInt(l)*binom(5./6.,l) ;
	/** The result is correctly multiplied with the constants of
	* \latexonly
	* \cite[(41)]{WangJOSA68} and \cite[(42)]{WangJOSA68}.
	* \endlatexonly
	*/
	return resul*FRIED_CONSTANT*4/M_PI*8. ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN
/*****************************
case q=0:

Digits := 26 ;
assume(x>0) ;
ikern := (x^2+x2p^2-2*x*x2p*cost)^(5/6) ;
las  := taylor(ikern,cost=0,8) ;
las := convert(las,polynom) ;
cost := cos(t) ;
las := int(las,t=0..2*Pi) ;
las := 3.44*las ;
las := expand(las) ;
las := simplify(las,power) ;
las := simplify(las,radical) ;
evalf(las) ;

Rsymm return R(x,x')/x' to support the symmetrized version
******************************/
double KarhLoeve::Rsymm(double x, double xprime, const int q)
{
	double resul;
	const double xoverxprime = -2.*x*xprime/(x*x+xprime*xprime) ;
	switch(q)
	{
	case 0:
		if( x==0. || xprime==0.)
			resul = FRIED_CONSTANT*2.*M_PI*pow(x*x+xprime*xprime,5./6.) ;
		else	// 3.44*int(0..2*pi, [x^2+x'^2-2*x*x'*cos(t)]^(5/6), binomial up to 'ord'
			// 3.44*(x^2+x'^2)^(5/6) int(0..2*pi( 1-xoverprime*cos(t))^5/6)
			// =3.44*(x^2+x'^2)^(5/6) sum(l=0..ord) (-xoverprime^l binom(5/6,l)* int(t=0..2pi,cos^l(t))
		{
			resul = cosInt(0) ;	// term l=0
			//for(int l=1; l <= ord ; l++)	// skip the l=1,3,5,7.. because then the cosInt()=0.
			for(int l=2; l <= ord ; l += 2)	// skip the l=1,3,5,7.. because then the cosInt()=0.
				resul += pow(xoverxprime,l)*cosInt(l)*binom(5./6.,l) ;
			resul *= FRIED_CONSTANT*pow(x*x+xprime*xprime,5./6.) ;
		}
		return -resul+2.*M_PI*(B(1,x)+B(1,xprime)-B(2,x)) ;
		break ;
	default:
		if( x==0. || xprime==0.)
			return 0. ; // then the integral over t=0..2*pi of cos(q*t) is zero.
		else	// 3.44*int(0..2*pi, [x^2+x'^2-2*x*x'*cos(t)]^(5/6)*cos(q*t), binomial up to 'ord'
			// 3.44*(x^2+x'^2)^(5/6) int(0..2*pi( 1-xoverprime*cos(t))^5/6)
			// =3.44*(x^2+x'^2)^(5/6) sum(l=0..ord) (-xoverprime^l binom(5/6,l)* int(t=0..2pi,cos^l(t)*cos(q*t))
		{
			resul = cosInt(0,q) ;	// term l=0
			for(int l=1; l <= ord ; l++)
				resul += pow(xoverxprime,l)*cosInt(l,q)*binom(5./6.,l) ;
			resul *= FRIED_CONSTANT*pow(x*x+xprime*xprime,5./6.) ;
		}
		return -resul ;
		break ;
	}
}
#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN

/******************
I := int(x=0..1/2,x''=0..1/2) x x'' (x^2+x''^2)^(5/6)[-2xx''/(x^2+x''^2)]^l
Substitute x=r cos(phi), x''=r sin(phi); because integrand is symmetric, integrate only
up to phi=0..Pi/4, and r=0..1/(2 cos phi) and multiply by 2. Jacobian is r.
I := 2*(-2)^l int(phi=0..Pi/4,r=0..1/(2 cos phi)) r* x^(l+1) x''^(l+1) (x^2+x''^2)^(5/6)[1/(x^2+x''^2)]^l
   = 2*(-2)^l int(phi=0..Pi/4,r=0..1/(2 cos phi)) r* x^(l+1) x''^(l+1) (x^2+x''^2)^(5/6-l)
   = 2*(-2)^l int(phi=0..Pi/4,r=0..1/(2 cos phi)) r^(2l+3) (x^2+x''^2)^(5/6-l) cos^(l+1)phi sin^(l+1)phi
   = 2*(-2)^l int(phi=0..Pi/4,r=0..1/(2 cos phi)) r^(2l+3) r^(5/3-2l) cos^(l+1)phi sin^(l+1)phi
   = 2*(-2)^l int(phi=0..Pi/4,r=0..1/(2 cos phi)) r^(4+2/3) cos^(l+1)phi sin^(l+1)phi
	underivative r^(5+2/3)/(17/3)=r^(17/3)/(17/3)
   = 2*(-2)^l/2^(17/3)*3/17 int(phi=0..Pi/4) cos^(l+1)phi sin^(l+1)phi/cos^(17/3)phi
	subsittute x=cos phi, dx =- dphi sin phi
   = 3/17*(-2)^l/2^(14/3) int(x=1..1/sqrt(2)) (-) x^(l+1) sin(phi)^l /x^(17/3)
   = 3/17*(-2)^l/2^(14/3) int(x=1/sqrt(2)..1) x^(l-14/3) sin(phi)^l
   = 3/17*(-2)^l/2^(14/3) int(x=1/sqrt(2)..1) x^(l-14/3) [1-x^2]^(l/2) (only needed for even l)
   = 3/17*(-)^l 2^(l-14/3) int(x=1/sqrt(2)..1) x^(l-14/3) sum(m=0..l/2) (l/2 choose m) (-x^2)^m
   = 3/17*(-)^l 2^(l-14/3) sum(m=0..l/2) (l/2 choose m) (-)^m int(x=1/sqrt(2)..1) x^(2m+l-14/3)
   = 3/17*(-)^l 2^(l-14/3) sum(m=0..l/2) (l/2 choose m) (-)^m x^(2m+l-11/3)/(2m+l-11/3) at 1 minus 1/sqrt(2)
*************************/
double KarhLoeve::B2int(const int l)
{
	if ( l < 0 || l % 2)
	{
		cerr << __LINE__ << " invalid choice l " << l << endl ;
		exit(EXIT_FAILURE) ;
	}
	double resul =0. ;
	for(int m=0; m <= l/2 ; m++)
		if ( m % 2)
			resul -= binom(0.5*l,m)*(B2intaux(2*m+l,1.)-B2intaux(2*m+l,M_SQRT1_2)) ;
		else
			resul += binom(0.5*l,m)*(B2intaux(2*m+l,1.)-B2intaux(2*m+l,M_SQRT1_2)) ;
	// (-)^l always +1 since l even
	return 3./17.*resul*pow(2.,l-14./3.) ;
}
double KarhLoeve::B2intaux(const int twoml, const double x)
{
	return pow(x,twoml-11./3.)/(twoml-11./3.) ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#undef M_2_POW_11_6
#undef FRIED_CONSTANT

/** Default ctor.
* The internal parameters are instantiated with invalid values.
* This only serves to have a default ctor for use in vector initializations.
*/
Kq::Kq() : q(99), p(-1), lambdasq(-1.), see(1)
{
}

/** Explicit ctor.
* @param[in] qval the azimuthal quantum number \f$q\f$
* @param[in] lambda the square of the KH eigenvalue
* @param[in] K a vector of the values on the radial mesh, with \f$N\f$ elements
* @param[in] N the length of the vector \c K
* @param[in] pval an initializer to indicate the index \c p
* @param[in] rseed an explicit seed value for the random number generator
*/
Kq::Kq(const int qval, const double lambda, double * K , const int N, const int pval, const int rseed ) 
	: q(qval), lambdasq(lambda), p(pval), see(rseed)
{
	for(int e=0; e< N; e++)
		kvect.push_back(K[e]) ;
}

/** Calculate the modal eigenfunction at a point on the screen given in terms of 2 polar coordinates.
* @param[in] usecos flags whether the azimuthal functions \f$\cos q\theta\f$ are called (\c true) or \f$\sin q\theta\f$
* @param[in] r the radial coordinate in the interval \f$[0,1/2]\f$.
* @param[in] theta the azimuthal coordinate [rad] .
* @return the value of the normalized mode, which is the product of the radial and azimuthal functions
*/
double Kq::K(const bool usecos, const double r, const double theta) const
{
	/** Outside the circular pupil return zero.
	*/
	if( r > 0.5)
		return 0. ;

	/** The normalized azimuthal functions are \f$1/\sqrt{2\pi}\f$ (\f$q=0\f$),
	* \f$\cos (q\theta)/\sqrt\pi\f$ and
	* \f$\sin (q\theta)/\sqrt\pi\f$  (\f$q\neq 0\f$).
	*/
	double resul = q ? ( usecos ? cos(q*theta)/M_SQRT_PI : sin(q*theta)/M_SQRT_PI ) : 1.0/(M_SQRT_PI*M_SQRT2) ;

#if defined KARHLOE_RADIAL_POLYN || defined  KARHLOE_RADIAL_PFIT
	resul *= zbase.at(r,kvect,q) ;

	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " " << r << " " << q << " " << resul << endl ;
	*/
#else
			// linear interpolation radial mesh
			// was r=(rindx+0.5)/(2*N), N=kvect.size(), in the middle
	int rindx = 2.*kvect.size()*r	; // round to zero
	//resul *= kvect[rindx]+ (2.*kvect.size()*r-rindx)*(kvect[rindx+1]-kvect[rindx]) ;

	/* no attempt is made to interpolate; just a rounding to the nearest pixel.
	* The somewhat complicated call of Kq::at() ensures that the case of scaling with
	* lambdasq (enabled by the cpp macro KARHLOE_RETURNWI_LAMBDA) works in all cases.
	*/
	resul *= at(rindx) ;
#endif
	return resul ;
}

/** Calculate the modal eigenfunction at a point on the screen given in terms of 2 Cartesian coordinates.
* @param[in] usecos flags whether the azimuthal functions \f$\cos q\theta\f$ are called (\c true) or \f$\sin q\theta\f$
*   (\c false)
* @param[in] x the cartesian coordinates on the screen, each in \f$[-1/2,1/2]\f$
* @return the value of the mode
*/
double Kq::K(const bool usecos, double x[2]) const
{
	/** The implementation decomposes the cartesian coordinates in polar/circular coordinates
	* and forwards the request to K(bool,const double,const double).
	*/
	const double r=hypot(x[0],x[1]) ;
	const double thet = atan2(x[1],x[0]) ;
	return K(usecos,r,thet) ;
}

#ifndef KARHLOE_RADIAL_POLYN

/** Calculate the modal eigenfunction at a point on the screen given in terms of the abscissa index.
* @param[in] idx the index in the range of 0 to N() .
* @return the value of the mode
* @author Richard J. Mathar
* @since 2006-11-29
*/
double Kq::at(const int idx) const
{
#ifdef DEBUG
	cout << __LINE__ << " " << idx << " " << kvect.size() << " " << kvect[idx] << endl ;
#endif
	if ( idx>=0 && idx < kvect.size())
#ifdef KARHLOE_RETURNWI_LAMBDA
		if ( lambdasq > 0.)
			return kvect[idx]*sqrt(lambdasq) ;
		else
			return 0. ;
#else
		return kvect[idx] ;
#endif
	else
		return 0. ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN

/** Calculate the Hankel transform (that is the FT of the basis function).
* The factor of the eigenvalue is not incorporated.
* @param[in] sigma the wavenumber
* @return the value of the radial integral, skipping the factor i^q.
* @author Richard J. Mathar
* @since 2007-04-29
*/
double Kq::fft(const double sigma) const
{
	double resul = 0. ;
	const double h = 0.5/N() ;
	for(int i=0; i < kvect.size(); i++)
	{
		const double x = h*(i+0.5) ;
#ifdef KARHLOE_RETURNWI_LAMBDA
		resul += x*kvect[i]* jn(q,2.*M_PI*sigma*x) ;
#else
		resul += x*at(i)* jn(q,2.*M_PI*sigma*x) ;
#endif
	}

	/** The approximation is the simple Riemann sum with interval lengths \f$\Delta x=0.5/N()\f$,
	* which yields a division through 2*N(). Another factor \f$ 2\pi \f$ comes from the integration over
	* the azimuths.
	*/
	return 2.*M_PI*h*resul ;
}

#endif /* KARHLOE_RADIAL_POLYN */

#ifndef KARHLOE_RADIAL_POLYN

/** Calculate the derivative of the modal eigenfunction at a point on the screen given in terms of the abscissa index.
* @param[in] idx the index in the range of 0 to N() .
* @return the value of the derivative \f$\partial_xK_p^q(x)\f$
* @author Richard J. Mathar
* @since 2006-11-29
*/
inline double Kq::gradat(const int idx) const
{
	if ( idx< 0 || idx >= kvect.size())
		return 0. ;
	else
	{
		/** Step size \c h is the total length of the radial interval, 1/2, divided
		* by the number of abscissa points.
		*/
		const double h = 0.5/N() ;
		/** At the end points we use  the two point approximation
		* \latexonly
		* \cite[25.3.21]{AS}
		* \endlatexonly
		* else the three-point approximation
		* \latexonly
		* \cite[25.3.4]{AS}
		* \endlatexonly
		*/
		if ( idx == 0 )
			return ( at(1)-at(0) )/h ;
		else if ( idx == kvect.size()-1 )
			return ( at(idx)-at(idx-1) )/h ;
		else
			return 0.5*( at(idx+1)-at(idx-1) )/h ;
	}
}

#endif /* KARHLOE_RADIAL_POLYN */

/** The overlap integral with another radial function.
* @param[in] oth the other radial function in the overlap integral
* @return the value of \f$\int_0^{1/2}dx this (x) other(x)\f$.
* @warning Note that this is not the ordinary overlap which would involve a factor \f$x\f$,
*  the Jacobian in the polar coordinate system.
* @author Richard J. Mathar
* @since 2006-11-29
*/
double Kq::Overl(const Kq & oth) const
{
#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
	return zbase.Overl(q,kvect,oth.q,oth.kvect) ;
#else
	double resul=0. ;
	for(int i=0; i < kvect.size(); i++)
		resul += at(i)*oth.at(i) ;

	/** The approximation is the simple Riemann sum with interval lengths \f$\Delta x=0.5/N()\f$.
	*/
	return resul/(2.*N()) ;
#endif
}

/** The overlap integral of the derivative with another radial function.
* @param[in] oth the other radial function in the overlap integral
* @return the value of \f$\int_0^{1/2}dx x [\partial_x this (x)] other(x)\f$.
* @author Richard J. Mathar
* @since 2006-11-29
* @warning unlike Overl(), this function here is not commutative in the two functions.
*/
double Kq::gradOverl(const Kq & oth) const
{
#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
	return zbase.gradOverl(q,kvect,oth.q,oth.kvect) ;
#else /* KARHLOE_RADIAL_POLYN */
	double resul=0. ;
	for(int i=0; i < kvect.size(); i++)
	{
		/* \f$x\f$ is actually (i+1/2)/(2N). The missing factor is recovered at the return statement
		* One could of course become more efficient by eliminating a factor \c h in gradat() which
		* cancels against one factor here.
		*/
		resul += (i+0.5)*gradat(i)*oth.at(i) ;

		/* debugging
		* cout << __LINE__ << " " << i << " " << gradat(i) << " " << oth.at(i) << endl ;
		*/
	}

	/** The approximation is the simple Riemann sum with interval lengths \f$\Delta x=0.5/N()\f$.
	*/
	return resul/(4*N()*N()) ;
#endif /* KARHLOE_RADIAL_POLYN */
}

/** Size of the representation vector.
* @return number of elements (abscissa points)
*/
int Kq::N() const
{
	return kvect.size() ;
}

/** Neumann factor associated with the azimuthal quantum number
* @return 1 if \f$q=0\f$, 2 else.
*/
int Kq::Neumann() const
{
	return (q ? 2 : 1) ;
}

/** The azimuthal overlap integral.
* @param[in] m1 azimuthal quantum number of first factor
* @param[in] m2 azimuthal quantum number of second factor
* @param[in] isCos flags this as \f$\int d\theta \cos m_1\theta\cos m_2\theta\f$ if \c true,
*    else as \f$\int_0^{2\pi} d\theta \sin m_1\theta\sin m_2\theta\f$ if \c false
* @return the angular integral,
*/
int Kq::azOverl(const int m1, const int m2, const bool isCos)
{
	/* The overlap between two azimuthal functions of the cosine type is 1
	* if the absolute quantum numbers are the same; for the sine case,
	* a different sign in the arguments may toggle the sine in the result.
	*/
	if(isCos)
		if ( abs(m1)==abs(m2) )
			return 1 ;
		else
			return 0 ;
	else
		if ( abs(m1)==abs(m2) )
			return (m1*m2 >0)? 1 : -1 ;
		else
			return 0 ;
}

#if 0
/** The azimuthal overlap integral.
* @param[in] oth
* @param[in] isCos flags this as \f$\int d\theta \cos m_1\theta\cos m_2\theta\f$ if \c true,
*    else as \f$\int_0^{2\pi} d\theta \sin m_1\theta\sin m_2\theta\f$ if \c false
* @return the angular integral
*/
int Kq::azOverl(const Kq & oth, const bool isCos) const
{
	return azOverl(q,oth.q,isCos) ;
}
#endif

#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
ZernikeBase Kq::zbase ;
#endif /* KARHLOE_RADIAL_POLYN, KARHLOE_RADIAL_PFIT */

/** Comparison operator.
* The radial functions are sorted with largest (dominant) eigenvalue first.
* @param[in] left the Kq left from the less sign
* @param[in] right the Kq right from the less sign
* @return true of \c left is smaller (has larger eigenvalue) than \c right
* @author Richard J. Mathar
* @since 2005-12-06
*/
bool operator< (const Kq & left, const Kq & right)
{
	return ( left.lambdasq > right.lambdasq) ;
}

/** Overloaded output operator.
* A simple ASCII representation of a single radial function
* @param[in,out] os the output stream
* @param[in] some the Kq instance to be printed
* @author Richard J. Mathar
* @since 2005-12-06
*/
ostream & operator<<(ostream &os, const Kq & some)
{
		// the values of Wang-Markey's Table III are obtained by dividing lambda^2 through pi
		// the values of Dai's Table 1 "Variance" are obtained by multiplying this in addition with 4
	os << "# lambda^2 " << some.lambdasq << " lambda^2/pi " << some.lambdasq/M_PI << " 4xlambda^2 " 
		<< 4.*some.lambdasq << endl ;
	os << "# p " << some.p << " q +-" << some.q << endl ;
	const int N= some.kvect.size() ;
	for(int j=0;j< N; j++)
	{
#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
		/* debugging
		* os << "# Zern mode " << j << " coef " << some.kvect[j] << endl ;
		*/
		const double r = (j+0.5)/(2.*N) ;
		cout << r << " " << some.K(true,r,0.) << endl ;
#else
		os << (j+0.5)/(2.*N) << " " << some.kvect[j] << endl ;
#endif
	}
	os << endl << endl ;
	return os ;
}

/** Ctor
* @param[in] see a seed value for the Gaussian random numbers of the eigenfunctions
* @param[in] qmax the maximum number of phase screen basis functions.
*      This puts a constraint on the degrees of freedom admitted in the
*      speckles.
* @param[in] N the number of points on the radial mesh used to generate
*       the radial eigenfunctions.
* @param[in] ord a cut-off parameter for the number of terms in the binomial
*       expansions for some kernel integrals,
* \latexonly
* \cite[(27)]{FriedJOSA68} and
* \cite[(39)--(42)]{WangJOSA68}. 
* \endlatexonly
*/
KqBasis::KqBasis(const int see, const int qmax, const int N, const int ord) : Ksort(), scrs()
{
	initseed(see) ;

	KarhLoeve KL(ord) ;
	double * lambda = new double[N] ;
	double * K = new double[N*N] ;
	for(int q=0; q < qmax ; q++)
	{
#if defined KARHLOE_RADIAL_POLYN
		KL.Kpq(q,K,lambda,N,zbase) ;
#else
		KL.Kpq(q,K,lambda,N) ;
#endif
		// scan all the newly generated radial functions; insert them
		// into the vector Ksort and optionally remove at the end with the too small
		// eigenvalues one, if the maximum number of requested modes is reached, so only
		// those with the largest eigenvalues survive.
		for(int e=0; e < N ; e++)
		{
			const int efrom0 = N-1-e ;
			Kq tmp(q,lambda[efrom0],&K[N*efrom0],N,e+1) ;
			vector<Kq>::iterator i= Ksort.begin();
			for( ; i != Ksort.end() ; i++)
				if( tmp < *i)
				{
					Ksort.insert(i,tmp)  ;
					break ;
				}
			if( i == Ksort.end() )
				Ksort.push_back(tmp) ;
			Ksort.resize(qmax) ;
		}

#ifdef KARHLOE_SHOW_RESIZE
		printf("# "__FILE__" %d q=%d lambda^2= %.6f %.6f %.6f %.6f\n",__LINE__,q,lambda[N-1],lambda[N-2],lambda[N-3],lambda[N-4]) ;
		double *DDbar = Dresize(K,N,2,lambda) ;
		delete [] DDbar ;
#endif
	}
	delete [] K ;
	delete [] lambda ;
}

/**
* @param[in] K the radial expansion of the KL eigenmodes of a constant q
* @param[in] shrink a small pupil resizing factor, a divisors of N
* @since 2007-03-29
*/
double * KqBasis::Dresize(const double *K, const int N, const int shrink, const double *lambda) const
{

#if KARHLOE_SHOW_RESIZE == 1
	integer lda = N ;
	integer nrhs = N/shrink ;
#else
	integer lda = N/shrink ;
	integer nrhs = N ;
#endif

	double *Dcoup = new double[lda*nrhs] ;

	/* need a temproray storage because dgesv_() overwrites its A matrix
	*/
	double *A = new double[lda*lda] ;

	/* rather simple solution that doesnot reverse the order (still smallest eigenvalues first
	* memcpy(A,K,N*N*sizeof(double)) ;
	*/

#if KARHLOE_SHOW_RESIZE == 1
	/* the pixels vary along the columns (changing row index), and this order is maintained,
	* the reversal of the eigenvectors (to largest eigenvalues first) means reversal of the column indices
	*/
	for(int row = 0 ; row < lda ; row++)
		for(int col =0 ; col < lda ; col++)
#if KARHLOE_SHOW_TYPE == 1
			// A[row+col*lda] = K[row+(lda-1-col)*lda] ;
			A[KARHLOE_F77_FLAT(row,col,lda)] = K[row+(N-1-col)*N] ;
#else
			/* perhaps solve for the gamma matrix defined to couple the radial values multiplied
			* by the lambda (eigenvalues square roots)
			*/
			if ( lambda[lda-1-col] > 0. )
				A[KARHLOE_F77_FLAT(row,col,lda)] = sqrt(lambda[lda-1-col])*K[row+(lda-1-col)*lda] ;
			else
				A[KARHLOE_F77_FLAT(row,col,lda)] = 0. ;
#endif

	for(int row = 0 ; row < lda ; row++)
		for(int col =0 ; col < nrhs ; col++)
			// Dcoup[row+col*lda] = A[shrink*row+col*lda] ;
			Dcoup[KARHLOE_F77_FLAT(row,col,lda)] = A[KARHLOE_F77_FLAT(shrink*row,col,lda)] ;
#else /* KARHLOE_SHOW_RESIZE */
	/* the pixels vary along the columns (changing row index), and this order is maintained,
	* the reversal of the eigenvectors (to largest eigenvalues first) means reversal of the column indices
	*/
	for(int row = 0 ; row < lda ; row++)
		for(int col =0 ; col < lda ; col++)

#if KARHLOE_SHOW_TYPE == 1
			A[KARHLOE_F77_FLAT(row,col,lda)] = K[shrink*row+(N-1-col)*N] ;
#else
			/* perhaps solve for the gamma matrix defined to couple the radial values multiplied
			* by the lambda (eigenvalues square roots)
			*/
			if ( lambda[N-1-col] > 0. )
				A[KARHLOE_F77_FLAT(row,col,lda)] = sqrt(lambda[N-1-col])*K[shrink*row+(N-1-col)*N] ;
			else
				A[KARHLOE_F77_FLAT(row,col,lda)] = 0. ;
#endif

	for(int row = 0 ; row < lda ; row++)
		for(int col =0 ; col < nrhs ; col++)
#if KARHLOE_SHOW_TYPE == 1
			Dcoup[KARHLOE_F77_FLAT(row,col,lda)] = K[row+(N-1-col)*N] ;
#else
			if ( lambda[N-1-col] >0.)
				Dcoup[KARHLOE_F77_FLAT(row,col,lda)] = sqrt(lambda[N-1-col])*K[row+(N-1-col)*N] ;
			else
				Dcoup[KARHLOE_F77_FLAT(row,col,lda)] = 0. ;
#endif

#endif /* KARHLOE_SHOW_RESIZE */

#if 0
	/* debugging input
	*/
	cout << "# solving A" << endl ;
	for(int row = 0 ; row < lda && row < 10; row++)
	{
		for(int col =0 ; col < lda && col < 10; col++)
			printf("%.4f ",A[KARHLOE_F77_FLAT(row,col,lda)]) ;
		printf("\n") ;
	}

	cout << "# solving B" << endl ;
	for(int row = 0 ; row < nrhs && row < 10 ; row++)
	{
		for(int col =0 ; col < lda && col < 10 ; col++)
			printf("%.4f ",Dcoup[KARHLOE_F77_FLAT(row,col,lda)]) ;
		printf("\n") ;
	}
#endif

	integer info ;
	integer *ipiv = new integer[N] ;

	/* solve the system of linear equations K(x=r/shrink)= K(r)*Dcoup by a call to LAPACK's dgesv().
	* A of the dgesv is represented by K(r), and B by K(x=r/shrink), and X by Dcoup (on return stored in B).
	*/

	dgesv_(&lda,&nrhs,A,&lda,ipiv,Dcoup,&lda,&info) ;
	if ( info)
		cerr << __LINE__ << " dgesv info " << info << endl ;

	delete ipiv ;
	delete A ;

#if 0
	/* debugging  ASCII output
	*/
	for(int row = 0 ; row < lda ; row++)
	{
		for(int col =0 ; col < nrhs ; col++)
		{
			printf("%.4f ",Dcoup[KARHLOE_F77_FLAT(row,col,lda)]) ;
		}
		if ( row == nrhs)
			printf("\n------") ;
		printf("\n") ;
	}
#endif

#if 1
	/* debugging  gnuplot output. lambda still contains the square here, which we convert
	* to the values themselves for display
	*/
	printf("\n\n") ;
	for(int col =0 ; col < nrhs ; col++)
	{
		for(int row = 0 ; row < lda ; row++)
		{
			/* lambda still contains the square here, which we would convert to lambda for display
			*/
			//printf("%.1f %d %e %e\n",row-0.5,col,fabs(sqrt(lambda[lda-1-row])*Dcoup[row+col*lda]),fabs(Dcoup[row+col*lda])) ;
			printf("%.1f %d %e\n",row-0.5,col,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
			printf("%.1f %d %e\n",row+0.5,col,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
		}
		printf("\n") ;
		for(int row = 0 ; row < lda ; row++)
		{
			printf("%.1f %.1f %e\n",row-0.5,col+0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
			printf("%.1f %.1f %e\n",row+0.5,col+0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
		}
		printf("\n") ;
	}
#endif

#if 0
	/* debugging  Maple output
	*/
	printf("dta :=[\n\n") ;
	for(int col =0 ; col < nrhs ; col++)
	{
		if( col != 0 )
			printf(",\n") ;
		printf("[") ;
		for(int row = 0 ; row < lda ; row++)
		{
			if( row != 0 )
				printf(",\n") ;
			printf("[%d,%f,%e],",col,row-0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
			printf("[%d,%f,%e]",col,row+0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
		}
		printf("],") ;
#if 0
		if( col != 0 )
			printf(",\n") ;
#endif
		printf("[") ;
		for(int row = 0 ; row < lda ; row++)
		{
			if( row != 0 )
				printf(",\n") ;
			printf("[%f,%f,%e],",ccol+0.5,crow-0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
			printf("[%f,%f,%e]",ccol+0.5,crow+0.5,fabs(Dcoup[KARHLOE_F77_FLAT(row,col,lda)])) ;
		}
		printf("]") ;
	}
	printf("] ;\n") ;
	printf("interface(plotdevice=\"ps\",plotoutput=\"tmp.ps\",plotoptions=\"noborder,portrait\");\n") ;
	printf("plots[surfdata](dta,axes=FRAMED,labels=[p,p,\"\"],orientation=[-148,54]);\n") ;
#endif

	/* save the eigenvectors into return area
	*/
	return Dcoup ;
}

#if 0
/** Ctor
* Not needed and incompatible with the other ctor because that one would
* generate the resizing matrix.
* @param[in] src
* @param[in] see
*/
KqBasis::KqBasis( vector<Kq> & src, const int see) : Ksort(src), scrs()
{
	initseed(see) ;
}
#endif

#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
ZernikeBase KqBasis::zbase ;
#endif /* KARHLOE_RADIAL_POLYN */

/** Dtor
*/
KqBasis::~KqBasis()
{
	/* debugging
	* cout << "# dtor KqBasis " << __FILE__ << " " << __LINE__ << endl ;
	*/

	/** The memory allocated in makeScreen() to keep the vector of the Taylor screens
	* is de-allocated.
	*/
	for(int s=0; s < scrs.size() ; s++)
		delete [] scrs[s] ;
}

/** Produce a gnuplot type of output for the FT of the radial functions
* @since 2007-04-29
* @author Richard J. Mathar
*/
void KqBasis::showFT() const
{
	const int Nradial = Ksort[0].N() ;
	for(int no=0; no < Ksort.size() ; no++)
	{
		cout << "# q= " << Ksort[no].q << " lambda^2=" << Ksort[no].lambdasq << " index " << no << endl ;
		/* the x have been tabulated along integers i=0,1,2... so the natural resolution
		* on the sigma axis is j=0,1,2,...  for 2*pi*sigma*x = 2*pi*j*i/Nradial 
		*/
		for (int s=0 ; s < 2*Nradial ; s++)
		{
			const double sigm = 0.05*s ;
			cout << sigm << " " << Ksort[no].fft(sigm) << endl ;
		}
		cout << endl << endl ;
	}
}

/** produce another single screen and add it to the scrs[]
* @param[in] Dr0 the ratio \f$D/r_0\f$ of pupil diameter and Fried parameter
* @param[in] withtipt a flag indicating whether the tip-tilt component ought remain
*	in the screen generated (\c true) or be removed (\c false)
*/
void KqBasis::makeScreen(double Dr0, const bool withtipt)
{
	if ( withtipt)
		/* only piston removed */
		makeScreen(Dr0,0) ;
	else
		/* piston and tip-tilt removed */
		makeScreen(Dr0,2) ;
}

/** produce another single screen and add it to the scrs[]
* @param[in] Dr0 the ratio \f$D/r_0\f$ of pupil diameter and Fried parameter
* @param[in] aoOrd the number of low-order modes skipped. 3 means only tip-tilt removed.
*/
void KqBasis::makeScreen(double Dr0, const int aoOrd)
{
	/** A cartesian coordinate pair in the pupil square (which has the side
	* length \f$D\f$), normalized to \f$0\le x[]\le 1/2\f$.
	*/
	double x[2] ={0.,0.} ;

	/** The number of radial absissa points in the numerical representation
	* of the radial eigenfunctions, or the number of different Zernike polynomials
	* if this has been chosen as the representation..
	*/
	const int Nradial = Ksort[0].N() ;

	/** The distance between two abscissa points on the radial mesh. This
	* is equivalent to the value used in the Riemann sense of integration over
	* the interval [0,1/2].
	*/
	const double xstep = 0.5/Nradial ;

	/** Further down we only need the fluctuation amplitude
	* \f$ (D/r_0)^{5/6}\f$. The factor \f$D^2\f$ of Fried's equation (23) is already
	* implicitly contained in the fact that our basis functions are normalized to be
	* 1 if integrated from x=0 to x=1/2, which needs to become an integral over the
	* r=0 to r=D/2 and implicitly divided thru \f$D^2\f$ to match the equation.
	*/
	Dr0 = pow(Dr0,5./6.) ;

	/** The flac \c idist indicates to the \c dlarnv function of LAPACK that we need Gaussian distributed
	* random variables.
	*/
	integer idist = 3;

	/** The number of independent eigenfunctions is kept in \c iN
	*/
	integer iN = modes() ;
	/** Each mode is modelled for each instance of the phase screen with an
	* individual randomized amplitude with Gaussian distribution. All these
	* are kept for one single screen in the vector \c rstd . 
	*/
	double * rstd = new double[iN] ;

	/** A call to dlarnv() of the LAPACK library produces \c iN Gaussian random numbers
	* with mean 0 and standard deviation 1.
	*/
	dlarnv_(&idist,lapsee,&iN,rstd) ;

	/** A square grid of width and height equal to the diameter needs two times
	* as much grid points along \f$x\f$ and two times as much grid points along \f$y\f$
	* as there were abscissa values in the radial meshes, if one wants to match
	* the  resolution in the KL eigenfunctions sensibly. For each point we allocate
	* one value (that is the real-valued phase at that point in the phase screen).
	*/
	double *scr = new double[4*Nradial*Nradial] ;

#if 0
	/* debugging. Explicitly document the modes and their eigenvalues that will
	* be skipped with this level of AO correction. The actual skipping is done
	* later below in the loop.
	*/
	for(int no=0; no < aoOrd ; no++)
	{
		double lambdasq ;
		K(no,x,lambdasq) ;
		cout << "# skipping " << no << " lambda^2 " << lambdasq << endl ;
	}
#endif

	/** The square grid is filled by a double loop over the two orthogonal axes.
	*/
	for(int xindx=0 ; xindx < 2*Nradial ; xindx++)
	{
		x[0] = -0.5+(0.5+xindx)*xstep ;
		for(int yindx=0 ; yindx < 2*Nradial ; yindx++)
		{
			x[1] = -0.5+(0.5+yindx)*xstep ;
			double resul=0. ;

			/** The value on the grid is a sum over all eigenmodes (which are basis functions)
			* scaled to the physically relevant units.
			*/
			for(int no=0; no< iN;no++)
			{
				double lambdasq ;

				/** For each point in the grid we call K() to recover the
				* normalized value of the eigenmode and to query its eigenvalue.
				* On return, K() has not yet multiplied these with the square root
				* of \f$\lambda^2\f$.
				*/
				double val = K(no,x,lambdasq) ;	// scaling with lambda not in here

				/** If this particular mode is the tip-tilt mode, we skip it. This
				* can be sensed by comparison of the squared eigenvalue with the threshold
				* value of 0.0178 for the next higher mode, in comparison to the value
				* \f$ 0.3521=0.1121\pi\f$ of the tip-tilt mode(s).
				* if ( !withtipt && lambdasq > 0.03)	// tip-tilt has lambdasquared = 0.3521=pi*0.1121
				*continue ;
				*/
				if ( no < aoOrd)
					continue ;
				if( lambdasq > 0.)
				{
#ifndef KARHLOE_RETURNWI_LAMBDA
					/* with the normalization of the basis function over x=0..1/2 and
					* theta =0 .. 2*pi equal to 1, we need a variance of the new random
					* number of (D/r0)^5/3*lambdasq. We upscale the values of rstd[] from
					* a standard deviation of 1 to this value by multiplying the entire set of
					* random numbers with the square root of the std dev.
					*/
					val *= sqrt(lambdasq) ;
#endif
					val *= rstd[no] ;
					val *= Dr0;	// already powered above with power 5/6
					resul += val ;
				}
			}
			/** The same over all these properly scaled eigenmodes is kept
			* in the linear vector \c scr . There is no sense of using the language
			* of ``row-order'' or ``column-order'' here; the stride is set by the side length
			* of the phase screen square.
			*/
			scr[xindx+2*Nradial*yindx] = resul ;
		}
	}

	/** As a clean-up procedure, the set of random numbers allocated above is de-allocated.
	*/
	delete [] rstd ;

	/** A pointer to this particular phase screen (but not the data themselves)
	* is stored in \c scr for future use.
	*/
	scrs.push_back(scr) ;	// store the pointer to the data (to be released by the dtor)

	scrsIsBasis=false ;
}

/** produce the basis function no and add it to the scrs[]
* @param[in] Dr0 the ratio \f$D/r_0\f$ of pupil diameter and Fried parameter
* @param[in] no the number of the basis function, starting at 0
*/
void KqBasis::makeBasis(double Dr0, const int no)
{
	const int iN = modes() ;
	if ( no >= iN)
	{
		cerr << __FILE__ << " " << __LINE__ << " no=" << no << ">=" << iN << endl ;
		return ;
	}
	/** A cartesian coordinate pair in the pupil square (which has the side
	* length \f$D\f$), normalized to \f$0\le x[]\le 1/2\f$.
	*/
	double x[2] ={0.,0.} ;

	/** The number of radial absissa points in the numerical representation
	* of the radial eigenfunctions.
	*/
	const int Nradial = Ksort[0].N() ;

	/** The distance between two abscissa points on the radial mesh. This
	* is equivalent to the value used in the Riemann sense of integration over
	* the interval [0,1/2].
	*/
	const double xstep = 0.5/Nradial ;

	/** Further down we only need the fluctuation amplitude
	* \f$ (D/r_0)^{5/6}\f$. The factor \f$D^2\f$ of Fried's equation (23) is already
	* implicitly contained in the fact that our basis functions are normalized to be
	* 1 if integrated from x=0 to x=1/2, which needs to become an integral over the
	* r=0 to r=D/2 and implicitly divided thru \f$D^2\f$ to match the equation.
	*/
	Dr0 = pow(Dr0,5./6.) ;

	/** A square grid of width and height equal to the diameter needs two times
	* as much grid points along \f$x\f$ and two times as much grid points along \f$y\f$
	* as there were abscissa values in the radial meshes, if one wants to match
	* the  resolution in the KL eigenfunctions sensibly. For each point we allocate
	* one value (that is the real-valued phase at that point in the phase screen).
	*/
	double *scr = new double[4*Nradial*Nradial] ;

	/** The square grid is filled by a double loop over the two orthogonal axes.
	*/
	for(int xindx=0 ; xindx < 2*Nradial ; xindx++)
	{
		x[0] = -0.5+(0.5+xindx)*xstep ;
		for(int yindx=0 ; yindx < 2*Nradial ; yindx++)
		{
			x[1] = -0.5+(0.5+yindx)*xstep ;
			double resul=0. ;

			/** The value on the grid is the basis function
			* scaled to the physically relevant units.
			*/

			double lambdasq ;

			/** For each point in the grid we call K() to recover the
			* normalized value of the eigenmode and to query its eigenvalue.
			* On return, K() has not yet multiplied these with the square root
			* of \f$\lambda^2\f$.
			*/
			double val = K(no,x,lambdasq) ;	// scaling with lambda not in here

			if( lambdasq > 0.)
			{
#ifndef KARHLOE_RETURNWI_LAMBDA
				/* with the normalization of the basis function over x=0..1/2 and
				* theta =0 .. 2*pi equal to 1, we need a variance of the new random
				* number of (D/r0)^5/3*lambdasq. We upscale the values of rstd[] from
				* a standard deviation of 1 to this value by multiplying the entire set of
				* random numbers with the square root of the std dev.
				*/
				val *= sqrt(lambdasq) ;
#endif
				val *= Dr0;	// already powered above with power 5/6
				resul = val ;
			}

			/** The same over all these properly scaled eigenmodes is kept
			* in the linear vector \c scr . There is no sense of using the language
			* of ``row-order'' or ``column-order'' here; the stride is set by the side length
			* of the phase screen square.
			*/
			scr[xindx+2*Nradial*yindx] = resul ;
		}
	}

	/** A pointer to this particular phase screen (but not the data themselves)
	* is stored in \c scr for future use.
	*/
	scrs.push_back(scr) ;	// store the pointer to the data (to be released by the dtor)
	scrsIsBasis=true ;
}

/** produce all basis functions and add it to the scrs[]
* @param[in] Dr0 the ratio \f$D/r_0\f$ of pupil diameter and Fried parameter
*/
void KqBasis::makeBasis(double Dr0)
{
	for(int no=0; no < modes() ; no++)
		makeBasis(Dr0,no) ;
}

/** Read one value of a previously generated phase screen pixel.
* @param[in] screenNo the index of the phase screen, 0 -based.
* @param[in] xindx the integer coordinate of the \f$x\f$ value
* @param[in] yindx the integer coordinate of the \f$y\f$ value
* @param[in] scrV a pointer to a vector of screens. If the parameter is NULL,
*       the collection of screens in KqBasis::scrs is used.
* @return the value that had been generated with the \c screenNo 'th call to screen() before.
* @todo include a check for out-of-bound indices
*/
double KqBasis::inScreen(const int screenNo, const int xindx, const int yindx,const vector<double*> *scrV) const
{
	/** This reads immediately, with no protection against index over/underflow,
	* one value of the screen set.
	*/
	if ( scrV)
		return (*scrV)[screenNo][xindx+2*Ksort[0].N()*yindx] ;
	else
		return scrs[screenNo][xindx+2*Ksort[0].N()*yindx] ;
}

/** Read one value of a previously generated phase screen pixel.
* @param[in] screenNo the index of the phase screen, 0 -based.
* @param[in] r the radial distance to the pupil center between 0 and 0.5
* @param[in] theta the azimuthal coordinate between 0 to \f$2\pi\f$
* @return the value that had been generated with the \c screenNo 'th call to screen() before.
*       Points outside the circular aperture are returned as zero.
*/
double KqBasis::inScreen(const int screenNo, const double r, const double theta) const
{
	/* in a Cartesian system, the values -1/2<=x<=1/2 and -1/2<=y<=1/2 are mapped
	* onto 0<=xindx < 2*N() and 0<=yindx < 2*N().
	* We invert this transformation to
	* get xindx and yindx: (2*x+1)/N() = 0.5+xindx; (2*y+1)/N()=0.5+yindx
	*/
	double x[2] ;
	x[0] = r*cos(theta) ;
	x[1] = r*sin(theta) ;

	/** The number of radial absissa points in the numerical representation
	* of the radial eigenfunctions.
	*/
	const int Nradial = Ksort[0].N() ;

#if 0
	/* this is a quasi reversed but unnecessarily complicated implementation of
	* the formula used before.
	*/
	const int xindx = round((2.*x[0]+1.)*Nradial-0.5) ;
	const int yindx = round((2.*x[1]+1.)*Nradial-0.5) ;
#else
	const int xindx = (2.*x[0]+1.)*Nradial ;
	const int yindx = (2.*x[1]+1.)*Nradial ;
#endif
	if ( xindx >=0 && xindx < 2*Nradial && yindx >=0 && yindx < 2*Nradial)
		return inScreen(screenNo,xindx,yindx) ;
	else
		return 0. ;
}

/** Reproduce the phase screens previously generated with makeScreen() in various output formats.
* @param[in] Dr0 the value of \f$D/r_0\f$. This value is not used again to scale the actual
*        phases, but only used to generate a FITS header line.
* @param[in] D the value of \f$D\f$. This value is not used again to scale the actual
*        phases, but only used to generate a FITS header line.
* @param[in] slrnseed a seed value for the random number generator. This value is not actually used
* 	to compute again any phase value, but only used to generate a FITS header line.
* @param[in] withtipt This parameter recalls whether tip-tilt is still in the phase screen (\c true)
*	or removed (\c false). Only used to generate a parameter in the FITS header.
* @param[in] mapletype flag to request that the ASCII output is suitable for display with Maple .
* @param[in] clip if \c true demands that the phase screens are confined to the circular
*	central piece of the pupil. This is the region to which the phase statistics has been
*	adapted.
* @param[in] fits if \c cfitsio support had been compiled in, and if this is not a \c NULL pointer,
*	the output is a FITS file, not an ASCII report to stdout .
* @param[in] scrV a pointer to a vector of screens. If the parameter is NULL,
*       the collection of screens in KqBasis::scrs is used.
* Since \c cfitsio does not declare the variables \c const where possible, we avoid using this here.
*/
void KqBasis::out(double D, double Dr0, int slrnseed, bool withtipt, bool mapletype, bool clip, char *fits,
	const vector<double*> * scrV) const
{
	double x[2] ;
	const int Nradial = Ksort[0].N() ;
	const double xstep = 0.5/Nradial ;
	const int scrCount = (scrV == 0)? scrs.size() : scrV->size() ;


	/* debugging
	* cerr << __FILE__ << " " << __LINE__ << " " << fits << " " << Nradial << " " << scrCount << endl ;
	*/

	if( fits == 0 )
	{
		for(int screenNo=0; screenNo < scrCount ;screenNo++)
		{
			if ( scrV == NULL)
			{
				bool usecos ;
				int j=mod2rad(screenNo,usecos) ;
				if ( usecos)
					cout << "# q=+" ;
				else
					cout << "# q=-" ;
				cout << Ksort[j].q << " p=" << Ksort[j].p << " lambdasq=" << Ksort[j].lambdasq << endl ;
			}
			if( mapletype)
				cout << "dat" << screenNo << " := [" ;
			bool firstx = true ;
			for(int xindx=0 ; xindx < 2*Nradial ; xindx++)
			{
				x[0] = -0.5+(0.5+xindx)*xstep ;
				if( mapletype)
					if(firstx)
						cout << "[" ;
					else
						cout << ",[" ;
				bool firsty = true ;
				for(int yindx=0 ; yindx < 2*Nradial ; yindx++)
				{
					x[1] = -0.5+(0.5+yindx)*xstep ;
					if ( clip && hypot(x[0],x[1]) > 0.5)
						continue ;
					const double resul= inScreen(screenNo,xindx,yindx,scrV) ;
					if( mapletype)
					{
						if ( !firsty)
							cout << "," ;
						cout << "[" << x[0] << "," << x[1] << "," << resul << "]\n" ;
					}
					else
						cout << x[0] << " " << x[1] << " " << resul << endl ;
					firsty= false ;
				}
				if( mapletype)
					cout << "]" << endl ;
				cout << endl ;
				firstx = false ;
			}
			if( mapletype)
				cout << "]:" << endl ;
			cout << endl ;
		}
		if ( mapletype)
		{
			cout << "interface(plotdevice=postscript) :\n" ;
			cout << "interface(plotoptions=`portrait,noborder`) :" << endl  ;
			for(int screenNo=0; screenNo < scrCount ;screenNo++)
			{
				cout << "interface(plotoutput=`tmp" << screenNo << ".ps`) :" << endl  ;
				bool usecos ;
				int j=mod2rad(screenNo,usecos) ;
				cout << "plots[surfdata](dat" << screenNo ;
				if ( scrV)
				{
					cout << ",title=\"" << screenNo << "\"" ;
				}
				else if ( scrsIsBasis)
				{
					cout << ",title=\"q="  ;
					if (usecos) 
						cout << Ksort[j].q ;
					else
						cout << -Ksort[j].q ;
					cout << " p=" << Ksort[j].p << "\"" ;
				}
				else
				{
					// TODO: generate time labels
					cout << ",title=\"" << screenNo << "\"" ;
				}
				//cout << ",axes=NORMAL,labels=[x,y,\"\"],tickmarks=[1,1,1],orientation=[60,45],scaling=CONSTRAINED) ;\n" ;
				cout << ",axes=NORMAL,labels=[x,y,\"\"],tickmarks=[1,1,1],orientation=[60,45]) ;\n" ;
			}
		}
	}
#ifdef HAVE_CFITSIO
	else
	{
		fitsfile * fptr;
		int status=0 ;
		/** Use fits_create_file() of \c cfitsio to create the file. This implies that we'll return
		* with an error condition if a file of the same name already exists.
		*/
		fits_create_file(&fptr,fits,&status) ;
		if( status)
			fits_report_error(stderr,status) ;
		else
		{

			/** This becomes a simple image cube, one phase screen per images.
			*/
			long naxes[3] ;
			naxes[0]=naxes[1]=2*Nradial ; naxes[2] = scrCount ;

			long fpixel[3] ;
			fpixel[0]=fpixel[1]=fpixel[2] = 1 ;

			char * srcfile = __FILE__ ;
			fits_create_img(fptr,FLOAT_IMG,3,naxes,&status) ;

			/** All the screens generated with makeScreen() are copied to the FITS file in the
			* order of creation.
			*/
			for(int screenNo=0; screenNo < scrCount ;screenNo++)
			{
				fpixel[2] = screenNo+1 ;
				if ( scrV )
					fits_write_pix(fptr,TDOUBLE,fpixel,4*Nradial*Nradial,(*scrV)[screenNo],&status) ;
				else
					fits_write_pix(fptr,TDOUBLE,fpixel,4*Nradial*Nradial,scrs[screenNo],&status) ;
			}
			fits_movabs_hdu(fptr,1,0,&status) ;

			/** We install primary header keywords \c CREATOR, \c DR0, \c SEED and \C TIPTILT
			* to indicate the parameters of creation.
			*/
			fits_write_key(fptr,TSTRING,"CREATOR",srcfile,"Richard J. Mathar",&status) ;
			fits_write_key(fptr,TDOUBLE,"D",&D,"[m] pupil diameter",&status) ;
			fits_write_key(fptr,TDOUBLE,"DR0",&Dr0,"D over r0",&status) ;
			fits_write_key(fptr,TINT,"SEED",&slrnseed,"seed of random number generator",&status) ;
			int ttilt = withtipt ;
			fits_write_key(fptr,TLOGICAL,"TIPTILT",&ttilt,"false if tip-tilt removed",&status) ;
			fits_write_date(fptr,&status) ;
			fits_close_file(fptr,&status) ;
			if( status)
				fits_report_error(stderr,status) ;
		}
	}
#endif
}

/** Return the number of available modes.
* @return The number of modes, counting those with nonzero azimuthal quantum number
*   with multiplicity (twice).
*/
int KqBasis::modes() const
{
	int iN=0 ;
	for(int i=0 ; i < Ksort.size() ; i++)
		iN += Ksort[i].q ? 2 : 1 ;
	return iN ;
}

/** Return the radial mode for each of the available pupil modes.
* @param[out] iscos \c true if this is a cosine in azimuth, \c false if this is a sine
* @return 0,1,1,2,2,3,3 etc if \c m = 0,1,2,3,...
* @author Richard J. Mathar
* @since 2006-11-29
*/
int KqBasis::mod2rad(int m, bool & iscos) const
{
	/* The basis assumption is that we did not remove any modes when the basis
	* functions were created, only the q=0 piston mode.
	*/
	int j=0 ;
	iscos= true ;
	while(m--)
	{
		if ( j >= Ksort.size() )
		{
			cerr << __FILE__ << " " << __LINE__ << " intrn error "  << Ksort.size() << endl ;
			exit(EXIT_FAILURE) ;
		}
		if ( Ksort[j].q) 	// have cos and sine components
		{
			if ( !iscos)
				// previous was a sine (2nd call) so move on to the next radial function
				j++ ;
			iscos= !iscos ;	// toggle the value
		}
		else
		{
			j++ ;
			iscos=true ;	// always report iscos=true for the q=0 components
		}
	}
	return j ;
}

/**
* @param[in] i the index of the basis function to be used, starting at 0 for the
*      one with the largest eigenvalue (tip-tilt).
* @param[in] x the 2 Cartesian coordinates of the point on the screen. Both in \f$[0,1/2]\f$.
* @param[out] lambdasq the normalized square of the KH eigenvalue
* @return the value of the eigenmode.
*   This value is not yet multiplied by the eigenvalue \f$\lambda\f$, nor by any
*   weight derived from the Gaussian statistics.
*/
double KqBasis::K(const int i, double x[2], double & lambdasq) const
{
	/** The point on the radial mesh is given by the Pythagoran distance
	* from the pupil center.
	*/
	const double r=hypot(x[0],x[1]) ;
	/** Initialize the result to 0. This will be the value returned if the
	* point selected by \c x falls outside the circular center of the pupil.
	*/
	double resul=0. ;
	if( r > 0.5 || i>= modes() )
		return resul ;

	/** expansions \f$\propto \cos q\theta\f$ or \f$\propto \sin q\theta\f$.
	*/
	bool usecos  ;

	/** \f$i=0\f$ selects the first eigenvector, \f$i=1,2\f$ use the \f$\cos \theta\f$
	* and \f$\sin\theta\f$ azimuthal modes, \f$i=3,4\f$ the \f$\cos 2\theta\f$ and
	* \f$\sin 2\theta\f$ azimuthal modes, etc.
	*/
	const int j= mod2rad(i,usecos) ;

	/** If the value of \c i has been larger than the number of eigenmodes, we return 0.
	* Otherwise the value is computed by a call to Kq::K() .
	*/
	if( j < Ksort.size())
	{
		resul = Ksort[j].K(usecos,x) ;
		lambdasq = Ksort[j].lambdasq ;
		/* already included or not included in the call to Kq::K(), depending on
		* the cpp macro KARHLOE_RETURNWI_LAMBDA.
		* resul *= sqrt(Ksort[j].lambdasq) ;
		*/
	}
	return resul ;
}

/** Feed a new seed into the random number generator of the screens.
* @param[in] see the new seed
*/
void KqBasis::initseed( const int see)
{
	lapsee[0]=lapsee[1]=lapsee[2]=0 ;
	lapsee[3] = (see & 4095) | 1 ;
}

/**
* A simple ASCII representation of all radial functions
* @param[in,out] os the output stream
* @param[in] some the KqBasis instance to be documented
* @author Richard J. Mathar
* @since 2005-12-06
*/
ostream & operator<<(ostream &os, const KqBasis & some)
{
#ifndef KARHLOE_SHOW_RESIZE
	for(int q=0 ; q < some.Ksort.size() ; q++)
	{
		os << "# index " << q << endl ;
		os << some.Ksort[q] << endl ;
	}
#endif
	return os ;
}
#undef KARHLOE_F77_FLAT

#endif /* KARHLOE_C */