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

/* 
 * tfln.cpp
 * Three layer TFLN slab waveguides, anisotropic core, 
 * X- and Z-cut configurations with varying propagation direction
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"matlvis.h"
#include"tfln.h"

/* error message */
void tflnerror(const char *s)
{
        fprintf(stderr, "\ntfln: %s.\n", s);
        exit(1);
}

/* --------------------------------------------------------------- 
   TFLN slab waveguide */

/* get interval index corresponding to position x */
int TflnWaveguide::layeridx(double x) const
{
	int l=0;
	while(l<=1 && x>hx(l)) ++l;
	return l;
}
/* get layer boundaries corresponding to index l */
Interval TflnWaveguide::layer(int l) const
{
	Interval i;
	if(l <= 0) i.x0 = -AWAY;
	else i.x0 = hx(l-1);
	if(l >= 2) i.x1 = AWAY;
	else i.x1 = hx(l);
	return i;
}

/* get layer boundaries corresponding to position x */ 
Interval TflnWaveguide::layer(double x) const
{
	return layer(layeridx(x));
}

/* initialize */
TflnWaveguide::TflnWaveguide()
{
	type = Xcut;
	hx = Dvector(2);
	hx(0) = 0.0; 
	hx(1) = 1.0;
	ns = 1.0;
	no = 2.0;
	ne = 2.0;
	nc = 1.0;
	theta = 0.0;
	lambda = 1.0;
	turnto(0.0);
}
TflnWaveguide::TflnWaveguide(TflnCut c, double d, double tns, double tno, double tne, double tnc, double wavel)
{
	type = c;
	hx = Dvector(2);
	hx(0) = 0.0; 
	hx(1) = d;
	ns = tns;
	no = tno;
	ne = tne;
	nc = tnc;
	theta = 0.0;
	lambda = wavel;
	turnto(0.0);
}

/* adjust permittivity for propagation at angle ttheta [deg] versus crystal Y-axis */
void TflnWaveguide::turnto(double ttheta)
{
	theta = ttheta;
	double ct = cos(theta*PI/180.0);
	double st = sin(theta*PI/180.0);
	if(type == Xcut)
	{
		epsx = no*no;
		epsy = ne*ne*ct*ct + no*no*st*st;
		epsz = ne*ne*st*st + no*no*ct*ct;
		delta = ct*st*(ne*ne-no*no);
		return;
	}
	// Zcut
	epsx = ne*ne;
	epsy = no*no;
	epsz = no*no;
	delta = 0.0;
	return;
}

/* an isotropic waveguide with the same layering */
Waveguide TflnWaveguide::toisotropic() const
{
	Waveguide g(1);
	g.hx = hx;
	g.n(0) = ns;
	g.n(1) = 0.5*(no+ne);
	g.n(2) = nc;
	g.lambda = lambda;
	return g;
}

/* --------------------------------------------------------------- 
   a mode of the TFLN slab */

/* initialize */
TflnMode::TflnMode()
{
	wg = TflnWaveguide();
	k0 = val_k0(wg.lambda);
	invommu0 = val_invommu0(wg.lambda);
	invomep0 = val_invomep0(wg.lambda);
	invsqrtmu0 = 1.0/sqrt(MU0);
	invsqrtep0 = 1.0/sqrt(EP0);
}
TflnMode::TflnMode(const TflnWaveguide& g)
{
	wg = g;
	k0 = val_k0(wg.lambda);
	invommu0 = val_invommu0(wg.lambda);
	invomep0 = val_invomep0(wg.lambda);
	invsqrtmu0 = 1.0/sqrt(MU0);
	invsqrtep0 = 1.0/sqrt(EP0);
}

/* principal field components */
Complex TflnMode::pfe(double x) const
{
	int l = wg.layeridx(x);
	switch(l)
	{
	case 0: 
		return FAes*exp(kappas*x);
		break;
	case 1: 
		return FAa1*CVe1*expmi(kappa1*x)
		      +FAa2*CVe2*expmi(kappa2*x)
		      +FAa3*CVe3*expmi(kappa3*x)
		      +FAa4*CVe4*expmi(kappa4*x);
		break;
	case 2: 
		return FAec*exp(-kappac*(x-(wg.hx(1)-wg.hx(0))));
		break;
	default: 
		tflnerror("TflnMode, pfe, l");
		return CC0;
		break;
	}
}
Complex TflnMode::pfh(double x) const
{
	int l = wg.layeridx(x);
	switch(l)
	{
	case 0: 
		return FAhs*exp(kappas*x);
		break;
	case 1: 
		return FAa1*CVh1*expmi(kappa1*x)
		      +FAa2*CVh2*expmi(kappa2*x)
		      +FAa3*CVh3*expmi(kappa3*x)
		      +FAa4*CVh4*expmi(kappa4*x);
		break;
	case 2: 
		return FAhc*exp(-kappac*(x-(wg.hx(1)-wg.hx(0))));
		break;
	default: 
		tflnerror("TflnMode, pfh, l");
		return CC0;
		break;
	}
}
/* derivatives of the principal field components */
Complex TflnMode::pfde(double x) const
{
	int l = wg.layeridx(x);
	switch(l)
	{
	case 0: 
		return kappas*FAes*exp(kappas*x);
		break;
	case 1: 
		return FAa1*CCmI*kappa1*CVe1*expmi(kappa1*x)
		      +FAa2*CCmI*kappa2*CVe2*expmi(kappa2*x)
		      +FAa3*CCmI*kappa3*CVe3*expmi(kappa3*x)
		      +FAa4*CCmI*kappa4*CVe4*expmi(kappa4*x);
		break;
	case 2: 
		return -kappac*FAec*exp(-kappac*(x-(wg.hx(1)-wg.hx(0))));
		break;
	default: 
		tflnerror("TflnMode, pfde, l");
		return CC0;
		break;
	}
}
Complex TflnMode::pfdh(double x) const
{
	int l = wg.layeridx(x);
	switch(l)
	{
	case 0: 
		return kappas*FAhs*exp(kappas*x);
		break;
	case 1: 
		return FAa1*CCmI*kappa1*CVh1*expmi(kappa1*x)
		      +FAa2*CCmI*kappa2*CVh2*expmi(kappa2*x)
		      +FAa3*CCmI*kappa3*CVh3*expmi(kappa3*x)
		      +FAa4*CCmI*kappa4*CVh4*expmi(kappa4*x);
		break;
	case 2: 
		return -kappac*FAhc*exp(-kappac*(x-(wg.hx(1)-wg.hx(0))));
		break;
	default: 
		tflnerror("TflnMode, pfdh, l");
		return CC0;
		break;
	}
}

/* the mode profile, 
   values of component cp: EX - HZ at transverse coordinate x */
