/*
 * WMM - Wave-matching method for mode analysis of dielectric waveguides
 * https://wmm.computational-photonics.eu/
 */

/*
 * wmmfld.cpp
 * WMM_Mode - procedures for evaluating the mode profile
 */

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmtrif.h"
#include"mode.h"
#include"wmmfld.h"

/* minimum distance between boundaries to be considered separated */
#define HDIST 1.0e-8

/* sum up trial functions */
double WMM_Mode::fieldsum(int cp, Derlev abl, int l, int m, double x, double y)
const
{
	int j;
	double f = 0.0;
	for(j=0; j<=phi.ntf(cp,l,m)-1; ++j)
		f += (*phi(cp,l,m,j)).val(abl, x, y);
	return f;
}

/* product of field values at point (x,y); cp1, cp2 != Sz */
double WMM_Mode::fieldprod(Fcomp cp1, Fcomp cp2, double x, double y) const
{
	int l, m, j; 
	double s1, s2, sp;
	Fieldrep id; 
	double nlm;
	wg.rectidx(x, y, l, m);
	nlm = wg.n(l,m);
	id = getfieldrep(cp1, pol, vform);
	s1 = 0.0;
	for(j=0; j<=id.numc-1; ++j)
	{
		sp = fieldsum(id.ctr[j].c, id.ctr[j].fod, l, m, x, y);
		sp *= evalvfac(id.ctr[j].sfid, beta, k0, nlm)*id.ctr[j].fac;
		s1 += sp;
	}
	s1 *= evalvfac(id.glsfid, beta, k0, nlm)*id.glfac;
	id = getfieldrep(cp2, pol, vform);
	s2 = 0.0;
	for(j=0; j<=id.numc-1; ++j)
	{
		sp = fieldsum(id.ctr[j].c, id.ctr[j].fod, l, m, x, y);
		sp *= evalvfac(id.ctr[j].sfid, beta, k0, nlm)*id.ctr[j].fac;
		s2 += sp;
	}
	s2 *= evalvfac(id.glsfid, beta, k0, nlm)*id.glfac;
	return ampfac*ampfac*s1*s2;
}

/* field values at point (x,y) */
double WMM_Mode::field(Fcomp cp, double x, double y) const
{
	int l, m, j; 
	double s, sp;
	Fieldrep id; 
	double nlm;
	wg.rectidx(x, y, l, m);
	nlm = wg.n(l,m);
	if(cp == SZ)
	{
		return 0.5*(fieldprod(EX, HY, x, y)
			   -fieldprod(HX, EY, x, y));
	}
	id = getfieldrep(cp, pol, vform);
	s = 0.0;
	for(j=0; j<=id.numc-1; ++j)
	{
		sp = fieldsum(id.ctr[j].c, id.ctr[j].fod, l, m, x, y);
		sp *= evalvfac(id.ctr[j].sfid, beta, k0, nlm)*id.ctr[j].fac;
		s += sp;
	}
	s *= evalvfac(id.glsfid, beta, k0, nlm)*id.glfac;
	return ampfac*s;
}

/* ##################################################################### */

#define LHOR 0
#define LVER 1

/* integration of field products along lines */

/* helper function */
double WMM_Mode::lint(int dir, Derlev fod0, int cp0, Derlev fod1, int cp1, 
		      int l, int m, double a, double b, double c) const
{
	double sum;
	int j, k;
	int jrj, jrk;

	sum = 0.0;
	jrj = phi.ntf(cp0,l,m);
	for(j=0; j<=jrj-1; ++j)
	{
		jrk = phi.ntf(cp1,l,m);
		for(k=0; k<=jrk-1; ++k)
		{
			switch(dir)
			{
				case LHOR: 
					sum += tfinty(*phi(cp0,l,m,j), fod0, 
					  	      *phi(cp1,l,m,k), fod1, 
						      c, a, b); 
					break;
				case LVER: 
					sum += tfintx(*phi(cp0,l,m,j), fod0, 
						      *phi(cp1,l,m,k), fod1, 
						      c, a, b); 
					break;
				default: break;
			}
		}
	}
	return ampfac*ampfac*sum;
}

