#ifndef KARHLOE_H
#define KARHLOE_H

/**
 @author Richard J. Mathar
 @since 2005-12-06
 <a href="http://www.strw.leidenuniv.nl/~mathar">Richard J. Mathar homepage</a>.
*/
#include <iostream>
#include <cmath>
#include <vector>
#include <string.h>

extern "C" {
/** The f2c header file f2c.h.
* A prerequisite to include clapack.h below.
* This is obtained from <a href="http://www.netlib.org/f2c">www.netlib.org/f2c</a>.
*/
#include <f2c.h>
}

#ifdef HAVE_CFITSIO
/** cfitsio support, if compiled in.
* If cfitsio is available and support has been selected with the \c HAVE_CFITSIO
* compiler macro, we need \c fitsio.h as from
* <a href="http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html">heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html</a>.
*/
#include <fitsio.h>
#endif

/* Define whether the radial functions are on pixel or polynomial
* basis functions. Uncomment for the former.
* @since 2007-03-22
*/
//#define KARHLOE_RADIAL_POLYN

/* Define whether the pixelized radial basis is converted to a 
* polynomial basis.
*/
//#define KARHLOE_RADIAL_PFIT


#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
#include "Zernike.h"
#endif

using namespace std ;

/* Define whether the calls to the interface of the 
* essentially return a value multiplied by the square root of the KL
* eigenvalue, or not. The latter means the values are such that the
* integrals over the pupil remain normalized to unity.
* @since 2007-03-09
*/
#define KARHLOE_RETURNWI_LAMBDA

/** An auxiliary class that builds various integrals associated with the
* Karhunen-Loeve integral equation for Kolmogorov statics of the phase
* structure function.
* @author Richard J. Mathar
* @since 2005-12-06
*/
class KarhLoeve {
public:

	/** The cut-off order for binomial expansions of integrals.
	*/
	int ord; 

	KarhLoeve(const int maxthetord) ;


#if defined KARHLOE_RADIAL_POLYN
	void Kpq(const int q, double *K, double *lambda, const int N, ZernikeBase & zbase) ;
#else
	void Kpq(const int q, double *K, double *lambda, const int N) ;
#endif

#ifndef KARHLOE_RADIAL_POLYN
	double B(const int n, const double x) ;
#endif

	static double binom(double x, int l) ;

#ifndef KARHLOE_RADIAL_POLYN
	static double cosInt(int l) ;

	static double cosInt(const int l, int q) ;
#endif

#if 0
        static double I2(const int j, const int q, const int o) ;

        static double I1(const int j, const int s, const int q, const int o) ;
#endif

#ifdef KARHLOE_RADIAL_POLYN
        static double Ialpha(int j, int s, const int q) ;

        static double S(const int j, const int s, const int q) ;

        static double S(const int j, const int s, const int q, ZernikeBase & z) ;
#endif

#ifndef KARHLOE_RADIAL_POLYN
	double B2() ;

	double Rsymm(double x, double xprime, const int q) ;
#endif

protected:
private:

#ifndef KARHLOE_RADIAL_POLYN
	static double B2int(const int l) ;

	static double B2intaux(const int twoml, const double x) ;
#endif

} ; /* KarhLoeve */

/**
* A single Karhunen-Loeve radial function. It represents one of the eigenvectors
* of the KH eigenvalue equation for a specific azimuthal quantum number. It does not explicitly state
* whether this is of the cosine or sine type in azimuthal terms.
* @author Richard J. Mathar
* @since 2005-12-06
*/
class Kq {
public:
#ifdef KARHLOE_RADIAL_POLYN
	/** The values on the radial mesh, one value per Zernike basis function
	*/
#else
	/** The values on the radial mesh, one value per mesh point.
	*/
#endif
	vector<double> kvect ;

#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
	/** the Zernike radial polynomials.
	* @todo it ought be superfluous to have such a static base here
	* and also in the KqBasis class.
	*/
	static ZernikeBase zbase ;
#endif

