#ifndef INTEGRAL_H
#define INTEGRAL_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * integral.h
 * Integrals of pairs of harmonic functions, zero counting, maxima
 */

/* integral_x0^x1 cos(p1 (x-r1))*cos(p2 (x-r2)) dx */
double intfCC(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 cos(p1 (x-r1))*sin(p2 (x-r2)) dx */
double intfCS(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 sin(p1 (x-r1))*cos(p2 (x-r2)) dx */
double intfSC(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 sin(p1 (x-r1))*sin(p2 (x-r2)) dx */
double intfSS(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 exp(p1 (x-r1))*exp(p2 (x-r2)) dx */
double intfEE(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 exp(p1 (x-r1))*cos(p2 (x-r2)) dx */
double intfEC(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 cos(p1 (x-r1))*exp(p2 (x-r2)) dx */
double intfCE(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 sin(p1 (x-r1))*exp(p2 (x-r2)) dx */
double intfSE(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 exp(p1 (x-r1))*sin(p2 (x-r2)) dx */
double intfES(const double& p1, const double& r1,
              const double& p2, const double& r2,
              const double& x0, const double& x1);

/* integral_x0^x1 cos(p (x-r))*cos(p (x-r)) dx */
double intfCCeq(const double& p, const double& r,  
                const double& x0, const double& x1);

/* integral_x0^x1 cos(p (x-r))*sin(p (x-r)) dx */
double intfCSeq(const double& p, const double& r,  
                const double& x0, const double& x1);

/* integral_x0^x1 sin(p (x-r))*cos(p (x-r)) dx */
double intfSCeq(const double& p, const double& r,  
                const double& x0, const double& x1);

/* integral_x0^x1 sin(p (x-r))*sin(p (x-r)) dx */
double intfSSeq(const double& p, const double& r,  
                const double& x0, const double& x1);

/* integral_x0^x1 exp(p (x-r))*exp(p (x-r)) dx */
double intfEEpp(const double& p, const double& r,
                const double& x0, const double& x1);

/* integral_x0^x1 exp(p (x-r1))*exp(-p (x-r2)) dx */
double intfEEpn(const double& p, const double& r1, const double& r2,
                const double& x0, const double& x1);

/* integral_x0^x1 (a sin(p (x-r)) + b cos(p (x-r)))^2 dx */
double intfqOSC(const double& a, const double& b, 
                const double& p, const double& r, 
                const double& x0, const double& x1);

/* integral_x0^x1 (a exp(p (x-ra)) + b exp(-p (x-rb)))^2 dx */
double intfqHYP(const double& a, const double& b, const double& p, 
		const double& ra, const double& rb, 
                const double& x0, const double& x1);

/* integral_x0^x1 exp(-i k (x-r)) dx */
Complex intcexp(const Complex& k, const double& r, 
                const double& x0, const double& x1);

/* Number of zeroes in 
   a sin(g (x-x0)) + b cos(g (x-x0)) in the interval xl < x < xr */
int harmon_nodes(const double& a, const double& b, const double& g, 
                 const double& x0, const double& xl, const double& xr);

/* Number of zeroes in 
   a sin(g (x-x0)) + b cos(g (x-x0)) in the interval xl < x < xr;
   zl and zr are assigned the positions of the smallest and largest node,
   if present */
int harmon_nodes(const double& a, const double& b, const double& g, 
                 const double& x0, const double& xl, const double& xr,
                 double& zl, double& zr);

/* Number of zeroes in 
   a exp(g (x-xa)) + b exp(-g(x-xb)) in the interval xl < x < xr */
int exp_nodes(const double& a, const double& b, const double& g, 
              const double& xa, const double& xb, 
              const double& xl, const double& xr);

/* Number of zeroes in 
   a exp(g (x-xa)) + b exp(-g(x-xb)) in the interval xl < x < xr; 
   z is assigned the position of the node, if present */
int exp_nodes(const double& a, const double& b, const double& g, 
              const double& xa, const double& xb, 
              const double& xl, const double& xr, double& z);

/* maximum absolute level of  
   a sin(g (x-x0)) + b cos(g (x-x0)) in the interval xl< x < xr */
double harmon_max(const double& a, const double& b, const double& g, 
                  const double& x0, const double& xl, const double& xr);

/* maximum absolute level of  
   a exp(g (x-xa)) + b exp(-g(x-xb)) in the interval xl < x < xr */
double exp_max(const double& a, const double& b, const double& g, 
               const double& xa, const double& xb, 
               const double& xl, const double& xr);

/* --------------------------------------------------------------------- */

/* integral_x0^x1 (exp(i ga (x-ra))*exp(i gb (x-rb))) dx */
Complex intfcEcE(const Complex& ga, const double& ra,                  
                 const Complex& gb, const double& rb,
                 const double& x0, const double& x1);

/* maximum absolute level of  
   a exp(-i g (x-xa)) + b exp(i g(x-xb)) in the interval xl < x < xr,
   ph: phase at maximum position, 
   rough numerical evaluation */
double cexp_max(const Complex& a, const Complex& b, const Complex& g, 
                const double& xa, const double& xb, 
                const double& xl, const double& xr,
                double &ph);

#endif // INTEGRAL_H