/* helper function */
double WMM_Mode::lineint(int dir, int l, int m, Fcomp cp1, Fcomp cp2, 
	           double a, double b, double c) const
{
	double s, sp;
	int j, k;
	Fieldrep idj, idk; 
	double vz, t;
	double nlm;
	nlm = wg.n(l,m);
	vz = 1.0; if(a>b) {t = a; a = b; b = t; vz = -1.0;
fprintf(stderr, "lineint: vz !\n"); 
			}

	idj = getfieldrep(cp1, pol, vform);
	idk = getfieldrep(cp2, pol, vform);
	s = 0.0;
	for(j=0; j<=idj.numc-1; ++j)
	{
		for(k=0; k<=idk.numc-1; ++k)
		{
			sp = lint(dir, idj.ctr[j].fod, idj.ctr[j].c,
			               idk.ctr[k].fod, idk.ctr[k].c, l, m,
				       a, b, c);

			sp *= evalvfac(idj.ctr[j].sfid, beta, k0, nlm)
			      *idj.ctr[j].fac;
			sp *= evalvfac(idk.ctr[k].sfid, beta, k0, nlm)
			      *idk.ctr[k].fac;
			s += sp;
		}
	}
	s *= evalvfac(idj.glsfid, beta, k0, nlm)*idj.glfac;
	s *= evalvfac(idk.glsfid, beta, k0, nlm)*idk.glfac;
// ORG: return s;
	return vz*s;
}

/* integration of field products along horizontal lines */
double WMM_Mode::horlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
		  double x, double ya, double yb) const
{
	double i;
	i = lineint(LHOR, l, m, cp1, cp2, ya, yb, x);
//	fprintf(stderr, "[%d,%d] <%c%c,%c%c> (%g) %g -> %g = %g\n",
//		l, m, fldchr(cp1), cpchr(cp1), fldchr(cp2), cpchr(cp2), 
//		x, ya, yb, i);
	return i;
}
/* integration of field products along vertical lines */
double WMM_Mode::verlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
		  double xa, double xb, double y) const
{
	return lineint(LVER, l, m, cp1, cp2, xa, xb, y);
}

/* ##################################################################### */

/* integrations over rectangles ... */

/* helper function */
double WMM_Mode::pint(Derlev fod0, int cp0, Derlev fod1, int cp1, int l, int m,
		      Rect r) const
{
	double sum;
	int j, k;
	int jrj, jrk;

	sum = 0.0;
	jrj = phi.ntf(cp0,l,m);
	for(j=0; j<=jrj-1; ++j)
	{
		jrk = phi.ntf(cp1,l,m);
		for(k=0; k<=jrk-1; ++k)
			sum += tfnrm(*phi(cp0,l,m,j), fod0, 
                                     *phi(cp1,l,m,k), fod1, r);
	}
	return ampfac*ampfac*sum;
}

/* helper function */
double WMM_Mode::prodint(int l, int m, Fcomp cp1, Fcomp cp2, Rect r) const
{
	double s, sp;
	int j, k;
	Fieldrep idj, idk; 
	double nlm;
	nlm = wg.n(l,m);
	idj = getfieldrep(cp1, pol, vform);
	idk = getfieldrep(cp2, pol, vform);
	s = 0.0;
	for(j=0; j<=idj.numc-1; ++j)
	{
		for(k=0; k<=idk.numc-1; ++k)
		{
			sp = pint(idj.ctr[j].fod, idj.ctr[j].c,
			          idk.ctr[k].fod, idk.ctr[k].c, l, m, r);

			sp *= evalvfac(idj.ctr[j].sfid, beta, k0, nlm)
			      *idj.ctr[j].fac;
			sp *= evalvfac(idk.ctr[k].sfid, beta, k0, nlm)
			      *idk.ctr[k].fac;
			s += sp;
		}
	}
	s *= evalvfac(idj.glsfid, beta, k0, nlm)*idj.glfac;
	s *= evalvfac(idk.glsfid, beta, k0, nlm)*idk.glfac;
	return s;
}