	/** The azimuth parameter. Since the kernel of the integral equation depends
	* only on the absolute value \f$|q|\f$, only positive values \f$0,1,2,\ldots\f$ are stored here.
	*/
	int q;

	/** Random number generator seed.
	*/
	int see;

	/** Index to differentiate eigenvectors for a single value of \f$q\f$.
	*/
	int p ;

	/** Normalized eigenvalue. Normalized means it
	* refers to the mesh on radial values from \f$x=0\f$ to \f$x=1/2\f$ and
	* is not multiplied with any power of \f$D/r_0\f$ or \f$D\f$. It has
	* been written as a variable squared in the literature to indicate
	* that it is a real-valued and positive.
	*/
	double lambdasq ;

	Kq() ;

	Kq(const int qval, const double lambda, double * K , const int N, const int pval= 99, const int rseed = 1) ;

	double K(const bool usecos, const double r, const double theta) const ;

	double K(const bool usecos, double x[2]) const ;

	double fft(const double sigma) const ;

#ifndef KARHLOE_RADIAL_POLYN
	double at(const int idx) const ;

	inline double gradat(const int idx) const ;
#endif


	double Overl(const Kq & oth) const ;

	double gradOverl(const Kq & oth) const ;

	int N() const ;

	int Neumann() const ;

	static int azOverl(const int m1, const int m2, const bool isCos) ;

#if 0
	int azOverl(const Kq & oth, const bool isCos) const ;
#endif

protected:
private:
} ;	/* Kq */

bool operator< (const Kq & left, const Kq & right) ;

ostream & operator<<(ostream &os, const Kq & some) ;

/**
 A set of Karhunen-Loeve radial functions, sorted according to 
 decreasing eigenvalue.
 @author Richard J. Mathar
 @since 2005-12-06
*/
class KqBasis {
public:
	/** a list of radial meshes. Each member of the list is a tabulation
	* of a KH eigenvector along the radial separation to the middle of the screen.
	*/
	vector<Kq> Ksort ;

#if defined KARHLOE_RADIAL_POLYN || defined KARHLOE_RADIAL_PFIT
	/** the Zernike radial polynomials.
	*/
	static ZernikeBase zbase ;
#endif /* KARHLOE_RADIAL_POLYN */

	/** A collection of phase screens on a square grid.
	* Each pointer represents one random draw with a statistics in compliance with
	* the Kolmogorov structure function, or a bare basis function.
	*/
	vector<double *> scrs ;

	/** denote whether scrs collected Kolmogorov screens or KL basis functions
	*/
	bool scrsIsBasis ;

	KqBasis(const int see, const int qmax, const int N, const int ord) ;

#if 0
	KqBasis( vector<Kq> & src, const int see) ;
#endif

	~KqBasis() ;

	void makeScreen(double Dr0, const int aoOrd) ;

	void makeScreen(double Dr0, const bool withtipt) ;

	void makeBasis(double Dr0, const int no) ;
	void makeBasis(double Dr0) ;

	double inScreen(const int screenNo, const int xindx, const int yindx, const vector<double*> *scrV=0) const ;
	double inScreen(const int screenNo, const double r, const double theta) const ;

	void out(double D, double Dr0, int slrnseed, bool withtipt, bool mapletype, bool clip, char *fits, const vector<double*> * scrV = 0 ) const ;

	void showFT() const ;

	int modes() const ;

	int mod2rad(const int m, bool & iscos) const ;

	double * Dresize(const double *K, const int N, const int shrink, const double *lambda) const ;

protected:

	double K(const int i, double x[2], double & lambdasq) const ;
	/** Seeds for LAPACK's dlarn()
	*/
	integer lapsee[4] ;
private:

	void initseed( const int see) ;

} ;

ostream & operator<<(ostream &os, const KqBasis & some) ;

/** undefine some defines of f2c.h because they'll collide with some complex<double>
* definitions of /usr/include/c++ which are used in fibmod.C
*/
#undef abs
#undef dabs
#undef min
#undef max
#undef dmin
#undef dmax

#endif /* KARHLOE_H */