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

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

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include "gengwed.h"
#include "integral.h"
#include "complex.h"

#define PQTolCC 1.0e-5
#define PQTolSS 1.0e-5
#define PQTolCS 1.0e-5
#define PQTolCE 1.0e-5
#define PQTolEE 1.0e-5

#define AWAY_2 5.0e+19       // smaller than [S_AWAY] in "structure.c"

/* error message */
void interror(char *m)
{
        fprintf(stderr, "\nintegral: %s.\n", m);
        exit(1);
}

/* 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)
{
	double pr, a, b, s1, s2;
	if(fabs(p1-p2)<PQTolCC)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p2+p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2+p1)*(r1-r2); 
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2-p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		s1 = cos(b+a)*sin(b-a)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCC)
	{
		if(p1 == -p2) 
			s2 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2); 
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2+p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		a = 0.5*(p1*(x0-r1)+p2*(x0-r2));
		b = 0.5*(p1*(x1-r1)+p2*(x1-r2));
		s2 = cos(b+a)*sin(b-a)/(p2+p1);
	}
	return s1+s2;
}

/* 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)
{
	double pr, a, b, s1, s2;
	if(fabs(p1-p2)<PQTolCS)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2-p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = -sin(a+b)*sin(a-b)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCS)
	{
		if(p1 == -p2)
			s2 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2+p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -sin(a+b)*sin(a-b)/(p2+p1);
	}
	return s1+s2;
}

/* 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)
{
	double pr, a, b, s1, s2;
	if(fabs(p1-p2)<PQTolCS)
	{
		if(p1 == p2)
			s1 = -0.25*(x1+x1-x0-x0)*sin(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = -0.25*(b-a)*sin(pr)
			     +0.0625*(p2-p1)*(b*b-a*a)*cos(pr);
		}	
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = sin(a+b)*sin(a-b)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCS)
	{
		if(p1 == -p2)
			s2 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2+p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -sin(a+b)*sin(a-b)/(p2+p1);
	}
	return s1+s2;
}

/* 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)
{
	double pr, a, b, s1, s2;
	if(fabs(p1-p2)<PQTolSS)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2-p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = cos(b+a)*sin(b-a)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolSS)
	{
		if(p1 == -p2)
			s2 = -0.25*(x1+x1-x0-x0)*cos(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = -0.25*(b-a)*cos(pr)
			     -0.0625*(p2+p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -cos(b+a)*sin(b-a)/(p2+p1);
	}
	return s1+s2;
}

/* 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)
{
	double a, b, s1, s2;
	if(fabs(p1+p2)<PQTolEE)
	{
		if(p1 == -p2)
			return 0.5*(x1+x1-x0-x0)*exp(0.5*(p2-p1)*(r1-r2)); 
		else
		{
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			return 0.5*(b-a+0.25*(p2+p1)*(b*b-a*a))
				  *exp(0.5*(p2-p1)*(r1-r2)); 
		}
	}
	else 
	{
		if(fabs(x1)>AWAY_2) s1 = 0.0;
		else s1 = exp(p1*(x1-r1)+p2*(x1-r2));
		if(fabs(x0)>AWAY_2) s2 = 0.0;
		else s2 = exp(p1*(x0-r1)+p2*(x0-r2));
		return (s1-s2)/(p1+p2);
	}
}

/* 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)
{
	double da2, db2, s1, s2;
	da2 = p2*(x0-r2);
	db2 = p2*(x1-r2);
	if(fabs(x1)>AWAY_2) s1 = 0.0;
	else s1 = exp(p1*(x1-r1))*(p2*sin(db2)+p1*cos(db2));	
	if(fabs(x0)>AWAY_2) s2 = 0.0;
	else s2 = exp(p1*(x0-r1))*(p2*sin(da2)+p1*cos(da2));	
	return (s1-s2)/(p1*p1+p2*p2);
}

/* 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)
{
	double da2, db2, s1, s2;
	da2 = p2*(x0-r2);
	db2 = p2*(x1-r2);
	if(fabs(x1)>AWAY_2) s1 = 0.0;
	else s1 = exp(p1*(x1-r1))*(p1*sin(db2)-p2*cos(db2));	
	if(fabs(x0)>AWAY_2) s2 = 0.0;
	else s2 = exp(p1*(x0-r1))*(p1*sin(da2)-p2*cos(da2));	
	return (s1-s2)/(p1*p1+p2*p2);
}

/* 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)
{
	double da1, db1, s1, s2;
	da1 = p1*(x0-r1);
	db1 = p1*(x1-r1);
	if(fabs(x1)>AWAY_2) s1 = 0.0;
	else s1 = exp(p2*(x1-r2))*(p1*sin(db1)+p2*cos(db1));	
	if(fabs(x0)>AWAY_2) s2 = 0.0;
	else s2 = exp(p2*(x0-r2))*(p1*sin(da1)+p2*cos(da1));	
	return (s1-s2)/(p1*p1+p2*p2);
}

/* 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)
{
	double da1, db1, s1, s2;
	da1 = p1*(x0-r1);
	db1 = p1*(x1-r1);
	if(fabs(x1)>AWAY_2) s1 = 0.0;
	else s1 = exp(p2*(x1-r2))*(p2*sin(db1)-p1*cos(db1));	
	if(fabs(x0)>AWAY_2) s2 = 0.0;
	else s2 = exp(p2*(x0-r2))*(p2*sin(da1)-p1*cos(da1));	
	return (s1-s2)/(p1*p1+p2*p2);
}

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

/* 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)
{
	double a, b, s1, s2;
	s1 = 0.25*(x1+x1-x0-x0);
	a = p*(x0-r);
	b = p*(x1-r);
	s2 = 0.5*cos(b+a)*sin(b-a)/p;
	return s1+s2;
}

/* 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)
{
	double a, b;
	b = p*(x1-r);
	a = p*(x0-r);
	return -0.5*sin(a+b)*sin(a-b)/p;
}

/* 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)
{
	double a, b;
	b = p*(x1-r);
	a = p*(x0-r);
	return -0.5*sin(a+b)*sin(a-b)/p;
}

/* 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)
{
	double a, b, s1, s2;
	s1 = 0.25*(x1+x1-x0-x0);
	b = p*(x1-r);
	a = p*(x0-r);
	s2 = -0.5*cos(b+a)*sin(b-a)/p;
	return s1+s2;
}

/* 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)
{
	double s1, s2;
	if(fabs(x1)>AWAY_2) s1 = 0.0;
	else s1 = exp(2.0*p*(x1-r));
	if(fabs(x0)>AWAY_2) s2 = 0.0;
	else s2 = exp(2.0*p*(x0-r));
	return 0.5*(s1-s2)/p;
}

/* 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)
{
	return (x1-x0)*exp(p*(r2-r1)); 
}

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

/* 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)
{
	double s;
	s  =     a*a*intfSSeq(p, r, x0, x1); 
	s += 2.0*a*b*intfSCeq(p, r, x0, x1); 
	s +=     b*b*intfCCeq(p, r, x0, x1); 
	return s;
} 

/* 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)
{
	double s;
	s  =     a*a*intfEEpp(p, ra, x0, x1); 
	s += 2.0*a*b*intfEEpn(p, ra, rb, x0, x1); 
	s +=     b*b*intfEEpp(-p, rb, x0, x1); 
	return s;
} 

/* integral_x0^x1 exp(-i k (x-r)) dx */
Complex intcexp(const Complex& k, const double& r, 
                const double& x0, const double& x1)
{
	double p = k.re;
	double q = k.im;
	return Complex( intfEC(q, r, p, r, x0, x1),
	               -intfES(q, r, p, r, x0, 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)
{
	double phi=0.0;
//	double amp=0.0;
	double ni, nii;

	if(a == 0.0 && b == 0.0) return 10000; 
	if(a == 0.0)
	{
		phi = PI/2.0;
//		amp = b;
	}
	else
	{
//		amp = sqrt(a*a+b*b);
		phi = atan(b/a);
		if(a < 0.0) phi += PI;
	}
//	if(fabs(a/amp - cos(phi))> 1.0e-12 || fabs(b/amp - sin(phi))>1.0e-12 )
//		interror("harmon_nodes: decomposition");
	ni  =  ceil((g*(xl-x0)+phi)/PI);
	nii = floor((g*(xr-x0)+phi)/PI);
	return ((int) nii) - ((int) ni) + 1;
}

/* 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) 
{
	double phi=0.0;
//	double amp=0.0;
	double ni, nii;
	int nz;

	if(a == 0.0 && b == 0.0) return 10000; 
	if(a == 0.0)
	{
		phi = PI/2.0;
//		amp = b;
	}
	else
	{
//		amp = sqrt(a*a+b*b);
		phi = atan(b/a);
		if(a < 0.0) phi += PI;
	}
//	if(fabs(a/amp - cos(phi))> 1.0e-12 || fabs(b/amp - sin(phi))>1.0e-12 )
//		interror("harmon_nodes: decomposition");
	ni  =  ceil((g*(xl-x0)+phi)/PI);
	nii = floor((g*(xr-x0)+phi)/PI);
	nz = ((int) nii) - ((int) ni) + 1;
	if(nz >= 1) 
	{
		zl = (ni *PI-phi)/g+x0; 
		zr = (nii*PI-phi)/g+x0; 
	}	
	return nz;
}

/* 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)
{
	double x;
	if(a == 0.0 && b == 0.0) return 10000;
	if(a == 0.0 || b == 0.0) return 0;
	if(a*b > 0.0) return 0;
	x = 0.5*(xb+xa+log(-b/a)/g);
	if(xl<x && x<xr) return 1;
	return 0;
}

/* 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)
{
	double x;
	if(a == 0.0 && b == 0.0) return 10000;
	if(a == 0.0 || b == 0.0) return 0;
	if(a*b > 0.0) return 0;
	x = 0.5*(xb+xa+log(-b/a)/g);
	if(xl<x && x<xr) { z = x; return 1; }
	return 0;
}

/* 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) 
{
	double phi=0.0;
//	double amp=0.0;
	int pl, pr;
	double p;

	if(a==0.0 && b==0.0) return 0.0;
	if(a == 0.0)
	{
		phi = PI/2.0;
//		amp = b;
	}
	else
	{
//		amp = sqrt(a*a+b*b);
		phi = atan(b/a);
		if(a < 0.0) phi += PI;
	}
//	if(fabs(a/amp - cos(phi))> 1.0e-12 || fabs(b/amp - sin(phi))>1.0e-12 )
//		interror("harmon_max: decomposition");
	pl = (int) ceil((g*(xl-x0)+phi)/PI-0.5);
	pr = (int) floor((g*(xr-x0)+phi)/PI-0.5);
	double m, mt; 
	m  = fabs(a*sin(g*(xl-x0))+b*cos(g*(xl-x0)));
	for(int l=pl; l<=pr; ++l)
	{
		p  = ((((double) l) + 0.5)*PI-phi)/g+x0;
		mt = fabs(a*sin(g*(p-x0))+b*cos(g*(p-x0)));
		if(mt > m) m = mt;
	}
	mt = fabs(a*sin(g*(xr-x0))+b*cos(g*(xr-x0)));
	if(mt > m) m = mt;
	return m;
}

/* 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)
{
	double x;
	if(a==0.0 && b==0.0) return 0.0;
	double m, mt;
	m  = fabs(a*exp(g*(xl-xa))+b*exp(g*(xb-xl)));
	mt = fabs(a*exp(g*(xr-xa))+b*exp(g*(xb-xr)));
	if(mt > m) m = mt;
	if(a == 0.0 || b == 0) return m;
	if(a*b > 0.0)
	{
		x = 0.5*(xb+xa+log(b/a)/g);
		if(xl<x && x<xr)
		{
			mt = fabs(a*exp(g*(x-xa))+b*exp(g*(xb-x)));
			if(mt > m) m = mt;
		}
	}
	return m;
}

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

/* 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)
{
	double e;
	Complex s1, s2;
	if((ga+gb).abs()<PQTolEE)
	{
		if((ga+gb).abs()<1.0e-15)
			return (x1-x0)*expi(0.5*(ga-gb)*(rb-ra)); 
		else
			return expi(0.5*(ga-gb)*(rb-ra))*(x1-x0+0.5*CCI*(ga+gb)*(x1*x1-x0*x0-(x1-x0)*(ra+rb)));
	}
	else 
	{
		e = (ga+gb).re*x1;
		if(e > 0 && fabs(x1)>AWAY_2) s1 = 0.0;
		else                         s1 = expi(ga*(x1-ra)+gb*(x1-rb));
		e = (ga+gb).re*x0;
		if(e > 0 && fabs(x0)>AWAY_2) s2 = 0.0;
		else                         s2 = expi(ga*(x0-ra)+gb*(x0-rb));
		return CCmI*(s1-s2)/(ga+gb);
	}
}

/* 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 the position of maximum field,
   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)
{
	double m = 0.0;
	Complex f;
	double s;
	if(fabs(g.re) > 0.0)
	{ 
		s = 2.0*PI/fabs(g.re)/20.0;
		if((xr-xl)/s < 10.0) s = (xr-xl)/10.0;
		if((xr-xl)/s > 2000.0) s = (xr-xl)/2000.0;
	}	
	else 
	{
		s = (xr-xl)/100.0;
	}
	f = a*expmi(g*(xl-xa)) + b*expi(g*(xl-xb));
	m = f.abs();
	ph = f.arg();
	double x = xl+s;
	double mt;
	while(x < xr)
	{
		f = a*expmi(g*(x-xa)) + b*expi(g*(x-xb));
		mt = f.abs();
		if(mt > m)
		{
			m = mt;
			ph = f.arg();		
		} 
		x += s;
	}
	f = a*expmi(g*(xr-xa)) + b*expi(g*(xr-xb));
	mt = f.abs();
	if(mt > m) 
	{
		m = mt;
		ph = f.arg();
	}
	return m;	
}