/* integration of a fieldproduct over the entire rectangle l, m */
double WMM_Mode::recintlm(int l, int m, Fcomp cp1, Fcomp cp2) const
{
	Rect r;
	r = wg.rectbounds(l,m);
	return prodint(l, m, cp1, cp2, r);
}

/* integration of a fieldproduct over an arbitrary rectangle */
double WMM_Mode::recint(Fcomp cp1, Fcomp cp2, Rect r) const
{
	double x, y;
	double xt, yr;
	double s;
	int l, m, dum;
	Rect rp;
	if(r.x0 > r.x1) { x=r.x1; r.x1=r.x0; r.x0=x; }
	if(r.y0 > r.y1) { y=r.y1; r.y1=r.y0; r.y0=y; }
	s = 0.0;
	y = r.y0;
	while(fabs(r.y1-y) > HDIST)
	{
		wg.rectidx(0.0, y+HDIST/2.0, dum, m);
		yr = wg.rectbounds(0,m).y1;
		if(r.y1<yr) yr = r.y1;
		x = r.x0;	
		while(fabs(r.x1-x) > HDIST)
		{
			wg.rectidx(x+HDIST/2.0, 0.0, l, dum);
			xt = wg.rectbounds(l,0).x1;
			if(r.x1<xt) xt = r.x1;
			rp.x0 = x; rp.y0 = y; rp.x1 = xt; rp.y1 = yr;
			s += prodint(l, m, cp1, cp2, rp); 
			x = xt;
		}
		y = yr;
	}
	return s;
}

/* translate mode by (dx, dy) */
void WMM_Mode::translate(double dx, double dy)
{
	wg.translate(dx, dy);
	phi.translate(dx, dy);
	sym = NOS;
	return;
}

/* output to file, representation depending part of WMM_Mode */
void WMM_Mode::writerdpid(FILE *dat)
{
	fputchar('W', dat);
	fputchar('M', dat);
	fputchar('M', dat);
	return;
}
void WMM_Mode::writerdp(FILE *dat)
{
	comment("WMM_MODE", dat);
	comment("vform", dat);
	fputint(vform, dat);
	phi.write(dat);
	return;
}

/* input from file, representation depending part of WMM_Mode */
void WMM_Mode::readrdpid(FILE *dat)
{
	int c1, c2, c3;
	c1 = fgetchar(dat);
	c2 = fgetchar(dat);
	c3 = fgetchar(dat);
	if(c1 != 'W' || c2 != 'M' || c3 != 'M')
	{
		fprintf(stderr, 
		"WMM_Mode: readrdpid: bad data file ! (%c,%c,%c)\n", c1,c2,c3);
		exit(1);
	}
	return;
}
void WMM_Mode::readrdp(FILE *dat)
{
	vform = Vecform(fgetint(dat));
	phi.read(dat);
	return;
}

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

/* scalar product 0.5\int\int(E_1x H_2y - E_1y H_2x)dxdy of two modes : */

/* helper function */
double twompint(const WMM_Mode& mode0, Derlev fod0, int cp0, int l0, int m0, 
		const WMM_Mode& mode1, Derlev fod1, int cp1, int l1, int m1,
		Rect r)
{
	double sum;
	int j, k;
	int jrj, jrk;

	sum = 0.0;
	jrj = mode0.phi.ntf(cp0,l0,m0);
	for(j=0; j<=jrj-1; ++j)
	{
		jrk = mode1.phi.ntf(cp1,l1,m1);
		for(k=0; k<=jrk-1; ++k)
			sum += tfnrm(*(mode0.phi)(cp0,l0,m0,j), fod0, 
                                     *(mode1.phi)(cp1,l1,m1,k), fod1, r);
	}
	return mode0.ampfac*mode1.ampfac*sum;
}