Complex TflnMode::field(Fcomp cp, double x) const
{
	Complex e = pfe(x);
	Complex h = pfh(x);
	Complex de = pfde(x);
	Complex dh = pfdh(x);
	double epx, epz, del;
	int l = wg.layeridx(x);
	switch(l) 
	{
	case 0: 
		epx = wg.ns*wg.ns;
		epz = epx;
		del = 0.0;
		break;
	case 1: 
		epx = wg.epsx;
		epz = wg.epsz;
		del = wg.delta;
		break;
	case 2: 
		epx = wg.nc*wg.nc;
		epz = epx;
		del = 0.0;
		break;
	default: tflnerror("TflnMode, field, l");
		break;
	}	
	switch(cp)
	{
	case EX: 
		return beta*invomep0*invsqrtmu0*h/epx;
		break;
	case EY: 
		return invsqrtep0*e;
		break;
	case EZ: 
		return CCmI*invomep0*invsqrtmu0*dh/epz - del*invsqrtep0*e/epz;
		break;
	case HX: 
		return -beta*invommu0*invsqrtep0*e;
		break;
	case HY: 
		return invsqrtmu0*h;
		break;
	case HZ: 
		return CCI*invommu0*invsqrtep0*de;
		break;
		
	default: tflnerror("TflnMode, field, cp");
		break;
	}
	return CC0;
}

/* adjust global phase of the mode, such that, 
   at the point in the core with maximum electric field strength, 
   the largest transverse electric field component becomes real */
void TflnMode::adjustphase()
{
	double x0 = wg.hx(0)+HDIST;
	double x1 = wg.hx(1)-HDIST;
	double xs = (x1-x0)/1000.0;
	double xm = 0.5*(x1+x0);	
	double em = field(EX, xm).sqabs()+field(EY, xm).sqabs()+field(EZ, xm).sqabs();
	
	for(double x=x0; x<=x1; x+=xs)
	{
		double e = field(EX, x).sqabs()+field(EY, x).sqabs()+field(EZ, x).sqabs();
		if(e > em) { xm = x; em = e; }
	}
	Complex ex = field(EX, xm);
	Complex ey = field(EY, xm);
	Complex f = CC1;
	if(ex.sqabs() > ey.sqabs()) f = ex.abs()/ex;
	else                        f = ey.abs()/ey;
	FAa1 *= f;
	FAa2 *= f;
	FAa3 *= f;
	FAa4 *= f;
	FAes *= f;
	FAhs *= f;
	FAec *= f;
	FAhc *= f;
	return;
}

/* integrate a product of field components cp1 (cc) and cp2
   over the interval [x0, x1]
   based on: Numerical Recipes in C --- The Art of Scientific Computing
   Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */
Complex TflnMode::integral(Fcomp cp1, Fcomp cp2, double x0, double x1) const
{
	static double x[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	Complex ovl = 0.0; 
	int numx = ((int) ceil((x1-x0)/wg.lambda*TFLN_SLAMS_NINTSPWL));
	double step = (x1-x0)/((double) numx);
	for(int i=0; i<=numx-1; ++i)
	{
		double a = x0+((double) i)*step;
		double b = a+step;
		double xm = 0.5*(b+a);
		double xr = 0.5*(b-a);
		Complex s = CC0;
		for(int j=1; j<=5; j++) 
		{
			double dx=xr*x[j];
			s += w[j]*( field(cp1, xm+dx).conj()*field(cp2, xm+dx)
				   +field(cp1, xm-dx).conj()*field(cp2, xm-dx));
		}
		ovl += s*xr;
	}
	return ovl;
}
/* ... over the domain of the mode */
Complex TflnMode::integral(Fcomp cp1, Fcomp cp2) const
{
	if(domain.x0 > wg.hx(0)) tflnerror("TflnMode::integral, domain / interface positions (0)");
	if(wg.hx(1) > domain.x1) tflnerror("TflnMode::integral, domain / interface positions (1)");
	double dx = wg.lambda/TFLN_SLAMS_NINTSPWL*1.0e-10;
	return integral(cp1, cp2, domain.x0+dx, wg.hx(0)-dx)
	     + integral(cp1, cp2, wg.hx(0)+dx,  wg.hx(1)-dx)
	     + integral(cp1, cp2, wg.hx(1)+dx, domain.x1-dx);
}


/* integrate a product of field components cp1 (cc) and cp2 related to modes m1 and m2  
   over the interval [x0, x1]
   based on: Numerical Recipes in C --- The Art of Scientific Computing
   Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */
Complex overlap(TflnMode m1, Fcomp cp1, TflnMode m2, Fcomp cp2, double x0, double x1)
{
	static double x[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	Complex ovl = 0.0; 
	int numx = ((int) ceil((x1-x0)/m1.wg.lambda*TFLN_SLAMS_NINTSPWL));
	double step = (x1-x0)/((double) numx);
	for(int i=0; i<=numx-1; ++i)
	{
		double a = x0+((double) i)*step;
		double b = a+step;
		double xm = 0.5*(b+a);
		double xr = 0.5*(b-a);
		Complex s = CC0;
		for(int j=1; j<=5; j++) 
		{
			double dx=xr*x[j];
			s += w[j]*( m1.field(cp1, xm+dx).conj()*m2.field(cp2, xm+dx)
				   +m1.field(cp1, xm-dx).conj()*m2.field(cp2, xm-dx));
		}
		ovl += s*xr;
	}
	return ovl;
}
/* ... over the intersection of the domains of both modes */
Complex overlap(TflnMode m1, Fcomp cp1, TflnMode m2, Fcomp cp2)
{
	double x0 = m1.domain.x0;
	if(m2.domain.x0 > x0) x0 = m2.domain.x0;
	double x1 = m1.domain.x1;
	if(m2.domain.x1 < x1) x1 = m2.domain.x1;
	if(x0 > x1) return CC0;
	Dvector ip(6);
	ip(0) = x0; 
	ip(1) = x1; 
	ip(2) = m1.wg.hx(0); 
	ip(3) = m1.wg.hx(1); 
	ip(2) = m2.wg.hx(0); 
	ip(3) = m2.wg.hx(1); 
	ip.sort(1);
	double dx = m1.wg.lambda/TFLN_SLAMS_NINTSPWL*1.0e-10;
	Complex ovl = CC0;
	for(int l=0; l<=4; ++l)
	{	
		double xa = ip(l)+dx;
		double xb = ip(l+1)-dx;
		if(xb>xa) ovl += overlap(m1, cp1, m2, cp2, xa, xb);
	}
	return ovl;
}


/* modal orthogonality: overlap of two TflnModes */
Complex overlap(TflnMode m1, TflnMode m2)
{
	return 0.25*(overlap(m1, EX, m2, HY) 
                   - overlap(m1, EY, m2, HX) 
                   + overlap(m1, HY, m2, EX) 
                   - overlap(m1, HX, m2, EY));
}

/* longitudinal component of the Poyntingvector, 
   integrated over the entire cross section */
/* ... the TE-associated part */
double TflnMode::tepower() const
{
	return -0.5*integral(EY, HX).re;
}
/* ... the TM-associated part */
double TflnMode::tmpower() const
{
	return  0.5*integral(EX, HY).re;
}
/* full power */
double TflnMode::power() const
{
	return tepower()+tmpower(); 
}

/* normalize mode to power() == 1 
   set polarizatio ratio & electric and magnetic field maxima */
void TflnMode::normalize()
{
	double tep = tepower();
	double tmp = tmpower();
	double f = 1.0/sqrt(tep+tmp);
	FAa1 *= f;
	FAa2 *= f;
	FAa3 *= f;
	FAa4 *= f;
	FAes *= f;
	FAhs *= f;
	FAec *= f;
	FAhc *= f;
	polratio = tep/(tep+tmp);
	maxE = 0.0;
	maxH = 0.0;
	double s = wg.lambda/TFLN_SLAMS_NINTSPWL;
	for(double x=domain.x0; x<=domain.x1; x+=s) 
	{
		double e = sqrt(field(EX, x).sqabs()+field(EY, x).sqabs()+field(EZ, x).sqabs());
		if(e > maxE) maxE = e;
		double h = sqrt(field(HX, x).sqabs()+field(HY, x).sqabs()+field(HZ, x).sqabs());
		if(h > maxH) maxH = h;
	}
	return;
}

/* field profile -> MATLAB .m file  
	cp: EX - HZ
	foa: MOD, ORG, SQR, REP, IMP 
	disp: output window on the x axis
	np: number of output points 
	ext0, ext1: filename id characters
	pltype: 'L': geometry information & plot commands (default)
	        'V': field, mesh, and plot comand, to be included into *L.m  */
void TflnMode::plot(Fcomp cp, Afo foa, Interval disp, int np, 
	  char ext0, char ext1, char pltype) const
{
	FILE *dat;
	double x, dx, t;
	double x0, x1;
	double xbd;
	double epsold, epsnew;
	double minf, maxf, f;
	int nsec, nbd;
	char name[13] = "tf______.m";
	name[2] = afochr(foa);
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	if(pltype != 'V') pltype = 'L';
	name[7] = pltype;

	Waveguide iwg =  wg.toisotropic();

	x0 = disp.x0;
	x1 = disp.x1;
	
	if(x0 > x1) { t = x0; x0 = x1; x1 = t; }

	switch(foa)
	{
	case ORG:
		fprintf(stderr, "   %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case MOD:
		fprintf(stderr, "   |%c%c| [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case SQR:
		fprintf(stderr, "   |%c%c|^2 [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case REP:
		fprintf(stderr, "   Re %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case IMP:
		fprintf(stderr, "   Im %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	}

	Cmatrix cfld(iwg.nx+2, np+1);
	Dmatrix pos(iwg.nx+2, np+1);
	Ivector nump(iwg.nx+2);
	Dvector of(np+1);
	Dvector op(np+1);
	Dvector bd(iwg.nx+1);
	Dvector ri(iwg.nx+2);
	nsec = 0;
	nbd  = 0;

	dx = x1-x0;
	nump(nsec) = 0;
	x = x0;
	epsold = iwg.eps(x);  	
	ri(nbd) = sgn_sqrt(epsold);
	pos(nsec, nump(nsec)) = x;  
	cfld(nsec, nump(nsec)) = field(cp, x);  
	++nump(nsec);
	for(int i=1; i<=np; ++i)
	{
		x = x0+((double) i)*dx/((double) np);
		epsnew = iwg.eps(x);
		if(fabs(epsnew-epsold)<1.0e-10)
		{
			pos(nsec, nump(nsec)) = x;  
			cfld(nsec, nump(nsec)) = field(cp, x);  
			++nump(nsec);
		}
		else
		{
			xbd = (wg.layer(x)).x0; 
			pos(nsec, nump(nsec)) = xbd;  
			cfld(nsec, nump(nsec)) = 	
				field(cp, xbd-dx*1.0e-8);  
			++nump(nsec);
			bd(nbd) = xbd; 
			++nbd;
			epsold = epsnew;
			ri(nbd) = sgn_sqrt(epsold);
			++nsec;
			nump(nsec) = 0;
			pos(nsec, nump(nsec)) = xbd;  
			cfld(nsec, nump(nsec)) = 	
				field(cp, xbd+dx*1.0e-8);  
			++nump(nsec);
			pos(nsec, nump(nsec)) = x;  
			cfld(nsec, nump(nsec)) = field(cp, x);  
			++nump(nsec);
		}
	}
	++nsec;

	dat = fopen(name, "w+");
	mlout_title(dat, name, "Mode profile");
	name[8] = 0;
	Dmatrix fld = realize(cfld, foa);
	minf = fld.min();
	maxf = fld.max();
	if(minf >= maxf) { minf -= 1.0;	maxf = minf+2.0; }

	if(pltype == 'L')
	{
		switch(foa)
		{
			case MOD: case SQR:
				minf = 0.0;
				break;
			case ORG: case REP: case IMP:
				f = fabs(minf); maxf = fabs(maxf);
				if(f > maxf) maxf = f;
				minf = -maxf;
				break;
		}
		mlout_geo(dat, iwg, minf, maxf); 
		mlout_Lsecgeo(dat, x0, x1, nbd, bd, ri);
	}

	for(int j=0; j<=nsec-1; ++j)
	{
		for(int i=0; i<=nump(j)-1; ++i)
		{
			of(i) = fld(j, i);
			op(i) = pos(j, i);
		}
		if(pltype == 'L')
			mlout_sec(dat, cp, dig10(j), dig1(j), ' ', ' ', 
				  nump(j), of, op); 
		else
			mlout_sec(dat, cp, dig10(j), dig1(j), ext0, ext1, 
				  nump(j), of, op); 
	}
	if(pltype == 'L') 
	{
		if(cp == EX || cp == EY || cp == EZ)
		{
			fprintf(dat, "maxf = %g;\n", maxE);
			fprintf(dat, "minf = %g;\n", -maxE);
		}
		if(cp == HX || cp == HY || cp == HZ)
		{
			fprintf(dat, "maxf = %g;\n", maxH);
			fprintf(dat, "minf = %g;\n", -maxH);
		}
		mlout_Lsecplot(name, dat, cp, foa, nbd, nsec); 
		mlout_print(dat, name, 'e');
	}
	else mlout_Vsecplot(dat, cp, nsec, ext0, ext1); 
	fclose(dat);
	return;
}

/* vectorial field profile -> MATLAB .m file  
	disp: output window on the x axis
	np: number of output points 
	ext0, ext1: filename id characters */
void TflnMode::plot(Interval disp, int np, char ext0, char ext1, char ext2, char ext3) const
{
	FILE *dat;
	double x, dx, t;
	double x0, x1;
	double xbd;
	double epsold, epsnew;
	int nsec, nbd;
	char name[13] = "vmp____.m";
	name[3] = ext0;
	name[4] = ext1;
	name[5] = ext2;
	name[6] = ext3;

	Waveguide iwg =  wg.toisotropic();

	x0 = disp.x0;
	x1 = disp.x1;
	
	if(x0 > x1) { t = x0; x0 = x1; x1 = t; }
	fprintf(stderr, "   EX--HZ [%g, %g] (%d) >> %s\n", x0, x1, np, name);

	dat = fopen(name, "w+");
	mlout_title(dat, name, "Vectorial mode profile");
	name[7] = 0;

	Cmatrix exfld(iwg.nx+2, np+1);
	Cmatrix eyfld(iwg.nx+2, np+1);
	Cmatrix ezfld(iwg.nx+2, np+1);
	Cmatrix hxfld(iwg.nx+2, np+1);
	Cmatrix hyfld(iwg.nx+2, np+1);
	Cmatrix hzfld(iwg.nx+2, np+1);
	Dmatrix pos(iwg.nx+2, np+1);
	Ivector nump(iwg.nx+2);
	Cvector of(np+1);
	Dvector op(np+1);
	Dvector bd(iwg.nx+1);
	Dvector ri(iwg.nx+2);
	nsec = 0;
	nbd  = 0;

	dx = x1-x0;
	nump(nsec) = 0;
	x = x0;
	epsold = iwg.eps(x);  	
	ri(nbd) = sgn_sqrt(epsold);
	pos(nsec, nump(nsec)) = x;  
	exfld(nsec, nump(nsec)) = field(EX, x);  
	eyfld(nsec, nump(nsec)) = field(EY, x);  
	ezfld(nsec, nump(nsec)) = field(EZ, x);  
	hxfld(nsec, nump(nsec)) = field(HX, x);  
	hyfld(nsec, nump(nsec)) = field(HY, x);  
	hzfld(nsec, nump(nsec)) = field(HZ, x);  
	++nump(nsec);
	for(int i=1; i<=np; ++i)
	{
		x = x0+((double) i)*dx/((double) np);
		epsnew = iwg.eps(x);
		if(fabs(epsnew-epsold)<1.0e-10)
		{
			pos(nsec, nump(nsec)) = x;  
			exfld(nsec, nump(nsec)) = field(EX, x);  
			eyfld(nsec, nump(nsec)) = field(EY, x);  
			ezfld(nsec, nump(nsec)) = field(EZ, x);  
			hxfld(nsec, nump(nsec)) = field(HX, x);  
			hyfld(nsec, nump(nsec)) = field(HY, x);  
			hzfld(nsec, nump(nsec)) = field(HZ, x);  
			++nump(nsec);
		}
		else
		{
			xbd = (wg.layer(x)).x0; 
			pos(nsec, nump(nsec)) = xbd;  
			exfld(nsec, nump(nsec)) = field(EX, xbd-dx*1.0e-8);
			eyfld(nsec, nump(nsec)) = field(EY, xbd-dx*1.0e-8);
			ezfld(nsec, nump(nsec)) = field(EZ, xbd-dx*1.0e-8);
			hxfld(nsec, nump(nsec)) = field(HX, xbd-dx*1.0e-8);
			hyfld(nsec, nump(nsec)) = field(HY, xbd-dx*1.0e-8);
			hzfld(nsec, nump(nsec)) = field(HZ, xbd-dx*1.0e-8);
			++nump(nsec);
			bd(nbd) = xbd; 
			++nbd;
			epsold = epsnew;
			ri(nbd) = sgn_sqrt(epsold);
			++nsec;
			nump(nsec) = 0;
			pos(nsec, nump(nsec)) = xbd;  
			exfld(nsec, nump(nsec)) = field(EX, xbd+dx*1.0e-8);
			eyfld(nsec, nump(nsec)) = field(EY, xbd+dx*1.0e-8);
			ezfld(nsec, nump(nsec)) = field(EZ, xbd+dx*1.0e-8);
			hxfld(nsec, nump(nsec)) = field(HX, xbd+dx*1.0e-8);
			hyfld(nsec, nump(nsec)) = field(HY, xbd+dx*1.0e-8);
			hzfld(nsec, nump(nsec)) = field(HZ, xbd+dx*1.0e-8);
			++nump(nsec);
			pos(nsec, nump(nsec)) = x;  
			exfld(nsec, nump(nsec)) = field(EX, x);  
			eyfld(nsec, nump(nsec)) = field(EY, x);  
			ezfld(nsec, nump(nsec)) = field(EZ, x);  
			hxfld(nsec, nump(nsec)) = field(HX, x);  
			hyfld(nsec, nump(nsec)) = field(HY, x);  
			hzfld(nsec, nump(nsec)) = field(HZ, x);  
			++nump(nsec);
		}
	}
	++nsec;

	double f;
	double minE, maxE, minH, maxH;
	maxE = 0.0;
	f = realize(exfld, MOD).max(); if(f > maxE) maxE = f;
	f = realize(eyfld, MOD).max(); if(f > maxE) maxE = f;
	f = realize(ezfld, MOD).max(); if(f > maxE) maxE = f;
	minE = -maxE;
	maxH = 0.0;
	f = realize(hxfld, MOD).max(); if(f > maxH) maxH = f;
	f = realize(hyfld, MOD).max(); if(f > maxH) maxH = f;
	f = realize(hzfld, MOD).max(); if(f > maxH) maxH = f;
	minH = -maxH;

	for(int j=0; j<=nsec-1; ++j)
	{
		for(int i=0; i<=nump(j)-1; ++i) op(i) = pos(j, i);
		for(int i=0; i<=nump(j)-1; ++i) of(i) = exfld(j, i); 
		mlout_cplxsec(dat, EX, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
		for(int i=0; i<=nump(j)-1; ++i) of(i) = eyfld(j, i); 
		mlout_cplxsec(dat, EY, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
		for(int i=0; i<=nump(j)-1; ++i) of(i) = ezfld(j, i); 
		mlout_cplxsec(dat, EZ, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
		for(int i=0; i<=nump(j)-1; ++i) of(i) = hxfld(j, i); 
		mlout_cplxsec(dat, HX, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
		for(int i=0; i<=nump(j)-1; ++i) of(i) = hyfld(j, i); 
		mlout_cplxsec(dat, HY, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
		for(int i=0; i<=nump(j)-1; ++i) of(i) = hzfld(j, i); 
		mlout_cplxsec(dat, HZ, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); 
	}

	mlout_geo(dat, iwg); 
	mlout_Lsecgeo(dat, x0, x1, nbd, bd, ri);

	fprintf(dat, "minE = %g;\n", minE);
	fprintf(dat, "maxE = %g;\n", maxE);
	fprintf(dat, "minH = %g;\n", minH);
	fprintf(dat, "maxH = %g;\n", maxH);

	mlout_vecsecplot(name, dat, nbd, nsec, neff); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}

/* - mode solver routines ---------------------------------- */

/* transverse resonance evaluation */
#define TFLN_Eigen true   // use Eigen subroutines, or not
double TFLN_TR_KROOT = 1.0e-10;
double TflnMode::travres(double bq, bool fin) 
{
	betaq = bq;
	beta = sqrt(bq);
	neff = beta/k0;
	kappas = sqrt(bq-k0*k0*wg.ns*wg.ns);
	kappac = sqrt(bq-k0*k0*wg.nc*wg.nc);

	double c1 = k0*k0*(wg.epsy-wg.delta*wg.delta/wg.epsz)-bq;
	double c2 = k0*k0*wg.epsz-wg.epsz/wg.epsx*bq;
	double b1 = k0*wg.delta/wg.epsz;
	double b2 = k0*wg.delta;

	Cmatrix Q(4, 4);
	Q.init(CC0);
	Q(0, 2) = CC1;
	Q(1, 3) = CC1;
	Q(2, 0) = c1;
	Q(3, 1) = c2;
	Q(2, 3) = -b1;
	Q(3, 2) = -b2;

	double b1b2pc1pc2 = b1*b2+c1+c2;
	double kroot = b1b2pc1pc2*b1b2pc1pc2-4.0*c1*c2;

/*	     if(fabs(kroot)/fabs(bbcc) < 1.0e-14) fprintf(stderr, "kroot < 1e-14\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-13) fprintf(stderr, "kroot < 1e-13\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-12) fprintf(stderr, "kroot < 1e-12\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-11) fprintf(stderr, "kroot < 1e-11\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-10) fprintf(stderr, "kroot < 1e-10\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-9) fprintf(stderr, "kroot < 1e-9\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-8) fprintf(stderr, "kroot < 1e-8\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-7) fprintf(stderr, "kroot < 1e-7\n");
	else if(fabs(kroot)/fabs(bbcc) < 1.0e-6) fprintf(stderr, "kroot < 1e-6\n");
*/
	if(fabs(kroot)/fabs(b1b2pc1pc2) < TFLN_TR_KROOT)
	{
		Complex kqr = sqrt_pos(0.5*b1b2pc1pc2);
		kappa1 =  kqr; CVe1 = CC1; CVh1 = CC0;
		kappa2 = -kqr; CVe2 = CC1; CVh2 = CC0;
		kappa3 =  kqr; CVe3 = CC0; CVh3 = CC1;
		kappa4 = -kqr; CVe4 = CC0; CVh4 = CC1;
	}
	else
	{
		Cvector kap(4);
		if(TFLN_Eigen)
		{
			int info;
			kap = Q.geneigenv(info);
			if(info != 0) tflnerror("travres: Q.geneigenv failed");
			Cvector qr;
			kappa1 = kap(0).tidy(); qr = Q.col(0); CVe1 = qr(0).tidy(); CVh1 = qr(1).tidy();
			kappa2 = kap(1).tidy(); qr = Q.col(1); CVe2 = qr(0).tidy(); CVh2 = qr(1).tidy();
			kappa3 = kap(2).tidy(); qr = Q.col(2); CVe3 = qr(0).tidy(); CVh3 = qr(1).tidy();
			kappa4 = kap(3).tidy(); qr = Q.col(3); CVe4 = qr(0).tidy(); CVh4 = qr(1).tidy();
		}
		else
		{
			double kqp = 0.5*(b1b2pc1pc2 + sqrt(kroot));
			double kqm = 0.5*(b1b2pc1pc2 - sqrt(kroot));
			kap(0) = sqrt_pos(kqp).tidy();
			kap(1) = sqrt_neg(kqp).tidy();
			kap(2) = sqrt_pos(kqm).tidy();
			kap(3) = sqrt_neg(kqm).tidy();

			for(int evi=0; evi<=3; ++evi)
			{
				Dmatrix U; 
				U.unit(4);
				Cmatrix M = sub(Q, mult(kap(evi), U));
				Dmatrix A(8, 8);
				A.setsubmatrix(     M.re(),        0, 0); 
				A.setsubmatrix(mult(M.im(), -1.0), 0, 4); 
				A.setsubmatrix(     M.im(),        4, 0); 
				A.setsubmatrix(     M.re(),        4, 4); 
				Dvector w;
				Dmatrix V;
				A.svd(w, V, true);
				int i;
				for(i=0; i<=7; ++i) w(i) = fabs(w(i));
				w.min(i);
				Dvector s = V.col(i);
				Cvector qr;
				qr = add(s.subvector(4, 0), mult(s.subvector(4, 4), CCI));
				     if(evi == 0) { kappa1 = kap(evi); CVe1 = qr(0).tidy(); CVh1 = qr(1).tidy(); }
				else if(evi == 1) { kappa2 = kap(evi); CVe2 = qr(0).tidy(); CVh2 = qr(1).tidy(); }
				else if(evi == 2) { kappa3 = kap(evi); CVe3 = qr(0).tidy(); CVh3 = qr(1).tidy(); }
				else              { kappa4 = kap(evi); CVe4 = qr(0).tidy(); CVh4 = qr(1).tidy(); }
			}
		}
	}

	Cmatrix M(8, 8);
	M.init(CC0);
	double d = wg.hx(1)-wg.hx(0);

	M(0, 4) = CC1; 
	M(0, 1-1) = -CVe1; 
	M(0, 2-1) = -CVe2; 
	M(0, 3-1) = -CVe3; 
	M(0, 4-1) = -CVe4;

	M(1, 4) = kappas; 
	M(1, 1-1) = CCI*kappa1*CVe1;
	M(1, 2-1) = CCI*kappa2*CVe2;
	M(1, 3-1) = CCI*kappa3*CVe3;
	M(1, 4-1) = CCI*kappa4*CVe4;

	M(2, 5) = CC1; 
	M(2, 1-1) = -CVh1; 
	M(2, 2-1) = -CVh2; 
	M(2, 3-1) = -CVh3; 
	M(2, 4-1) = -CVh4;

	M(3, 5) = kappas/(wg.ns*wg.ns); 
	M(3, 1-1) = CCI/wg.epsz*kappa1*CVh1 + CCI*k0*wg.delta/wg.epsz*CVe1;
	M(3, 2-1) = CCI/wg.epsz*kappa2*CVh2 + CCI*k0*wg.delta/wg.epsz*CVe2;
	M(3, 3-1) = CCI/wg.epsz*kappa3*CVh3 + CCI*k0*wg.delta/wg.epsz*CVe3;
	M(3, 4-1) = CCI/wg.epsz*kappa4*CVh4 + CCI*k0*wg.delta/wg.epsz*CVe4;

	M(4, 6) = CC1; 
	M(4, 1-1) = -CVe1*expmi(kappa1*d); 
	M(4, 2-1) = -CVe2*expmi(kappa2*d); 
	M(4, 3-1) = -CVe3*expmi(kappa3*d); 
	M(4, 4-1) = -CVe4*expmi(kappa4*d);

	M(5, 6) = -kappac; 
	M(5, 1-1) = CCI*kappa1*CVe1*expmi(kappa1*d);
	M(5, 2-1) = CCI*kappa2*CVe2*expmi(kappa2*d);
	M(5, 3-1) = CCI*kappa3*CVe3*expmi(kappa3*d);
	M(5, 4-1) = CCI*kappa4*CVe4*expmi(kappa4*d);

	M(6, 7) = CC1; 
	M(6, 1-1) = -CVh1*expmi(kappa1*d); 
	M(6, 2-1) = -CVh2*expmi(kappa2*d); 
	M(6, 3-1) = -CVh3*expmi(kappa3*d); 
	M(6, 4-1) = -CVh4*expmi(kappa4*d);

	M(7, 7) = -kappac/(wg.nc*wg.nc); 
	M(7, 1-1) = (CCI/wg.epsz*kappa1*CVh1 + CCI*k0*wg.delta/wg.epsz*CVe1)*expmi(kappa1*d);
	M(7, 2-1) = (CCI/wg.epsz*kappa2*CVh2 + CCI*k0*wg.delta/wg.epsz*CVe2)*expmi(kappa2*d);
	M(7, 3-1) = (CCI/wg.epsz*kappa3*CVh3 + CCI*k0*wg.delta/wg.epsz*CVe3)*expmi(kappa3*d);
	M(7, 4-1) = (CCI/wg.epsz*kappa4*CVh4 + CCI*k0*wg.delta/wg.epsz*CVe4)*expmi(kappa4*d);

	Dvector w;
	Dmatrix V;
	Cvector mv;
	double tr;
	if(TFLN_Eigen)
	{
		int info;
		mv = M.geneigenv(info);
		if(info != 0) tflnerror("travres: M.geneigenv failed");
		tr = (mv(0)*mv(1)*mv(2)*mv(3)*mv(4)*mv(5)*mv(6)*mv(7)).abs();
	}
	else
	{
		Dmatrix A(16, 16);
		A.setsubmatrix(     M.re(),        0, 0); 
		A.setsubmatrix(mult(M.im(), -1.0), 0, 8); 
		A.setsubmatrix(     M.im(),        8, 0); 
		A.setsubmatrix(     M.re(),        8, 8); 
		A.svd(w, V, true);
		tr = 1.0; 
		for(int i=0; i<=15; ++i) tr *= w(i);
		tr = sqrt(tr);
	}

	if(fin == false)
	{
	//		fprintf(stderr, "TR: bq: %g, b: %g, N: %g -> %g\n", bq, beta, neff, tr);
		return tr;
	}

	Cvector av;
	if(TFLN_Eigen)
	{
		int svi = 0;
		double sva = mv(svi).abs();
		for(int j=1; j<=7; ++j)
		{
			if(mv(j).abs() < sva)
			{
				svi = j;
				sva = mv(svi).abs();
			}
		}
		av = M.col(svi);
	}
	else
	{
		int svi;
		w.min(svi);
		Dvector s = V.col(svi);
		av = add(s.subvector(8, 0), mult(s.subvector(8, 8), CCI));
	}

	FAa1 = av(0);
	FAa2 = av(1);
	FAa3 = av(2);
	FAa4 = av(3);
	FAes = av(4);
	FAhs = av(5);
	FAec = av(6);
	FAhc = av(7);
	return tr;
}

/* final mode: TR evaluation, residuum */
Complex residuum;

/* check mode: interface conditions */
void TflnMode::check()
{
	double x;
	double dx = 1.0e-13;
	double epsxp, epsxm;

	x = wg.hx(0);
	epsxp = wg.epsx;
	epsxm = wg.ns*wg.ns;
	fprintf(stderr, "    x = %g, (eE)x: %g\n", x, (epsxp*field(EX, x+dx)-epsxm*field(EX, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Ey: %g\n", x, (field(EY, x+dx)-field(EY, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Ez: %g\n", x, (field(EZ, x+dx)-field(EZ, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Hx: %g\n", x, (field(HX, x+dx)-field(HX, x-dx)).abs()/maxH);
	fprintf(stderr, "    x = %g,    Hy: %g\n", x, (field(HY, x+dx)-field(HY, x-dx)).abs()/maxH);
	fprintf(stderr, "    x = %g,    Hz: %g\n", x, (field(HZ, x+dx)-field(HZ, x-dx)).abs()/maxH);
	x = wg.hx(1);
	epsxp = wg.nc*wg.nc;
	epsxm = wg.epsx;
	fprintf(stderr, "    x = %g, (eE)x: %g\n", x, (epsxp*field(EX, x+dx)-epsxm*field(EX, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Ey: %g\n", x, (field(EY, x+dx)-field(EY, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Ez: %g\n", x, (field(EZ, x+dx)-field(EZ, x-dx)).abs()/maxE);
	fprintf(stderr, "    x = %g,    Hx: %g\n", x, (field(HX, x+dx)-field(HX, x-dx)).abs()/maxH);
	fprintf(stderr, "    x = %g,    Hy: %g\n", x, (field(HY, x+dx)-field(HY, x-dx)).abs()/maxH);
	fprintf(stderr, "    x = %g,    Hz: %g\n", x, (field(HZ, x+dx)-field(HZ, x-dx)).abs()/maxH);

	dx = 1.0e-6;
	double rE1=0.0, rE2=0.0, rE3=0.0, rH1=0.0, rH2=0.0, rH3=0.0;
	for(x = wg.hx(0)-1.0; x <= wg.hx(1)+1.0; x+=(wg.hx(1)-wg.hx(0)+2.0)/20.0)
	{
		if((x+dx < wg.hx(0)) || (x-dx > wg.hx(1)) || ((wg.hx(0)<x-dx) && (x+dx<wg.hx(1))))
		{
			Complex ex = field(EX, x);
			Complex ey = field(EY, x);
			Complex ez = field(EZ, x);
			Complex hx = field(HX, x);
			Complex hy = field(HY, x);
			Complex hz = field(HZ, x);
			Complex dey = (field(EY, x+dx)-field(EY, x-dx))/(2.0*dx);
			Complex dez = (field(EZ, x+dx)-field(EZ, x-dx))/(2.0*dx);
			Complex dhy = (field(HY, x+dx)-field(HY, x-dx))/(2.0*dx);
			Complex dhz = (field(HZ, x+dx)-field(HZ, x-dx))/(2.0*dx);
			
			double delta, epsx, epsy, epsz;
			if(x < wg.hx(0))      { epsx = wg.ns*wg.ns; epsy = wg.ns*wg.ns; epsz = wg.ns*wg.ns; delta = 0; }
			else if(x > wg.hx(1)) { epsx = wg.nc*wg.nc; epsy = wg.nc*wg.nc; epsz = wg.nc*wg.nc; delta = 0; }
			     else             { epsx = wg.epsx;     epsy = wg.epsy;     epsz = wg.epsz;     delta = wg.delta; }
		
			Complex lhs, rhs;
			double err;

			lhs =  CCI*beta*ey;      rhs = CCmI/invommu0*hx; err = (lhs-rhs).abs()/maxE; if(err > rE1) rE1 = err;
			lhs =  CCmI*beta*ex-dez; rhs = CCmI/invommu0*hy; err = (lhs-rhs).abs()/maxE; if(err > rE2) rE2 = err;
			lhs =  dey;              rhs = CCmI/invommu0*hz; err = (lhs-rhs).abs()/maxE; if(err > rE3) rE3 = err;

			lhs =  CCI*beta*hy;      rhs = CCI/invomep0*epsx*ex;            err = (lhs-rhs).abs()/maxH; if(err > rH1) rH1 = err;
			lhs =  CCmI*beta*hx-dhz; rhs = CCI/invomep0*(epsy*ey+delta*ez); err = (lhs-rhs).abs()/maxH; if(err > rH2) rH2 = err;
			lhs =  dhy;              rhs = CCI/invomep0*(delta*ey+epsz*ez); err = (lhs-rhs).abs()/maxH; if(err > rH3) rH3 = err;
		}
	}
	fprintf(stderr, "    rE1: %g\n", rE1);
	fprintf(stderr, "    rE2: %g\n", rE2);
	fprintf(stderr, "    rE3: %g\n", rE3);
	fprintf(stderr, "    rH1: %g\n", rH1);
	fprintf(stderr, "    rH2: %g\n", rH2);
	fprintf(stderr, "    rH3: %g\n", rE3);
	return;
}

/* --------------------------------------------------------------- 
an array of TFLN modes */

/* initialize */
TflnModeArray::TflnModeArray()
{
	num = 0;
	mvec = NULL;
}
/* destroy */
TflnModeArray::~TflnModeArray()
{
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0;
}
/* copy constructor */
TflnModeArray::TflnModeArray(const TflnModeArray& ma)
{	
	num = ma.num;
	mvec = new TflnMode[num];
	TflnMode* ap = ma.mvec;
	TflnMode* mp =    mvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}

/* assignment */
TflnModeArray& TflnModeArray::operator=(const TflnModeArray& ma)
{
	if(this != &ma)
	{
		if(mvec != NULL) delete[] mvec;
		num = ma.num;
		mvec = new TflnMode[num];
		TflnMode *ap = ma.mvec;
		TflnMode *mp =    mvec;
		for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
	}
	return *this;
}

/* free allocated memory */
void TflnModeArray::clear()
{
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0;
}

/* add a mode */
void TflnModeArray::add(TflnMode m)
{
	TflnMode *avec;
	avec = new TflnMode[num+1];
	TflnMode* ap = avec;
	TflnMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	*ap = m;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	++num;
	return;
}
/* delete a mode entry */
void TflnModeArray::remove(int i)
{
	if(i<0 || i>=num) tflnerror("TflnModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	TflnMode *avec;
	avec = new TflnMode[num-1];
	TflnMode* ap = avec;
	TflnMode* mp = mvec;
	for(int j=0; j<=i-1; ++j) *ap++ = *mp++;
	mp++;
	for(int j=i+1; j<=num-1; ++j) *ap++ = *mp++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	--num;
	return;
}
/* add an entire TflnModeArray nma, nma is cleared ! */
void TflnModeArray::merge(TflnModeArray& ma)
{
	TflnMode *avec;
	avec = new TflnMode[num+ma.num];
	TflnMode* ap = avec;
	TflnMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	TflnMode* np = ma.mvec;
	for(int i=0; i<=ma.num-1; ++i) *ap++ = *np++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	num += ma.num;
	return;
}
/* sort the array by squared propagation constants, highest first */
void TflnModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	TflnMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		maxi = j; 
		maxb = mvec[j].betaq;
		for(k=j+1; k<=num-1; ++k)
		{
			if(maxb<mvec[k].betaq-(1.0e-13))
			{
				maxb = mvec[k].betaq;
				maxi = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t;
	}
	return;
}

/* --------------------------------------------------------------- 
mode analysis of the TFLN waveguide */

/* computational parameters */
double TFLN_SLAMS_BQSEP = 1.0e-8;
int TFLN_SLAMS_STEPS = 500;
double TFLN_SLAMS_BQTOL = 1.0e-13;
double TFLN_SLAMS_TRERROR = 1.0e-9;
double TFLN_SLAMS_MINKAPPA = 1.0e-5;
double TFLN_SLAMS_DOMLVL = 1.0e-8;
double TFLN_SLAMS_NINTSPWL = 100.0;

#define GSR 0.38197
int modeanalysis(TflnWaveguide wg, TflnModeArray& ma, int quiet)
{
	double nmin = wg.ns; if(wg.nc > nmin) nmin = wg.nc;
	double nmax = wg.no; if(wg.ne > nmax) nmax = wg.ne;
	TflnMode m(wg);
	double k0 = m.k0;

if(quiet != 1)
{
fprintf(stderr, "\n------------- TFLN - MA ------------------------------------------- '23 ---\n");
fprintf(stderr, "n_eff: (%.10g, %.10g)  ", nmin, nmax);
fprintf(stderr, "Tol: %g  ", TFLN_SLAMS_BQTOL);
fprintf(stderr, "Sep: %g  ", TFLN_SLAMS_BQSEP);
fprintf(stderr, "#steps: %d\n", TFLN_SLAMS_STEPS);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "hx: ");
for(int j=0; j<=1; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "ns: %g  no: %g  ne: %g  nc: %g\n", wg.ns, wg.no, wg.ne, wg.nc);
if(wg.type == Xcut) fprintf(stderr, "Xcut, theta: %g\n", wg.theta);
else                fprintf(stderr, "Zcut (theta: %g)\n", wg.theta);
fprintf(stderr, "epsx: %g  epsy: %g  epsz: %g  delta: %g\n", wg.epsx, wg.epsy, wg.epsz, wg.delta);
fprintf(stderr, "lambda: %.10g  k0: %g\n", wg.lambda, k0);
fprintf(stderr, "----------------------------------------------------------------------------\n");
}
	ma.clear();

/*
for(double n = nmin; n <= nmax; n += (nmax-nmin)/100000)
{
	double bq = k0*k0*n*n;
	double tr = m.travres(bq, false);
	apptoxyf("tr", n, tr);
}
exit(1);
*/

	double bqmin = nmin*nmin*k0*k0+TFLN_SLAMS_BQSEP;
	double bqmax = nmax*nmax*k0*k0-TFLN_SLAMS_BQSEP;

	double bqs = (bqmax - bqmin)/((double) TFLN_SLAMS_STEPS);
	Dvector trv0(TFLN_SLAMS_STEPS+1);
	Dvector bqv0(TFLN_SLAMS_STEPS+1);
	for(int s=0; s<=TFLN_SLAMS_STEPS; ++s)
	{
		double bq = bqmin+((double) s)*bqs;
		bqv0(s) = bq;
		trv0(s) = m.travres(bq, false);
	}

	Dvector trv; trv.strip();
	Dvector bqv; bqv.strip();
	int ls = 0;
	for(int s=0; s<=bqv0.nel-3; ++s)
	{
		if(trv0(s) > trv0(s+1) && trv0(s+1) < trv0(s+2))
		{
			int a0 = ls;
			int a1 = s;
			trv.append(trv0.subvector(a1-a0+1, a0));
			bqv.append(bqv0.subvector(a1-a0+1, a0));
			double tbqmin = bqv0(s);
			double tbqmax = bqv0(s+2);
			double tbqs = (tbqmax - tbqmin)/((double) TFLN_SLAMS_STEPS);
			Dvector tbqv(TFLN_SLAMS_STEPS-1);
			Dvector ttrv(TFLN_SLAMS_STEPS-1);
			for(int k=1; k<=TFLN_SLAMS_STEPS-1; ++k)
			{
				double bq = tbqmin+((double) k)*tbqs;
				tbqv(k-1) = bq;
				ttrv(k-1) = m.travres(bq, false);
			}
			trv.append(ttrv);
			bqv.append(tbqv);
			ls = s+2;
			++s;
		}
	}
// toxyf("tr", bqv, trv); exit(1);
/*
trv = trv0;
bqv = bqv0;
*/
// toxyf("tr0", bqv, trv); exit(1);
	for(int s=0; s<=bqv.nel-3; ++s)
	{
		if(trv(s) > trv(s+1) && trv(s+1) < trv(s+2))
		{
			double bq0 = bqv(s); 
			double bq1 = bqv(s+1); 
			double bq2 = bqv(s+2);
			double tr1 = trv(s+1); 
			while(fabs(bq2-bq0) > TFLN_SLAMS_BQTOL)
			{
				double bq, tr;
				if(bq1-bq0 >= bq2-bq1)
				{
					bq = bq1-GSR*(bq1-bq0);
					tr = m.travres(bq, false);
					if(tr > tr1) { bq0 = bq; }
					else         { bq2 = bq1; bq1 = bq; tr1 = tr; }
				}
				else 	
				{
					bq = bq1+GSR*(bq2-bq1);
					tr = m.travres(bq, false);
					if(tr > tr1) { bq2 = bq; }
					else         { bq0 = bq1; bq1 = bq; tr1 = tr; }
				}
			}
			m.residuum = m.travres(0.5*(bq0+bq2), true);
			double minkappa = 1.0e+200;
			if(m.kappa1.abs() < minkappa) minkappa = m.kappa1.abs();
			if(m.kappa2.abs() < minkappa) minkappa = m.kappa2.abs();
			if(m.kappa3.abs() < minkappa) minkappa = m.kappa3.abs();
			if(m.kappa4.abs() < minkappa) minkappa = m.kappa4.abs();
			if(quiet != 1)
			{
fprintf(stderr, "neff = %g\n", m.neff);
/*
fprintf(stderr, "    kappac = %g\n", m.kappac);
fprintf(stderr, "    kappa1 = %g+i%g, |.| = %g\n", m.kappa1.re, m.kappa1.im, m.kappa1.abs());
fprintf(stderr, "    kappa2 = %g+i%g, |.| = %g\n", m.kappa2.re, m.kappa2.im, m.kappa2.abs());
fprintf(stderr, "    kappa3 = %g+i%g, |.| = %g\n", m.kappa3.re, m.kappa3.im, m.kappa3.abs());
fprintf(stderr, "    kappa4 = %g+i%g, |.| = %g\n", m.kappa4.re, m.kappa4.im, m.kappa4.abs());
fprintf(stderr, "    kappas = %g\n", m.kappas);
*/
			}
			if(m.residuum.abs() < TFLN_SLAMS_TRERROR && minkappa > TFLN_SLAMS_MINKAPPA)
			{
				m.domain = Interval(m.wg.hx(0)-(-log(TFLN_SLAMS_DOMLVL)/m.kappas), 
						    m.wg.hx(1)+(-log(TFLN_SLAMS_DOMLVL)/m.kappac)); 
				m.normalize();
				m.adjustphase();
				ma.add(m);
				if(quiet != 1)
				{
fprintf(stderr, "    domain: [%g, %g]\n", m.domain.x0, m.domain.x1);
fprintf(stderr, "    polratio: %g", m.polratio);
     if(m.polratio < 0.25) fprintf(stderr, " - TM\n");
else if(m.polratio > 0.75) fprintf(stderr, " - TE\n");
else                       fprintf(stderr, " - XX\n");
//					m.check();
fprintf(stderr, "    TR, residuum = %g+i%g, |.| = %g --- OK.\n", m.residuum.re, m.residuum.im, m.residuum.abs());
				}
			}
			else
			{
				if(quiet != 1)
				{
fprintf(stderr, "    TR, residuum = %g+i%g, |.| = %g --- discarded.\n", m.residuum.re, m.residuum.im, m.residuum.abs());
				}
			}
		}	
	}
	ma.sort();
	if(quiet != 1)
	{
fprintf(stderr, "---------------------------------------------------------------------------\n\n");
	}
	return ma.num;
}