/* helper function */
double twomprodint(const WMM_Mode& mode1, int l1, int m1, Fcomp cp1, 
	           const WMM_Mode& mode2, int l2, int m2, Fcomp cp2,
                   Rect r)
{
	double s, sp;
	int j, k;
	Fieldrep idj, idk; 
	double nlm1, nlm2;
	nlm1 = mode1.wg.n(l1,m1);
	nlm2 = mode2.wg.n(l2,m2);
	idj = getfieldrep(cp1, mode1.pol, mode1.vform);
	idk = getfieldrep(cp2, mode2.pol, mode2.vform);
	s = 0.0;
	for(j=0; j<=idj.numc-1; ++j)
	{
		for(k=0; k<=idk.numc-1; ++k)
		{
			sp = twompint(mode1, 
				      idj.ctr[j].fod, idj.ctr[j].c, l1, m1, 
			              mode2, 
				      idk.ctr[k].fod, idk.ctr[k].c, l2, m2, r);

			sp *= evalvfac(idj.ctr[j].sfid, 
				       mode1.beta, mode1.k0, nlm1)
			      *idj.ctr[j].fac;
			sp *= evalvfac(idk.ctr[k].sfid, 
				       mode2.beta, mode2.k0, nlm2)
			      *idk.ctr[k].fac;
			s += sp;
		}
	}
	s *= evalvfac(idj.glsfid, mode1.beta, mode1.k0, nlm1)*idj.glfac;
	s *= evalvfac(idk.glsfid, mode2.beta, mode2.k0, nlm2)*idk.glfac;
	return s;
}

/* integration of a fieldproduct between two modes 
   over an arbitrary rectangle */
double twomrecint(const WMM_Mode& mode1, Fcomp cp1, 
                  const WMM_Mode& mode2, Fcomp cp2, Rect r)
{
	double x, y;
	double xt, yr;
	double yr1, yr2;
	double xt1, xt2;
	double s, z;
	int l1, m1;
	int l2, m2;
	int dum;
	Rect rp;
	if(r.x0 > r.x1) { x=r.x1; r.x1=r.x0; r.x0=x; }
	if(r.y0 > r.y1) { y=r.y1; r.y1=r.y0; r.y0=y; }
	rp = mode1.wg.xyplane;
	if(r.x0 < rp.x0) r.x0=rp.x0;
	if(r.x1 > rp.x1) r.x1=rp.x1;
	if(r.y0 < rp.y0) r.y0=rp.y0;
	if(r.y1 > rp.y1) r.y1=rp.y1;
	rp = mode2.wg.xyplane;
	if(r.x0 < rp.x0) r.x0=rp.x0;
	if(r.x1 > rp.x1) r.x1=rp.x1;
	if(r.y0 < rp.y0) r.y0=rp.y0;
	if(r.y1 > rp.y1) r.y1=rp.y1;
	s = 0.0;
	y = r.y0;
	while(fabs(r.y1-y) > HDIST)
	{
		mode1.wg.rectidx(0.0, y+HDIST/2.0, dum, m1);
		mode2.wg.rectidx(0.0, y+HDIST/2.0, dum, m2);
		yr1 = mode1.wg.rectbounds(0,m1).y1;
		yr2 = mode2.wg.rectbounds(0,m2).y1;
		yr = yr1; 
		if(yr2<yr) yr = yr2;
		if(r.y1<yr) yr = r.y1;
		x = r.x0;	
		while(fabs(r.x1-x) > HDIST)
		{
			mode1.wg.rectidx(x+HDIST/2.0, 0.0, l1, dum);
			mode2.wg.rectidx(x+HDIST/2.0, 0.0, l2, dum);
			xt1 = mode1.wg.rectbounds(l1,0).x1;
			xt2 = mode2.wg.rectbounds(l2,0).x1;
			xt = xt1; 
			if(xt2<xt) xt = xt2;
			if(r.x1<xt) xt = r.x1;
			rp.x0 = x; rp.y0 = y; rp.x1 = xt; rp.y1 = yr;
			z = twomprodint(mode1, l1, m1, cp1,
			    	        mode2, l2, m2, cp2, rp);
			s += z;
			x = xt;
		}
		y = yr;
	}
	return s;
}

/* scalar product 0.5\int\int(E_1x H_2y - E_1y H_2x)dxdy of two modes */
double scalprod(const WMM_Mode& m1, const WMM_Mode& m2) 
{
	Rect r;
	r = m1.wg.xyplane; 	  
	if(m2.wg.xyplane.x0 > r.x0) r.x0 = m2.wg.xyplane.x0;
	if(m2.wg.xyplane.y0 > r.y0) r.y0 = m2.wg.xyplane.y0;
	if(m2.wg.xyplane.x1 < r.x1) r.x1 = m2.wg.xyplane.x1;
	if(m2.wg.xyplane.y1 < r.y1) r.y1 = m2.wg.xyplane.y1;
	if(r.x1 <= r.x0) return 0.0;
	if(r.y1 <= r.y0) return 0.0;
/*
	int l, m;
	double s;
	if(m1.pol == QTE && m2.pol == QTE)
		return 0.5*m2.beta*val_invommu0(m1.wg.lambda)
		       *twomrecint(m1, EY, m2, HY, r);	

	if(m1.pol == QTM && m2.pol == QTM)
	{
		s = 0.0;
		for(l=0; l<=m2.wg.nx+1; ++l)
		{
			for(m=0; m<=m2.wg.ny+1; ++m)
				s += twomrecint(m1, HY, m2, HY, 
				     m2.wg.rectbounds(l, m))/m2.wg.eps(l,m);
		}
		return 0.5*m2.beta*val_invomep0(m2.wg.lambda)*s;
	}
*/

	return 0.5*( twomrecint(m1, EX, m2, HY, r) 	
               	    -twomrecint(m1, EY, m2, HX, r));
}

/* scalar product 
   0.25\int\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dxdy 
   of two modes */
double lscalprod(const WMM_Mode& m1, const WMM_Mode& m2) 
{
	return 0.5*(scalprod(m1,m2)+scalprod(m2,m1));
}

/* for a superposition of two WMM_Modes m1, m2 with coefficients c1, c2:
   intensity behind a polarizer at an angle alpha versus TE-polarization */
double polarizeroutput(Complex c1, const WMM_Mode& m1, 
                       Complex c2, const WMM_Mode& m2, 
		       double alpha)
{
	double p, pt;
	p = 0.0;

	pt = 0.0;
	if(sin(alpha) != 0.0)
	    pt +=     sin(alpha)*sin(alpha)*m1.recint(EX, EX, m1.wg.xyplane); 
	if(cos(alpha) != 0.0)
	    pt +=     cos(alpha)*cos(alpha)*m1.recint(EY, EY, m1.wg.xyplane); 
	if(sin(alpha) != 0.0 && cos(alpha) != 0.0)
	    pt += 2.0*sin(alpha)*cos(alpha)*m1.recint(EX, EY, m1.wg.xyplane); 
	pt *= c1.sqabs();
	p += pt;

	pt = 0.0;
	if(sin(alpha) != 0.0)
	    pt +=     sin(alpha)*sin(alpha)*m2.recint(EX, EX, m2.wg.xyplane); 
	if(cos(alpha) != 0.0)
	    pt +=     cos(alpha)*cos(alpha)*m2.recint(EY, EY, m2.wg.xyplane); 
	if(sin(alpha) != 0.0 && cos(alpha) != 0.0)
	    pt += 2.0*sin(alpha)*cos(alpha)*m2.recint(EX, EY, m2.wg.xyplane); 
	pt *= c2.sqabs(); 
	p += pt;

	pt = 0.0;
	if(sin(alpha) != 0.0)
	 pt += sin(alpha)*sin(alpha)*twomrecint(m1, EX, m2, EX, m1.wg.xyplane); 
	if(cos(alpha) != 0.0)
	 pt += cos(alpha)*cos(alpha)*twomrecint(m1, EY, m2, EY, m1.wg.xyplane); 
	if(sin(alpha) != 0.0 && cos(alpha) != 0.0)
	 pt += sin(alpha)*cos(alpha)*twomrecint(m1, EX, m2, EY, m1.wg.xyplane); 
	if(sin(alpha) != 0.0 && cos(alpha) != 0.0)
	 pt += sin(alpha)*cos(alpha)*twomrecint(m1, EY, m2, EX, m1.wg.xyplane); 
	pt *= 2.0*(c1.re*c2.re + c1.im*c2.im);
	p += pt;

	p *= 0.5*INVSQRTMU0/INVSQRTEP0;
	return p;
}