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

/* 
 * cplxwg.cpp
 * Modes of waveguides with complex permittivity
 */

#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"slamode.h"
#include"integral.h"
#include"slamarr.h"
#include"slams.h"
#include"bepfld.h"
#include"cplxwg.h"

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

/* ------------------------------------------------------------------- */
// -> matlvis.h.cpp

void mlout_geo(FILE *dat, cWaveguide wg, double minf, double maxf)
{
        int i;
        fprintf(dat, "%% cWaveguide geometry\n");
        fprintf(dat, "nlX = %d;\n", wg.nx);
        fprintf(dat, "bdX = [");
        for(i=0; i<=wg.nx; ++i) fprintf(dat, "%g ", wg.hx(i));
        fprintf(dat, "];\n");
        fprintf(dat, "\n");
        fprintf(dat, "%% Refractive indices\n");
        fprintf(dat, "n_re = [");
       	for(i=0; i<=wg.nx+1; ++i) fprintf(dat, "%g ", wg.n(i).re);
        fprintf(dat, "];\n");
        fprintf(dat, "n_im = [");
       	for(i=0; i<=wg.nx+1; ++i) fprintf(dat, "%g ", wg.n(i).im);
        fprintf(dat, "];\n");
	fprintf(dat, "n = sqrt(n_re.*n_re+n_im.*n_im);\n");
        fprintf(dat, "\n");
        fprintf(dat, "%% Bounds on field values\n");
        fprintf(dat, "minf = %g;\n", minf);
        fprintf(dat, "maxf = %g;\n", maxf);
        fprintf(dat, "\n");
        return;
}

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

/* waveguide geometry + permittivity profile; vacuum wavelength ...,
   complex media, including losses or gain */

/* initialize */
cWaveguide::cWaveguide()
{
	nx = 0;
	hx = Dvector();
	lambda = 1.0;
	n = Cvector();
}
cWaveguide::cWaveguide(int vnx) 
{
	nx = vnx;
	hx = Dvector(nx+1);
	lambda = 1.0;
	n = Cvector(nx+2);
}
cWaveguide::cWaveguide(int vnx, Dvector vhx, double l, Cvector vn)
{
	nx = vnx;
	hx = vhx;
	lambda = l;
	n = vn;
}
cWaveguide::cWaveguide(Waveguide wg)
{
	nx = wg.nx;
	hx = wg.hx;
	lambda = wg.lambda;
	n = Cvector(wg.n);
	if(wg.special)
	{
		double e;	
		for(int j=0; j<=n.nel-1; ++j)
		{
			e = wg.n(j);
			if(e >= 0.0) n(j) = sqrt(e);
			else n(j) = CCmI*sqrt(fabs(e));
		}
	}
}

/* free allocated memory */
void cWaveguide::strip()
{
	nx = 0;
	hx.strip();
	n.strip();
}

/* get Interval index corresponding to position x */ 
int cWaveguide::layeridx(double x) const
{
	int l=0;
	while(l<=nx && x>hx(l)) ++l;
	return l;
}

/* get layer boundaries corresponding to index l */  
Interval cWaveguide::layer(int l) const
{
	Interval i;
	if(l <= 0) i.x0 = -AWAY;
	else i.x0 = hx(l-1);
	if(l >= nx+1) i.x1 = AWAY;
	else i.x1 = hx(l);
	return i;
}

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

/* permittivity on layer l */
Complex cWaveguide::eps(int l) const
{
	Complex e;
	e = n(l);
	return e*e;
}

/* permittivity at position x */
Complex cWaveguide::eps(double x) const
{
	Complex e;
	int l;
	l = layeridx(x);
	e = n(l);
	return e*e;
}

/* translate: hx -> hx+dx */
void cWaveguide::translate(double dx)
{
	int l;
	for(l=0; l<=nx; ++l) hx(l) += dx;
	return;
}

/* output to FILE dat */
void cWaveguide::write(FILE *dat) const
{
	int i;
	if(dat == stderr || dat == stdout)
	{
		comment("cWaveguide", dat);
		fprintf(dat, "Nx: %d, Hx: ", nx);
		for(i=0; i<=nx; ++i) fprintf(dat, "%g ", hx(i));
		fprintf(dat, "\n");
		fprintf(dat, "lambda: %g\n", lambda);
		fprintf(dat, "n:\n");
		for(i=0; i<=nx+1; ++i) { fprintf(dat, "  "); n(i).nice(dat); }
		fprintf(dat, "\n");
		fprintf(dat, "eps:\n");
		for(i=0; i<=nx+1; ++i) { fprintf(dat, "  "); eps(i).nice(dat); }
		fprintf(dat, "\n");
	}
	else
	{
		comment("cWaveguide", dat);
		comment("nx", dat);
		fputint(nx, dat);
		comment("hx[0..nx]", dat);
		hx.write(dat);
		comment("lambda", dat);
		fputdouble(lambda, dat);
		comment("n[0..nx+1]", dat);
		n.write(dat);
	}
	return;
}

/* input from FILE dat */
void cWaveguide::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		cplxwgerror("cWaveguide: read");
	}
	nx = fgetint(dat);
	hx.read(dat);
	lambda = fgetdouble(dat);
	n.read(dat);
	return;
}


/* compare the present object to wg, return 1 if equal, otherwise return 0 */
int cWaveguide::equal(const cWaveguide& wg) const 
{
	if(this == &wg) return 1;
	if(nx != wg.nx) return 0;
	if(fabs(lambda-wg.lambda) > COMPTOL_LAMBDA) return 0; 
	for(int l=0; l<=nx+1; ++l)
		if((n(l)-wg.n(l)).abs() > COMPTOL_N) return 0; 
	for(int l=0; l<=nx; ++l)
		if(fabs(hx(l)-wg.hx(l)) > COMPTOL_HX) return 0; 
	return 1;
}

/* remove unnecessary boundaries from waveguide */
void cWaveguide::optimize()
{
	int done = 0;
	int rl = 0;
	while(nx >= 2 && done != 1)	
	{	
		done = 1;
		for(int l=0; l<=nx; ++l)
		{
			if(done == 1 && (n(l)-n(l+1)).abs() < COMPTOL_N)
			{
				done = 0;
				rl = l; 
			}
		}
		if(done != 1)
		{
			for(int l=rl; l<=nx; ++l) n(l) = n(l+1);
			for(int l=rl; l<=nx-1; ++l) hx(l) = hx(l+1);
			--nx;
		}
	}
	Cvector n_n = n.subvector(nx+2, 0);
	Dvector n_hx = hx.subvector(nx+1, 0);
	n = n_n;
	hx = n_hx;
	return;
}

/* check consistency of the waveguide object: error&exit, if not ok */
void cWaveguide::consistency() const
{
	if(nx < 0 || nx >= 1000) cplxwgerror("cWaveguide: consistency: nx");
	if(hx.nel != nx+1) cplxwgerror("cWaveguide: consistency, hx: #elements");
	for(int b=0; b<=nx-1; ++b)	
	{
		if(hx(b) > hx(b+1) - HDIST)
		cplxwgerror("cWaveguide: consistency, hx: order/distance");
	}	
	if(lambda < 0.001 && lambda > 100.0)
		cplxwgerror("cWaveguide: consistency, lambda: magnitude");
	if(n.nel != nx+2) cplxwgerror("cWaveguide: consistency, n: #elements");
	return;
}

/* check the symmetry of the waveguide, return value l
   l == 0: no symmetry
   l >= 1: symmetric structure with respect to layer l */
int cWaveguide::checksymmetry() const
{
	if(nx % 2 == 0) return 0;
	double c = (hx(nx)+hx(0))/2.0;
	int cl = layeridx(c);
	for(int l=0; l<=nx; ++l)
	{
		if(fabs((c-hx(l)) - (hx(nx-l)-c))> COMPTOL_HX) return 0;
	}
	for(int l=0; l<=nx+1; ++l)
	{
		if((n(l)-n(nx-l+1)).abs() > COMPTOL_N) return 0;
	}
	return cl;
}

/* split the present waveguide into two, at layer l
   the pieces are stored in wg1 and wg2 */
void cWaveguide::split(int sl, cWaveguide& wg1, cWaveguide& wg2) const
{
	if(sl == 1)
	{
		wg1 = cWaveguide(1);
		wg1.hx(0) = hx(0)-LIMBDSEP/2.0;
		wg1.hx(1) = hx(0);
		wg1.n(0) = n(0);
		wg1.n(1) = n(0);
		wg1.n(2) = n(1);
		wg1.lambda = lambda;
	}
	else
	{
		int nnx = sl-1;
		wg1 = cWaveguide(nnx);
		for(int l=0; l<=nnx; ++l) wg1.hx(l) = hx(l);
		for(int l=0; l<=nnx+1; ++l) wg1.n(l) = n(l);
		wg1.lambda = lambda;
	}
	if(sl == nx)
	{
		wg2 = cWaveguide(1);
		wg2.hx(0) = hx(nx);
		wg2.hx(1) = hx(nx)+LIMBDSEP/2.0;
		wg2.n(0) = n(nx);
		wg2.n(1) = n(nx+1);
		wg2.n(2) = n(nx+1);
		wg2.lambda = lambda;
	}
	else
	{
		int nnx = nx-sl;
		wg2 = cWaveguide(nnx);
		for(int l=0; l<=nnx; ++l) wg2.hx(l) = hx(sl+l);
		for(int l=0; l<=nnx+1; ++l) wg2.n(l) = n(sl+l);
		wg2.lambda = lambda;
	}
	return;
}

/* expand the waveguide towards a structure 
   that is symmetric with respect to xs; xs > hx(nx) */
cWaveguide cWaveguide::expand(double xs)
{
	if(xs < hx(nx)+HDIST) cplxwgerror("cWaveguide: expand: xs");
	int nnx = 2*nx+1;
	cWaveguide wg(nnx);
	for(int b=0; b<=nx; ++b) 
	{
		wg.hx(b) = hx(b);
		wg.hx(nnx-b) = xs+(xs-hx(b));
	}
	for(int l=0; l<=nx; ++l)
	{
		 wg.n(l) = n(l);
		 wg.n(nnx+1-l) = n(l);
	}
	wg.n(nx+1) = n(nx+1);
	wg.lambda = lambda;
	return wg;
}

/* mirror the waveguide at xs; xs > hx(nx) */
cWaveguide cWaveguide::mirror(double xs)
{
	if(xs < hx(nx)+HDIST) cplxwgerror("cWaveguide: mirror: xs");
	cWaveguide wg(nx);
	for(int b=0; b<=nx; ++b) 
	{
		wg.hx(nx-b) = xs+(xs-hx(b));
	}
	for(int l=0; l<=nx+1; ++l)
	{
		 wg.n(nx+1-l) = n(l);
	}
	wg.lambda = lambda;
	return wg;
}

/* helper function: match interface positions of two waveguides */
Dvector itfmatch(const Dvector a, const Dvector b)
{
	Dvector c(a.nel+b.nel);
	c.setsubvector(a, 0);
	c.setsubvector(b, a.nel);
	c.sort();
	int done = 0;
	int ri, i; 
	while(c.nel > 1 && done == 0)
	{
		ri = -1;
		i = 0;
		while(i<=c.nel-2 && ri < 0)
		{
			if(fabs(c(i)-c(i+1)) < HDIST) ri = i;
			++i;
		}
		if(ri >= 0) c.remove(ri);
		else done = 1;
	}
	return c;
}


/* distance in permittivity between this and twg, 0: equal permittivities */
double cWaveguide::distance(const cWaveguide twg) const
{
	double d = 0.0;
	Dvector thx = itfmatch(hx, twg.hx);
	for(int l=0; l <=thx.nel; ++l) 
	{
		double x;
		if(l == 0)            x = thx(0)-1.0;
		else if(l == thx.nel) x = thx(thx.nel-1)+1.0;
		else                  x = 0.5*(thx(l-1)+thx(l));
		d += (eps(x)-twg.eps(x)).sqabs();
	}
	return sqrt(d);
}

/* an intermediate between this and twg, at ratio r */
cWaveguide cWaveguide::intermediate(const cWaveguide twg, double r) const
{
	Dvector ihx = itfmatch(hx,twg.hx);
	cWaveguide iw(ihx.nel-1);
	iw.hx = ihx;
	for(int l=0; l <=ihx.nel; ++l) 
	{
		double x;
		if(l == 0)            x = ihx(0)-1.0;
		else if(l == ihx.nel) x = ihx(ihx.nel-1)+1.0;
		else                  x = 0.5*(ihx(l-1)+ihx(l));
		iw.n(l) = cepston(eps(x)*(1.0-r)+twg.eps(x)*r);
	}
	iw.lambda = 0.5*(lambda+twg.lambda);
	return iw;
}

/* refractive index profile -> MATLAB .m file  
   disp: output region on the x axis
   ext0, ext1: filename id characters */
void cWaveguide::plot(Interval disp, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "n__.m";
	name[1] = ext0;
	name[2] = ext1;
	fprintf(stderr, "n [%g, %g] >> %s\n", disp.x0, disp.x1, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "cWaveguide profile");
	name[3] = 0;
	double minn, maxn;
	minn = n.re().min();
	if(minn > n.im().min()) minn = n.im().min();
	maxn = n.re().max();
	if(maxn < n.im().max()) maxn = n.im().max();
	mlout_geo(dat, (*this), minn, maxn); 
	mlout_crefindplot(name, dat, disp); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}
void cWaveguide::plot(char ext0, char ext1) const
{
	Interval disp;
	double w = (hx(nx)-hx(0))/((double) nx);
	disp.x0 = hx(0)-w; 
	disp.x1 = hx(nx)+w; 
	plot(disp, ext0, ext1); 
	return;
}

/* the real waveguide, with permittivity Re(eps), special (!) */ 
Waveguide cWaveguide::re() const
{
	Waveguide wg;
	if(nx == 0)
	{
		wg.nx = 1;
		wg.lambda = lambda;
		wg.special = 1;
		wg.hx = Dvector(2);
		wg.n  = Dvector(3);
		
		double ndd0 = n(0).im;
		double ndd1 = n(1).im;
		if(fabs(ndd1) < fabs(ndd0))
		{
			wg.hx(0) = hx(0);
			wg.hx(1) = hx(0)+LIMBDSEP/2.0;
			wg.n(0) = eps(0).re;
			wg.n(1) = eps(1).re;
			wg.n(2) = eps(1).re;
		}
		else 
		{
			wg.hx(0) = hx(0)-LIMBDSEP/2.0;
			wg.hx(1) = hx(0);
			wg.n(0) = eps(0).re;
			wg.n(1) = eps(0).re;
			wg.n(2) = eps(1).re;
		}
		return wg;
	}
	wg.nx = nx;
	wg.lambda = lambda;
	wg.special = 1;
	wg.hx = hx;
	wg.n = Dvector(nx+2);
	for(int i=0; i<=nx+1; ++i) wg.n(i) = eps(i).re;
	return wg;
}

/* check, whether this is a (real) Waveguide */
bool cWaveguide::isreal() const
{
	if(nx < 1) return false;
	for(int l=0; l<=nx+1; ++l)
	{
		if(fabs(eps(l).im) > COMPTOL_N) return false;
	}
	return true;
}

/* potential SPsolutions: increase in search range for eps_eff */
double CPLXs_EeInc = 10.0;

/* upper bound for effective mode permittivities, default value */
double cWaveguide::defaultepseffmax() const
{
	Waveguide w0 = re();
        double em0 = -AWAY;
        for(int l=0; l<=w0.nx+1; ++l) 
	{ 
		if(em0 < w0.eps(l)) em0 = w0.eps(l);
	}
	double em = em0;
        for(int l=0; l<=w0.nx; ++l) 
	{ 
		double e = w0.eps(l)*w0.eps(l+1)/(w0.eps(l)+w0.eps(l+1));
		if(em < e) em = e; 
	}
	if(em0 < em) return em0+CPLXs_EeInc*(em-em0);
	return em0;
}

/* lower bound for effective mode permittivities, default value */
double cWaveguide::defaultepseffmin() const
{
	Waveguide w0 = re();
	return w0.defaultepseffmin();
}

/* modify the outer regions for guidance, if possible */
void cWaveguide::unleak(double emin)
{
	Waveguide w0 = re();
	if(w0.eps(0) > emin) n(0) = cepston(Complex(emin, 0.0)); 
	if(w0.eps(w0.nx+1) > emin) n(w0.nx+1) = cepston(Complex(emin, 0.0)); 
	return;
}

/* ------------------------------------------------------------------ 
 * modes of complex multilayer slab waveguides
 * ------------------------------------------------------------------ */

/* - internal mode representation: field on a homogeneous layer ------*/

/* initiate: set polarization, interfaces, permittivity */
void cSmPiece::init(Polarization p, Interval i, Complex e, double lambda)
{
	pol = p;
	bif = i.x0; 
	tif = i.x1; 
	k0 = val_k0(lambda);
	eps = e;
	krootrev = 0;
	bdktype = 0;
	return;
}

#define KAPPAQSHIFT 1.0e-16
#define KAPPAHYPOSC 1.0e-10

// translate longitudinal to lateral wavenumbers
void cSmPiece::gqtokappa(Complex gq)
{
	Complex kq = k0*k0*eps-gq;
	if(kq.abs() < KAPPAQSHIFT) kq = Complex(KAPPAQSHIFT, 0.0);	
	if(fabs(kq.im) < KAPPAQSHIFT)  
	{
		kq.im = 0;
		if(krootrev == 1) 
		{
			bdktype = 1;
			if(kq.re >= 0.0)	
			{
				kappa = Complex(-sqrt(kq.re), 0.0);
				return;	
			}
			kappa = Complex(0.0, sqrt(-kq.re));
			return;	
		}
		bdktype = 0;
		if(kq.re >= 0.0)	
		{
			kappa = Complex(sqrt(kq.re), 0.0);
			return;	
		}
		kappa = Complex(0.0, -sqrt(-kq.re));
		return;	
	}
	kappa = sqrt_pos(kq);
	if(krootrev == 1) 
	{
		bdktype = 1;
		if(fabs(kappa.re) > KAPPAQSHIFT)
		{
			if(kappa.re < 0.0) kappa.neg();
		} 
		else 
		{
			kappa.re = 0.0;
			if(kappa.im > 0.0) kappa.neg();
		}
		return;
	}
	bdktype = 0;
	if(fabs(kappa.im) > KAPPAQSHIFT)
	{
		if(kappa.im > 0.0) kappa.neg();
	} 
	else 
	{
		kappa.im = 0.0;
		if(kappa.re < 0.0) kappa.neg();
	}
	return;
}

/* mode solver, transverse resonance evaluation: start */
void cSmPiece::start_tre(Complex gq, Complex& phi, Complex &phis, 
                        Boundary_type bdt, double bp)
{
	gqtokappa(gq);
	Complex C = CC1;
	switch(bdt)
	{
	case OPEN:
		xa = tif; xb = tif;
		A = CC0; 
		B = C;
		phi = C;
		phis = CCI*kappa*C;
		break;
	case LIMD:
		xa = bif; xb = tif;
		A = C*expi( kappa*(bp-xa))/(CCmI*2.0*kappa);
		B = C*expi(-kappa*(bp-xb))/(CCI *2.0*kappa);
		phi  = A*expi(-kappa*(tif-xa))+B*expi(kappa*(tif-xb));
		phis = CCI*kappa*(-A*expi(-kappa*(tif-xa))+B*expi(kappa*(tif-xb)));
		break;
	case LIMN:
		xa = bif; xb = tif;
		A = C*expi( kappa*(bp-xa))*0.5;
		B = C*expi(-kappa*(bp-xb))*0.5;
		phi  = A*expi(-kappa*(tif-xa))+B*expi(kappa*(tif-xb));
		phis = CCI*kappa*(-A*expi(-kappa*(tif-xa))+B*expi(kappa*(tif-xb)));
		break;
	} 
	if(pol == TM) phis /= eps;
	return;
}

/* mode solver, transverse resonance evaluation: continue */
void cSmPiece::cont_tre(Complex gq, Complex& phi, Complex &phis)
{
	if(pol == TM) phis *= eps;
	gqtokappa(gq);
	Complex qm = expi(-kappa*(tif-bif));
	Complex qp = expi( kappa*(tif-bif));
	Complex p = phi;
	Complex ps = phis;
	xa = bif; xb = tif;
	A = (p - ps/(CCI*kappa))*expi( kappa*(bif-xa))*0.5;
	B = (p + ps/(CCI*kappa))*expi(-kappa*(bif-xb))*0.5;
	phi  = ((qm+qp)*p - CCI/kappa*(-qm+qp)*ps)*0.5;
	phis = (CCmI*kappa*(qm-qp)*p + (qm+qp)*ps)*0.5;
	if(pol == TM) phis /= eps;
	return;
}	

/* mode solver, transverse resonance evaluation: finish */
Complex cSmPiece::finish_tre(Complex gq, Complex& phi, Complex &phis, 
                           Boundary_type bdt, double bp)
{
	if(pol == TM) phis *= eps;
	gqtokappa(gq);
	Complex trs = CC0;
	switch(bdt)
	{
	case OPEN:
		xa = bif; xb = bif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(bif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(bif-xb))*0.5;
		trs = B;
		break;
	case LIMD: 
		xa = bif; xb = tif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(bif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(bif-xb))*0.5;
		trs = A*expi(-kappa*(bp-xa))+B*expi(kappa*(bp-xb));
		break;
	case LIMN:
		xa = bif; xb = tif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(bif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(bif-xb))*0.5;
		trs = CCI*kappa*(-A*expi(-kappa*(bp-xa))+B*expi(kappa*(bp-xb)));
		break;
	} 
	return trs;
}

/* set parameters in an upper segment for a modal solution */
void cSmPiece::polish_ub(Boundary_type bdt, double bp)
{
	Complex C;
	switch(bdt)
	{
	case OPEN:
		B = CC0;
		break;
	case LIMD: 
		C = CCmI*2.0*kappa*A*expi(-kappa*(bp-xa)) 
	          + CCI *2.0*kappa*B*expi( kappa*(bp-xb));
		C = C*0.5;
		A = CCI /(2.0*kappa)*C*expi( kappa*(bp-xa));
		B = CCmI/(2.0*kappa)*C*expi(-kappa*(bp-xb));
		break;
	case LIMN:
		C = A*2.0*expi(-kappa*(bp-xa)) 
	          + B*2.0*expi( kappa*(bp-xb));
		C = C*0.5;
		A = C*0.5*expi( kappa*(bp-xa));
		B = C*0.5*expi(-kappa*(bp-xb));
		break;
	} 
	return;
}

/* mode solver, reverse transverse resonance evaluation: start */
void cSmPiece::start_revtre(Complex gq, Complex& phi, Complex &phis, 
                            Boundary_type bdt, double bp)
{
	gqtokappa(gq);
	Complex C = CC1;
	switch(bdt)
	{
	case OPEN:
		xa = bif; xb = bif;
		A = C; 
		B = CC0;
		phi = C;
		phis = CCmI*kappa*C;
		break;
	case LIMD:
		xa = bif; xb = tif;
		A = C*expi( kappa*(bp-xa))/(CCmI*2.0*kappa);
		B = C*expi(-kappa*(bp-xb))/(CCI *2.0*kappa);
		phi  = A*expi(-kappa*(bif-xa))+B*expi(kappa*(bif-xb));
		phis = CCI*kappa*(-A*expi(-kappa*(bif-xa))+B*expi(kappa*(bif-xb)));
		break;
	case LIMN:
		xa = bif; xb = tif;
		A = C*expi( kappa*(bp-xa))*0.5;
		B = C*expi(-kappa*(bp-xb))*0.5;
		phi  = A*expi(-kappa*(bif-xa))+B*expi(kappa*(bif-xb));
		phis = CCI*kappa*(-A*expi(-kappa*(bif-xa))+B*expi(kappa*(bif-xb)));
		break;
	} 
	if(pol == TM) phis /= eps;
	return;
}


/* mode solver, reverse transverse resonance evaluation: continue */
void cSmPiece::cont_revtre(Complex gq, Complex& phi, Complex &phis)
{
	if(pol == TM) phis *= eps;
	gqtokappa(gq);
	Complex qm = expi(-kappa*(tif-bif));
	Complex qp = expi( kappa*(tif-bif));
	Complex p = phi;
	Complex ps = phis;
	xa = bif; xb = tif;
	A = (p - ps/(CCI*kappa))*expi( kappa*(tif-xa))*0.5;
	B = (p + ps/(CCI*kappa))*expi(-kappa*(tif-xb))*0.5;
	phi  = ((qm+qp)*p + CCI/kappa*(-qm+qp)*ps)*0.5;
	phis = ( CCI*kappa*(qm-qp)*p + (qm+qp)*ps)*0.5;
	if(pol == TM) phis /= eps;
	return;
}	

/* mode solver, reverse transverse resonance evaluation: finish */
Complex cSmPiece::finish_revtre(Complex gq, Complex& phi, Complex &phis, 
                                Boundary_type bdt, double bp)
{
	if(pol == TM) phis *= eps;
	gqtokappa(gq);
	Complex trs = CC0;
	switch(bdt)
	{
	case OPEN:
		xa = tif; xb = tif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(tif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(tif-xb))*0.5;
		trs = A;
		break;
	case LIMD: 
		xa = bif; xb = tif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(tif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(tif-xb))*0.5;
		trs = A*expi(-kappa*(bp-xa))+B*expi(kappa*(bp-xb));
		break;
	case LIMN:
		xa = bif; xb = tif;
		A = (phi - phis/(CCI*kappa))*expi( kappa*(tif-xa))*0.5;
		B = (phi + phis/(CCI*kappa))*expi(-kappa*(tif-xb))*0.5;
		trs = CCI*kappa*(-A*expi(-kappa*(bp-xa))+B*expi(kappa*(bp-xb)));
		break;
	} 
	return trs;
}

/* set parameters in a lower segment for a modal solution */
void cSmPiece::polish_lb(Boundary_type bdt, double bp)
{
	Complex C;
	switch(bdt)
	{
	case OPEN:
		A = CC0;
		break;
	case LIMD: 
		C = CCmI*2.0*kappa*A*expi(-kappa*(bp-xa)) 
	          + CCI *2.0*kappa*B*expi( kappa*(bp-xb));
		C = C*0.5;
		A = CCI /(2.0*kappa)*C*expi( kappa*(bp-xa));
		B = CCmI/(2.0*kappa)*C*expi(-kappa*(bp-xb));
		break;
	case LIMN:
		C = A*2.0*expi(-kappa*(bp-xa)) 
	          + B*2.0*expi( kappa*(bp-xb));
		C = C*0.5;
		A = C*0.5*expi( kappa*(bp-xa));
		B = C*0.5*expi(-kappa*(bp-xb));
		break;
	} 
	return;
}

/* field evaluation within a specific layer: principal component */
Complex cSmPiece::bfld(double x) const
{
	Complex f = CC0;
	if(A.abs()>NEGLECT_FLD_CONTRIB) f += A*expmi(kappa*(x-xa));
	if(B.abs()>NEGLECT_FLD_CONTRIB) f += B*expi(kappa*(x-xb));
	return f;
}

/* field evaluation within a specific layer: 
   the derivative of the principal component with respect to x */
Complex cSmPiece::bder(double x) const
{
	Complex f = CC0;
	if(A.abs()>NEGLECT_FLD_CONTRIB) f -= A*expmi(kappa*(x-xa));
	if(B.abs()>NEGLECT_FLD_CONTRIB) f += B*expi(kappa*(x-xb));
	return CCI*f*kappa;
}

/* maximum absolute level of the principal field component within a layer */
double cSmPiece::maxf(Interval i, double &ph) const
{
        return cexp_max(A, B, kappa, xa, xb, i.x0, i.x1, ph);
}

/* integrate the absolute squared principal field component over the entire layer */
double cSmPiece::intbfq() const
{
	return intbfq(bif, tif);
}

#define MININTCONTAMP 1.0e-16

/* integrate the absolute squared principal field component over Interval x0, x1 */
double cSmPiece::intbfq(double x0, double x1) const
{
	Complex s = CC0;
	if(A.sqabs() > MININTCONTAMP)
		s += A.conj()*A*intfcEcE( kappa.conj(), xa, -kappa, xa, x0, x1);
	if((A*B).sqabs() > MININTCONTAMP)
	{
		s += A.conj()*B*intfcEcE( kappa.conj(), xa,  kappa, xb, x0, x1);
		s += B.conj()*A*intfcEcE(-kappa.conj(), xb, -kappa, xa, x0, x1);
	}
	if(B.sqabs() > MININTCONTAMP)
		s += B.conj()*B*intfcEcE(-kappa.conj(), xb,  kappa, xb, x0, x1);
	return s.re;
}

/* integrate the cc of the present principal field component 
   times the principal component of modepiece mp over interval iv */
Complex cSmPiece::intff(const cSmPiece& mp, const Interval& iv) const
{
	Complex s = CC0;
	if((A.conj()*mp.A).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.A*intfcEcE( kappa.conj(), xa, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((A.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.B*intfcEcE( kappa.conj(), xa,  mp.kappa, mp.xb, iv.x0, iv.x1);
	if((B.conj()*mp.A).sqabs() > MININTCONTAMP)
		s += B.conj()*mp.A*intfcEcE(-kappa.conj(), xb, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((B.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += B.conj()*mp.B*intfcEcE(-kappa.conj(), xb,  mp.kappa, mp.xb, iv.x0, iv.x1);
	return s;
}

/* integrate the present principal field component times the derivative
   of the principal component of modepiece mp over interval iv */
Complex cSmPiece::intfd(const cSmPiece& mp, const Interval& iv) const
{
	Complex s = CC0;
	if((A.conj()*mp.A).sqabs() > MININTCONTAMP)
		s -= A.conj()*mp.A*intfcEcE( kappa.conj(), xa, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((A.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.B*intfcEcE( kappa.conj(), xa,  mp.kappa, mp.xb, iv.x0, iv.x1);
	if((B.conj()*mp.A).sqabs() > MININTCONTAMP)
		s -= B.conj()*mp.A*intfcEcE(-kappa.conj(), xb, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((B.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += B.conj()*mp.B*intfcEcE(-kappa.conj(), xb,  mp.kappa, mp.xb, iv.x0, iv.x1);
	return s*CCI*mp.kappa;
}

/* integrate the derivative of the present principal field component
   times the principal component of modepiece mp over interval iv */
Complex cSmPiece::intdf(const cSmPiece& mp, const Interval& iv) const
{
	Complex s = CC0;
	if((A.conj()*mp.A).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.A*intfcEcE( kappa.conj(), xa, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((A.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.B*intfcEcE( kappa.conj(), xa,  mp.kappa, mp.xb, iv.x0, iv.x1);
	if((B.conj()*mp.A).sqabs() > MININTCONTAMP)
		s -= B.conj()*mp.A*intfcEcE(-kappa.conj(), xb, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((B.conj()*mp.B).sqabs() > MININTCONTAMP)
		s -= B.conj()*mp.B*intfcEcE(-kappa.conj(), xb,  mp.kappa, mp.xb, iv.x0, iv.x1);
	return s*CCI*kappa.conj();
}

/* integrate the derivative of the present principal field component times 
   the derivative of the principal component of modepiece mp over interval iv */
Complex cSmPiece::intdd(const cSmPiece& mp, const Interval& iv) const
{
	Complex s = CC0;
	if((A.conj()*mp.A).sqabs() > MININTCONTAMP)
		s -= A.conj()*mp.A*intfcEcE( kappa.conj(), xa, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((A.conj()*mp.B).sqabs() > MININTCONTAMP)
		s += A.conj()*mp.B*intfcEcE( kappa.conj(), xa,  mp.kappa, mp.xb, iv.x0, iv.x1);
	if((B.conj()*mp.A).sqabs() > MININTCONTAMP)
		s += B.conj()*mp.A*intfcEcE(-kappa.conj(), xb, -mp.kappa, mp.xa, iv.x0, iv.x1);
	if((B.conj()*mp.B).sqabs() > MININTCONTAMP)
		s -= B.conj()*mp.B*intfcEcE(-kappa.conj(), xb,  mp.kappa, mp.xb, iv.x0, iv.x1);
	return s*(-1)*kappa.conj()*mp.kappa;
}

/* translate layer representation by dx */
void cSmPiece::translate(double dx)
{ 
	bif += dx;
	tif += dx;
	xa += dx;
	xb += dx;
	return;
}

/* output to file */
void cSmPiece::write(FILE *dat)
{
	comment("(( cSmPiece))", dat);	
	comment("pol", dat); fputint(pol, dat);
	comment("bif", dat); fputdouble(bif, dat);
	comment("tif", dat); fputdouble(tif, dat);
	comment("k0", dat); fputdouble(k0, dat);
	comment("A", dat); A.write(dat);
	comment("B", dat); B.write(dat);
	comment("kappa", dat); kappa.write(dat);
	comment("xa", dat); fputdouble(xa, dat);
	comment("xb", dat); fputdouble(xb, dat);
	comment("eps", dat); eps.write(dat);
	comment("krootrev", dat); fputint(krootrev, dat);
	comment("bdktype", dat); fputint(bdktype, dat);
	return;
}

/* input from file */
void cSmPiece::read(FILE *dat)
{
	pol = Polarization(fgetint(dat));
	bif = fgetdouble(dat);
	tif = fgetdouble(dat);
	k0 = fgetdouble(dat);
	A.read(dat);
	B.read(dat);
	kappa.read(dat);
	xa = fgetdouble(dat);
	xb = fgetdouble(dat);
	eps.read(dat);
	krootrev = fgetint(dat);
	bdktype = fgetint(dat);
	return;
}

/* - slab mode objects ---------------------------------------------------- */

/* initialize */
cMode::cMode()
{
	N = 0;
	piece = NULL;
	special = 0;
	cst = -1;
	orgM.piece = NULL;
}
cMode::cMode(const cWaveguide& g, 
             Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt,
             Polarization p)
{
	wg = g;
	pol = p;
	k0 = val_k0(wg.lambda);
	invommu0 = val_invommu0(wg.lambda);
	invomep0 = val_invomep0(wg.lambda);
	N = wg.nx+2;
	bdt_b = bdtb;
	bdt_t = bdtt;
	bp_b = bpb;
	bp_t = bpt;
	Interval l;
	double dx, wl = wg.lambda;
	piece = new cSmPiece[N];
	switch(bdtb)
	{
	case OPEN:
		l.x0 = -AWAY; l.x1 = wg.hx(0);
		piece[0].init(p, l, wg.eps(0), wl);
		break;
	case LIMD: case LIMN:
		dx = wg.hx(0)-bpb;
		l.x0 = wg.hx(0)-dx-dx; l.x1 = wg.hx(0);
		piece[0].init(p, l, wg.eps(0), wl);
		break;
	}
	for(int i=1; i<=N-2; ++i) piece[i].init(p, wg.layer(i), wg.eps(i), wl);	
	switch(bdtt)
	{
	case OPEN:
		l.x0 = wg.hx(wg.nx); l.x1 = AWAY;
		piece[N-1].init(p, l, wg.eps(wg.nx+1), wl);
		break;
	case LIMD: case LIMN:
		dx = bpt-wg.hx(wg.nx);
		l.x0 = wg.hx(wg.nx); l.x1 = wg.hx(wg.nx)+dx+dx;
		piece[N-1].init(p, l, wg.eps(wg.nx+1), wl);
		break;
	}
	special = 0;
	cst = -1;
	orgM.piece = NULL;
}
cMode::cMode(const Mode m)
{ 
	cMode mc(cWaveguide(m.wg), m.bdt_b, m.bp_b, m.bdt_t, m.bp_t, m.pol);
	Complex gq; 
	gq = Complex(m.betaq, 0.0);
	mc.ord = m.ord;
	mc.special = m.special;
	mc.polish(gq);
	mc.orgM = m;
	wg = mc.wg;	
	pol = mc.pol;	
	k0 = mc.k0;	
	invommu0 = mc.invommu0;	
	invomep0 = mc.invomep0;	
	gammaq = mc.gammaq;
	gamma = mc.gamma;
	neff = mc.neff;
	beta = mc.beta;
	alpha = mc.alpha;
	ord = mc.ord;
	bdt_b = mc.bdt_b;
	bdt_t = mc.bdt_t;
	bp_b = mc.bp_b;
	bp_t = mc.bp_t;
	maxF = mc.maxF;
	N = mc.N;
	piece = new cSmPiece[N];
	for(int i=0; i<=N-1; ++i) piece[i] = mc.piece[i];
	for(int i=0; i<=7; ++i) ids[i] = mc.ids[i];
	for(int i=0; i<=7; ++i) btids[i] = mc.btids[i];
	special = mc.special;
	orgM = mc.orgM;
	cst = 0;
}

/* destroy */
cMode::~cMode()
{
	if(piece != NULL) delete[] piece;
	piece = NULL;
	N = 0;
	cst = -1;	
	if(orgM.piece != NULL) delete[] orgM.piece;
	orgM.piece = NULL;
}

/* copy constructor */
cMode::cMode(const cMode& m)
{
	wg = m.wg;	
	pol = m.pol;	
	k0 = m.k0;	
	invommu0 = m.invommu0;	
	invomep0 = m.invomep0;	
	gammaq = m.gammaq;
	gamma = m.gamma;
	neff = m.neff;
	beta = m.beta;
	alpha = m.alpha;
	ord = m.ord;
	bdt_b = m.bdt_b;
	bdt_t = m.bdt_t;
	bp_b = m.bp_b;
	bp_t = m.bp_t;
	maxF = m.maxF;
	N = m.N;
	piece = new cSmPiece[N];
	for(int i=0; i<=N-1; ++i) piece[i] = m.piece[i];
	for(int i=0; i<=7; ++i) ids[i] = m.ids[i];
	for(int i=0; i<=7; ++i) btids[i] = m.btids[i];
	special = m.special;
	orgM = m.orgM;
	cst = m.cst;
}

/* assignment */
cMode& cMode::operator=(const cMode& m)
{
	if(this != &m)
	{
		wg = m.wg;	
		pol = m.pol;	
		k0 = m.k0;	
		invommu0 = m.invommu0;	
		invomep0 = m.invomep0;	
		gammaq = m.gammaq;
		gamma = m.gamma;
		neff = m.neff;
		beta = m.beta;
		alpha = m.alpha;
		ord = m.ord;
		bdt_b = m.bdt_b;
		bdt_t = m.bdt_t;
		bp_b = m.bp_b;
		bp_t = m.bp_t;
		maxF = m.maxF;
		if(piece != NULL) delete[] piece;
		N = m.N;
		piece = new cSmPiece[N];
		for(int i=0; i<=N-1; ++i) piece[i] = m.piece[i];
		for(int i=0; i<=7; ++i) ids[i] = m.ids[i];
		for(int i=0; i<=7; ++i) btids[i] = m.btids[i];
		special = m.special;
		orgM = m.orgM;
		cst = m.cst;
	}
	return *this;
}

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

/* mode solver: transverse resonance evaluation */
Complex cMode::travres(Complex tgq)
{
	Complex phi, phis;
	piece[0].start_tre(tgq, phi, phis, bdt_b, bp_b); 
	for(int l=1; l<=N-2; ++l) piece[l].cont_tre(tgq, phi, phis); 
	Complex r = piece[N-1].finish_tre(tgq, phi, phis, bdt_t, bp_t); 
// fprintf(stderr, "\n   TR("); tgq.nice(stderr); fprintf(stderr, ") = %g\n", r.sqabs());
	return r;
}

/* mode solver: reverse transverse resonance evaluation */
Complex cMode::revtravres(Complex tgq)
{
	Complex phi, phis;
	piece[N-1].start_revtre(tgq, phi, phis, bdt_t, bp_t); 
	for(int l=N-2; l>=1; --l) piece[l].cont_revtre(tgq, phi, phis); 
	return piece[0].finish_revtre(tgq, phi, phis, bdt_b, bp_b); 
}

/* mode solver: set the relevant variables for an actual modal solution */ 
void cMode::polish(Complex gq)
{
	gammaq = gq;
	gamma = sqrt_pos(gq);
	if(fabs(gamma.re) < CPLXs_GqTol)
	{
		gamma.re = 0.0;
		if(gamma.im > 0.0) gamma *= -1.0;
	}
	else 
	{
		if(gamma.re < 0.0) gamma *= -1.0;
	}
	neff = gamma/k0;
	beta = gamma.re;
	alpha = -gamma.im;
	
	gammaq = gamma*gamma;

	travres(gammaq);
	piece[N-1].polish_ub(bdt_t, bp_t);
	normalize();
	classify();
	setids();

	cMode m = (*this);
	m.revtravres(gammaq);
	m.piece[0].polish_lb(bdt_b, bp_b);
	m.normalize();
	m.classify();
	m.setids();

 	if(m.checkcontinuity() < checkcontinuity()) (*this) = m;

	return;
}

/* propagation length */
double cMode::Lp()
{
	if(alpha <= 0.0) return INFINITY;
	return 0.5/alpha;
}

// secant method, terminal distance of root prototypes
double CPLXs_RootAcc = 1.0e-12;
// number of trials for convergence 
int CPLXs_LimNumIts = 20;
// maximum growth of trial value distance permitted per iteration step
double CPLXs_ConvGrowth = 2.0;
// mimimum distance in gammaq for resonances to be considered different
double CPLXs_MinGqDiff = 1.0e-8;
// maximum acceptable level of absolute TRC violation 
double CPLXs_MaxTRC = 1.0e-6;

/* mode solver: trace propagation constants in the complex plane,
   converge a single step, return 1: success, 0: failure,
   kr_b, kr_t; reverse sign of kappaq-root, at bottom / top boundary, if applicable, 
   if OK: ftrE: final transverse resonance evaluation, remaining error */ 
int cMode::traceto(Complex egq, int kr_b, int kr_t, int quiet, double &ftrE)
{
	ftrE = 1.0;
	piece[0].krootrev = 0;
	if(bdt_b == OPEN && kr_b != 0) piece[0].krootrev = 1;
	piece[N-1].krootrev = 0;
	if(bdt_t == OPEN && kr_t != 0) piece[N-1].krootrev = 1;
	Complex z1 = gammaq;
	Complex z2 = egq;
	Complex f1, f2, f;
	Complex z = (z1+z2)*0.5;
if(quiet < 1)
{
	fprintf(stderr, "\n[%s_%d^%d, gq=", ids, kr_b, kr_t); z1.nice(stderr); fprintf(stderr, " egq="); z2.nice(stderr); fprintf(stderr, " |z2-z1| = %g ]\n", (z2-z1).abs());  
}
	int trcE = 0;
	f1 = travres(z1); if(f1.isnotOK()) ++trcE;
	f2 = travres(z2); if(f2.isnotOK()) ++trcE;
	if(f1.abs() < f2.abs())
	{
		f  = f2; z  = z2; 
		f2 = f1; z2 = z1; 
		f1 = f;  z1 = z; 
	}
if(quiet < 1) fprintf(stderr, "   i: "); 
	int nit = 0;
	while(trcE == 0 && ((z1-z2).abs() > CPLXs_RootAcc) && nit < CPLXs_LimNumIts)
	{
		z = (z2*f1-z1*f2)/(f1-f2);
if(quiet < 1) fprintf(stderr, "c"); 
		if((z-z2).abs() > CPLXs_ConvGrowth*(z2-z1).abs())
		{
			z = z2 + CPLXs_ConvGrowth*((z2-z1).abs()/(z-z2).abs())*(z-z2);
if(quiet < 1) fprintf(stderr, "g"); 
		}
		f = travres(z);  if(f.isnotOK()) ++trcE;
		f1 = f2; z1 = z2; 
		f2 = f;  z2 = z; 
		if(f1.abs() < f2.abs())
		{
if(quiet < 1) fprintf(stderr, "r"); 
			f  = f2; z  = z2; 
			f2 = f1; z2 = z1; 
			f1 = f;  z1 = z; 
		}
if(quiet < 1) fprintf(stderr, " "); 
		++nit;
	}
if(quiet < 1) fprintf(stderr, "\n   "); 
	if(trcE > 0)
	{
if(quiet < 1) fprintf(stderr, "F.\n"); 
	}
// if(quiet < 1) fprintf(stderr, "trcE: %d, |dz| = %g, Im(z) = %g, |TRC| = %g\n", trcE, (z1-z2).abs(), z2.im, f2.sqabs()); 
	if(trcE == 0 && ((z1-z2).abs() <= CPLXs_RootAcc) && f2.sqabs() < CPLXs_MaxTRC)
	{
		polish(z2);
if(quiet < 1)
{
		fprintf(stderr, "[-> gq="); gammaq.nice(stderr); fprintf(stderr, " errTRC: %.14g ] ", f2.sqabs());
}
		if(trcE == 0)
		{
if(quiet < 1) fprintf(stderr, ".\n");
			ftrE = f2.abs();
			return 1;
		}
		else
		{
if(quiet < 1) fprintf(stderr, " (-).\n");
			return 0;
		}
	}
if(quiet < 1) fprintf(stderr, "[-> -.]\n"); 
	return 0;
}

/* mode solver: perturbationally estimate the phase shift of the real basic mode, 
   using the permittivity difference to a complex target waveguide */
Complex estimategammaq(Mode mr, cWaveguide twg)
{
	cMode m = cMode(mr);
	cWaveguide iwg = m.wg.intermediate(twg, 0.5);
	Complex g = m.gamma;
	for(int l=0; l<=iwg.nx+1; ++l)
	{
		Cmatrix deps;
		deps.unit(3);
		double x = 0.5*(iwg.layer(l).x0+iwg.layer(l).x1);
		deps = mult(deps, (twg.eps(x)-mr.wg.eps(x)));
		Perturbation p;
		p.pr = iwg.layer(l);
		p.e = deps;
		g += mr.phaseshift(p, FORW);
	}
	return g*g;
}

/* complexification & leakyfication: 
   maximum number of steps in waveguide change */
double CPLXs_WgSteps = 16;

/* mode solver: trace propagation constants in the complex plane */
cModeArray complexify(const cMode o, const cWaveguide twg, int quiet)
{
	cModeArray ma;
	ma.clear();
	for(int nst=1; nst<=CPLXs_WgSteps; nst*=2)
        {
		cModeArray t = complexify(o, twg, nst, quiet);
                ma.merge(t);
        }
	ma.prune();
        return ma;
}
cModeArray complexify(const cMode o, const cWaveguide twg, int nst, int quiet)
{
	cModeArray ma;
	ma.clear();
// fprintf(stderr, "\ncplxfy: %d\n", nst);
	cWaveguide iwg;
	if(o.wg.equal(twg) == 1)
	{
		ma.add(o); 
		return ma;
	}
	double rst = 1.0/((double) nst);
	Complex ogq = o.gammaq; 
	Complex egq = estimategammaq(o.orgM, twg);
	Complex tgq = ogq+(egq-ogq)*rst;
	if((tgq-ogq).abs() < CPLXs_RootAcc)
	{
		Complex g = o.gamma+Complex(0.0, -0.01)*o.k0;	
		tgq = ogq+(g*g-ogq)*rst;
	}
	cMode m;
	for(int i=1; i<=nst; ++i)
	{
		double r = ((double) i)*rst;
		iwg = o.wg.intermediate(twg, r);
		m = cMode(iwg, o.bdt_b, o.bp_b, o.bdt_t, o.bp_t, o.pol);
		m.special = o.special;
		m.ord = o.ord;
		m.orgM = o.orgM;
		m.polish(ogq); 
// fprintf(stderr, "\ntrace: "); ogq.nice(stderr);
		double dummy;
		int ok = m.traceto(tgq, 0, 0, 1, dummy); 
		if(ok == 0)
		{
// fprintf(stderr, " -> NOK.");
			 return ma;
		}
// fprintf(stderr, " -> "); m.gammaq.nice(stderr); fprintf(stderr, "\n");
		tgq = ogq+(m.gammaq-ogq)*2.0;
		if((tgq-m.gammaq).abs() < CPLXs_RootAcc) 
		{
			Complex g = m.gamma+Complex(0.001, -0.0001)*m.k0;	
			tgq = m.gammaq+(g*g-m.gammaq)*rst;
		}
		ogq = m.gammaq;
	}
	if(ogq.abs() < CPLXs_RootAcc)
	{ 
// fprintf(stderr, " -> NOK.");
		return ma;
	}
	if(fabs(ogq.re) < CPLXs_RootAcc) ogq.re = 0.0;
	if(fabs(ogq.im) < CPLXs_RootAcc) ogq.im = 0.0;
	m.polish(ogq);
	m.cst = 1;
	ma.add(m);
	return ma;
} 

cModeArray leakify(const cMode o, const cWaveguide twg, int quiet)
{
	cModeArray ma;
	ma.clear();
	int nst = 1;
	while(ma.num <= 0 && nst <= CPLXs_WgSteps)
	{
		ma = leakify(o, twg, nst, quiet);
		nst = nst*2;
	}
	return ma;
}

/* leaky mode solver: trace propagation constants in the complex plane */
cModeArray leakify(const cMode o, const cWaveguide twg, int nst, int quiet)
{
	double k0 = val_k0(twg.lambda);
	cModeArray tma; 
	tma.clear();
	tma.add(o);
	if(o.wg.equal(twg) == 1)
	{
		return tma;
	}
	Complex ogq = o.gammaq; 
if(quiet == 0) fprintf(stderr, "lkf(%s) #%d:\n", o.ids, nst); 
	Complex egq = estimategammaq(o.orgM, twg);
	if((egq-ogq).abs() < CPLXs_RootAcc) 
	{
if(quiet == 0) 
{
		fprintf(stderr, "(%d) %s, neff = ", 0, tma(0).ids); tma(0).neff.nice(stderr); fprintf(stderr, " (d))\n"); 
}
		return tma;
	}
	Complex dgq = egq-ogq;
	double rst = 1.0/((double) nst);
	dgq = dgq*rst;

	for(int i=1; i<=nst; ++i)
	{
		double r = ((double) i)*rst;
		cWaveguide iwg = o.wg.intermediate(twg, r);

		cMode m;
		m = cMode(iwg, o.bdt_b, o.bp_b, o.bdt_t, o.bp_t, o.pol);
		m.special = o.special;
		m.ord = o.ord;
		m.orgM = o.orgM;

		for(int j=0; j<=tma.num-1; ++ j)
		{
			cModeArray nma; 
			nma.clear();

			int ok; 
			double ftrE;
			ogq = tma(j).gammaq;
			int nss = 0;
			double sof[10] = {0.0, 1.0, -1.0, 2.0, -2.0, 5.0, -5.0, 10.0, -10.0, 20.0};
			Complex dgqi(0.0, -0.00001);
			while(nma.num <= 0 && nss <= 9)
			{
				Complex gqa, gqb;
				gqa = ogq+sof[nss]*(dgq+dgqi);
				gqb = gqa+dgq+dgqi;
				m.polish(gqa); 
				ok = m.traceto(gqb, 0, 0, 1, ftrE); 
				if(ok == 1) nma.add(m);
				int ktb, ktt;
				double eef;

				eef = gqa.re/k0/k0;
				if(eef >= iwg.eps(0).re) ktb = 0;
				else                     ktb = 1;
				if(eef >= iwg.eps(iwg.nx+1).re) ktt = 0;
				else                            ktt = 1;
				if(ktb != 0 || ktt != 0)
				{
					m.polish(gqa); 
					ok = m.traceto(gqb, ktb, ktt, 1, ftrE); 
					if(ok == 1) nma.add(m);
				}

				eef = gqb.re/k0/k0;
				if(eef >= iwg.eps(0).re) ktb = 0;
				else                     ktb = 1;
				if(eef >= iwg.eps(iwg.nx+1).re) ktt = 0;
				else                            ktt = 1;
				if(ktb != 0 || ktt != 0)
				{
					m.polish(gqa); 
					ok = m.traceto(gqb, ktb, ktt, 1, ftrE);
					if(ok == 1) nma.add(m);
				}
				++nss;
			}	
			nma.sort();
			nma.prune();
			tma = nma;
		}	
		if(tma.num <= 0)
		{
			tma.clear();
			return tma;
		}
	}
	for(int j=0; j<=tma.num-1; ++j)
	{
		int i=0;
		tma(j).ids[i++] = 'T';
		if(tma(j).pol == TE) tma(j).ids[i++] = 'E';
		else                 tma(j).ids[i++] = 'M';
		if(tma(j).special != 0) tma(j).ids[i++] = 'x';
		tma(j).ids[i++] = '(';
		if(tma(j).piece[0].bdktype == 0) tma(j).ids[i++] = 'b';
		else                             tma(j).ids[i++] = 'l';
		if(tma(j).piece[tma(j).wg.nx+1].bdktype == 0) tma(j).ids[i++] = 'b';
		else                                          tma(j).ids[i++] = 'l';
		tma(j).ids[i++] = ')';
		while(i<=7) tma(j).ids[i++] = 0;
		i=0;
		if(tma(j).piece[0].bdktype == 0) tma(j).btids[i++] = 'b';
		else                             tma(j).btids[i++] = 'l';
		if(tma(j).piece[tma(j).wg.nx+1].bdktype == 0) tma(j).btids[i++] = 'b';
		else                                          tma(j).btids[i++] = 'l';
		while(i<=7) tma(j).btids[i++] = 0;
		tma(j).ord = -1;
		tma(j).cst = 2;
if(quiet == 0)
{
	fprintf(stderr, "(%d) %s, neff = ", j, tma(j).ids); tma(j).neff.nice(stderr); fprintf(stderr, "\n"); 
}
	}
	return tma;
} 

/* expand the present mode object with LIMN or LIMD condition
   at the upper boundary bp_t towards a symmetric (LIMN) or
   antisymmetric (LIMD) mode on the doubled waveguide */
void cMode::expand()
{
	if(bdt_t == OPEN) cplxwgerror("cMode: expand: bdt");
	if(orgM.piece != NULL)
	{
		if(orgM.bdt_t != bdt_t) cplxwgerror("cMode: orgM bdt mismatch");
		orgM.expand();
	}
	cWaveguide nwg;	
	nwg = wg.expand(bp_t);
	wg = nwg;
	int sym = 1;
	if(bdt_t == LIMD) sym = -1;
	double p = bp_t;
	bdt_t = bdt_b;
	bp_t = p+(p-bp_b);
	if(sym == 1) ord *= 2;
	else         ord = 2*ord+1;

	int nN = N+N-1;
	cSmPiece *npiece;
	npiece = new cSmPiece[nN];
	for(int l=0; l<=N-1; ++l) npiece[l] = piece[l];	
	for(int l=0; l<=N-2; ++l) 
	{
		int tl = nN-1-l;
		npiece[tl].pol   = piece[l].pol;
		npiece[tl].bif   = p+(p-piece[l].tif);
		npiece[tl].tif   = p+(p-piece[l].bif);
		npiece[tl].k0    = piece[l].k0;
		npiece[tl].kappa = piece[l].kappa;
		npiece[tl].eps   = piece[l].eps;
		npiece[tl].krootrev = piece[l].krootrev;
		npiece[tl].bdktype = piece[l].bdktype;
		if(sym == 1)
		{
			npiece[tl].A  = piece[l].B;
			npiece[tl].xa = p+(p-piece[l].xb);
			npiece[tl].B  = piece[l].A;
			npiece[tl].xb = p+(p-piece[l].xa);
		}
		else
		{
			npiece[tl].A  = -piece[l].B;
			npiece[tl].xa = p+(p-piece[l].xb);
			npiece[tl].B  = -piece[l].A;
			npiece[tl].xb = p+(p-piece[l].xa);
		}
	}
	delete[] piece;
	N = nN;
	piece = npiece;

	normalize();
	setids();
	return;
}

/* mirror the present mode object with OPEN condition
   at the upper boundary with respect to the boundary position p */
void cMode::mirror(double p)
{
	if(bdt_t != OPEN) cplxwgerror("cMode: mirror: bdt");
	if(p < wg.hx(wg.nx)+HDIST) cplxwgerror("cMode: mirror: p");

	cWaveguide nwg;	
	nwg = wg.mirror(p);
	wg = nwg;
	bdt_t = bdt_b;
	bp_t = p+(p-bp_b);
	bdt_b = OPEN;
	bp_b = -AWAY;

	cSmPiece *npiece;
	npiece = new cSmPiece[N];
	for(int l=0; l<=N-1; ++l) 
	{
		int tl = N-1-l;
		npiece[tl].pol   = piece[l].pol;
		npiece[tl].bif   = p+(p-piece[l].tif);
		npiece[tl].tif   = p+(p-piece[l].bif);
		npiece[tl].k0    = piece[l].k0;
		npiece[tl].kappa = piece[l].kappa;
		npiece[tl].eps   = piece[l].eps;
		npiece[tl].krootrev = piece[l].krootrev;
		npiece[tl].bdktype = piece[l].bdktype;

		npiece[tl].A  = piece[l].B;
		npiece[tl].xa = p+(p-piece[l].xb);
		npiece[tl].B  = piece[l].A;
		npiece[tl].xb = p+(p-piece[l].xa);
	}
	delete[] piece;
	piece = npiece;
	return;
}

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

/* classify: determine mode order, if possible,
             otherwise: set the special switch */
#define NONODESLIM 1.0e-6
void cMode::classify()
{
	return;
}

/* set the mode id string */
void cMode::setids()
{
	int i=0;
	ids[i++] = 'T';
	if(pol == TE) ids[i++] = 'E';
	else          ids[i++] = 'M';
	if(ord >=100) ids[i++] = dig100(ord); 
	if(ord >=10)  ids[i++] = dig10(ord); 
	ids[i++] = dig1(ord); 
	if(special != 0) ids[i++] = 'x';
	while(i<=7) ids[i++] = 0;
	i=0;
	btids[i++] = '-';
	while(i<=7) btids[i++] = 0;
	return;
}
/* ... for known mode order */
void cMode::setids(int o)
{
	ord = o;
	setids();
	return;
}

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

/* minimum distance between boundaries to be considered separated */
#define S_HDIST 1.0e-10

/* mode profile,
   the principal field component at transverse coordinate x */
Complex cMode::field(double x) const
{
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return CC0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return CC0;
	return piece[wg.layeridx(x)].bfld(x); 
}
/* and its derivative */
Complex cMode::d_field(double x) const
{
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return CC0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return CC0;
	return piece[wg.layeridx(x)].bder(x); 
}

/* field values at point x, cp: EX - HZ, SX, SY, SZ, the forward field */
Complex cMode::field(Fcomp cp, double x) const
{
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return CC0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return CC0;
	int l = wg.layeridx(x);
	Complex f, d;
	switch(pol)
	{
	case TE:
		switch(cp)
		{
		case EX: case EZ: case HY:
			return CC0;
			break;
		case EY:
			return piece[l].bfld(x);
			break;
		case HX: 
			return -gamma*piece[l].bfld(x)*invommu0;
			break;
		case HZ:
			return CCI*piece[l].bder(x)*invommu0;
			break;
		case SX:
			f = piece[l].bfld(x);
			d = piece[l].bder(x);
			return Complex((f.conj()*CCI*d).re*0.5*invommu0, 0.0);
			break;
		case SY:
			return CC0;
			break;
		case SZ:
			f = piece[l].bfld(x);
			return Complex((f.conj()*gamma*f).re*0.5*invommu0, 0.0);
			break;
		default:
			cplxwgerror("cMode: field, cp not implemented");
			return CC0;
			break;
		}
		break;
	case TM:
		switch(cp)
		{
		case HX: case HZ: case EY:
			return CC0;
			break;
		case HY:
			return piece[l].bfld(x);
			break;
		case EX: 
			return gamma*piece[l].bfld(x)*invomep0/wg.eps(l);
			break;
		case EZ:
			return CCmI*piece[l].bder(x)*invomep0/wg.eps(l);
			break;
		case SX:
			f = piece[l].bfld(x);
			d = piece[l].bder(x);
			return Complex(-(CCI*d.conj()*f/wg.eps(l).conj()).re*0.5*invomep0, 0.0);
			break;
		case SY:
			return CC0;
			break;
		case SZ:
			f = piece[l].bfld(x);
			return Complex((gamma.conj()*f.conj()*f/wg.eps(l).conj()).re*0.5*invomep0, 0.0);
			break;
		default:
			cplxwgerror("cMode: field, cp not implemented");
			return CC0;
			break;
		}
		break;
	default: 
		cplxwgerror("cMode: field, pol?");
		return CC0;
	}
	cplxwgerror("cMode: field, confused.");
	return CC0;
}

/* field values for forward or backward propagating modes, Fcomp: EX - HZ */
Complex cMode::field(Fcomp cp, Propdir d, double x) const
{
	if(((int) cp) > ((int) HZ)) cplxwgerror("cMode: field, Propdir: cp");
	Complex f = field(cp, x);
	if(d == FORW) return f;
	switch(cp) 
	{
		case EX: case EY: case HZ: 
			return f;
			break;
		case EZ: case HX: case HY:
			return -f;
			break;
		default:
			cplxwgerror("cMode: field, dir: cp not implemented");
			return CC0;
			break;
	}
	cplxwgerror("cMode: field, Propdir: ?");
	return f; 
}

/* factor of phase/amplitude change along a propagation distance z
   d: FORW, BACK selects the forward / backward traveling field */
Complex cMode::efac(Propdir d, double z) const
{
	if(d == FORW) return expmi(gamma*z);
	return expi(gamma*z);
}

/* mode field evolution along z, given a unit amplitude at z=0
   d: FORW, BACK selects the forward / backward traveling field 
   Fcomp: EX - HZ */
Complex cMode::field(Fcomp cp, Propdir d, double x, double z) const
{
	return field(cp, d, x)*efac(d, z);
}


/* check the continuity of the principal field components,
   returns a maximum error quantity */
double cMode::checkcontinuity() const
{
	double xm, xp;
	Complex p, m;
	double d, dmax = 0.0;
	for(int a=0; a<=wg.nx; ++a)
	{
		xp = wg.hx(a)+1.0e-12;
		xm = wg.hx(a)-1.0e-12;
		m = piece[a].bfld(xm);
		p = piece[a+1].bfld(xp);
		d = (p-m).abs()/maxF;
		if(d > dmax) dmax = d;
	} 
	return dmax;
}

/* check the validity of the Maxwell-equations at a number of points */
#define DERIVSTEP 1.0e-6
void cMode::checkmode() const
{
	Complex ex, ey, ez, hx, hy, hz, gd;
	Complex s1, s2, s3, s4, s5, s6, s7, s8;
	double s1sum, s2sum, s3sum, s4sum, s5sum, s6sum, s7sum, s8sum, s9sum;
	Complex dxex, dxey, dxez, dxhx, dxhy, dxhz;
	double x0, x1, xp, s, ds;
	double ommu, omep;
	Complex eps;
	Complex ddxxbf, s9, bf, dxbf;

	ommu = 1.0/val_invommu0(wg.lambda);
	omep = 1.0/val_invomep0(wg.lambda);
	fprintf(stderr, "Modecheck:\n");
	fprintf(stderr, "gamma: %g +i %g,  neff: %g +i %g\n", gamma.re, gamma.im, neff.re, neff.im);
	fprintf(stderr, "beta: %g, alpha: %g\n", beta, alpha);
	fprintf(stderr, "pol: T%c,  ", polCHR(pol) );

	if(bdt_b == OPEN) fprintf(stderr, "bdt_b: OPEN,  ");
	if(bdt_b == LIMD) fprintf(stderr, "bdt_b: LIMD,  ");
	if(bdt_b == LIMN) fprintf(stderr, "bdt_b: LIMN,  ");
	fprintf(stderr, "bp_b: %g\n", bp_b);
	if(bdt_t == OPEN) fprintf(stderr, "bdt_t: OPEN,  ");
	if(bdt_t == LIMD) fprintf(stderr, "bdt_t: LIMD,  ");
	if(bdt_t == LIMN) fprintf(stderr, "bdt_t: LIMN,  ");
	fprintf(stderr, "bp_t: %g\n", bp_t);

	fprintf(stderr, "Absolute field maximum, principal component:  %g\n", maxF);
	
	for(int d=0; d<=1; ++d)
	{
		Propdir dir;
		if(d==0) { dir = FORW; fprintf(stderr, "--- F ---\n"); }
		else     { dir = BACK; fprintf(stderr, "--- B ---\n"); }
	if(bdt_b == OPEN) x0 = wg.hx(0)-10.0;
	else              x0 = bp_b;
	ex = field(EX, dir, x0);
	ey = field(EY, dir, x0);
	ez = field(EZ, dir, x0);
	hx = field(HX, dir, x0);
	hy = field(HY, dir, x0);
	hz = field(HZ, dir, x0);
	if(bdt_b == OPEN) fprintf(stderr, "Absolute fields at x = %g:\n", x0);
	else fprintf(stderr, "Absolute fields at lower boundary:\n");
	fprintf(stderr, "Ex = %g,  Ey = %g, Ez = %g\n", 
		ex.abs(), ey.abs(), ez.abs());
	fprintf(stderr, "Hx = %g,  Hy = %g, Hz = %g\n", 
		hx.abs(), hy.abs(), hz.abs());

	if(bdt_t == OPEN) x1 = wg.hx(wg.nx)+10.0;
	else              x1 = bp_t;
	ex = field(EX, dir, x1);
	ey = field(EY, dir, x1);
	ez = field(EZ, dir, x1);
	hx = field(HX, dir, x1);
	hy = field(HY, dir, x1);
	hz = field(HZ, dir, x1);
	if(bdt_t == OPEN) fprintf(stderr, "Absolute fields at x = %g:\n", x1);
	else fprintf(stderr, "Absolute fields at upper boundary:\n");
	fprintf(stderr, "Ex = %g,  Ey = %g, Ez = %g\n", 
		ex.abs(), ey.abs(), ez.abs());
	fprintf(stderr, "Hx = %g,  Hy = %g, Hz = %g\n", 
		hx.abs(), hy.abs(), hz.abs());

	fprintf(stderr, "Check of Maxwells equations:\n");
	gd = gamma;
	if(dir == BACK) gd = -gamma;
	if(bdt_b == OPEN) x0 = wg.hx(0)-2.0; 
	if(bdt_t == OPEN) x1 = wg.hx(wg.nx)+2.0; 
	s = (x1-x0)/200.0;	
	ds = DERIVSTEP/2.0;
	s1sum = s2sum = s3sum = s4sum = s5sum = s6sum = s7sum = s8sum = s9sum = 0.0;
	for(xp=x0+s; xp<=x1-s; xp+=s)
	{
		if((wg.eps(xp+ds)-wg.eps(xp-ds)).abs()<1.0e-12)
		{
		eps = wg.eps(xp); 
		ex = field(EX, dir, xp);
		ey = field(EY, dir, xp);
		ez = field(EZ, dir, xp);
		hx = field(HX, dir, xp);
		hy = field(HY, dir, xp);
		hz = field(HZ, dir, xp);

		dxex = (field(EX, dir, xp+ds)-field(EX, dir, xp-ds))/(2.0*ds);
		dxey = (field(EY, dir, xp+ds)-field(EY, dir, xp-ds))/(2.0*ds);
		dxez = (field(EZ, dir, xp+ds)-field(EZ, dir, xp-ds))/(2.0*ds);
		dxhx = (field(HX, dir, xp+ds)-field(HX, dir, xp-ds))/(2.0*ds);
		dxhy = (field(HY, dir, xp+ds)-field(HY, dir, xp-ds))/(2.0*ds);
		dxhz = (field(HZ, dir, xp+ds)-field(HZ, dir, xp-ds))/(2.0*ds);

		s1 = CCmI*gamma*ey-CCI*hx*ommu; 
		s2 = dxez+CCI*gamma*ex-CCI*hy*ommu; 
		s3 = -dxey-CCI*hz*ommu; 
		s4 = CCmI*gamma*hy+CCI*ex*omep*eps;
		s5 = dxhz+CCI*gamma*hx+CCI*ey*omep*eps;
		s6 = -dxhy+CCI*ez*omep*eps;
		s7 = dxex-CCI*gamma*ez;
		s8 = dxhx-CCI*gamma*hz;
		
		if(s1.abs() > s1sum) s1sum = s1.abs();
		if(s2.abs() > s2sum) s2sum = s2.abs();
		if(s3.abs() > s3sum) s3sum = s3.abs();
		if(s4.abs() > s4sum) s4sum = s4.abs();
		if(s5.abs() > s5sum) s5sum = s5.abs();
		if(s6.abs() > s6sum) s6sum = s6.abs();
		if(s7.abs() > s7sum) s7sum = s7.abs();
		if(s8.abs() > s8sum) s8sum = s8.abs();

		if(pol == TE)
		{
			bf = field(EY, xp);
			ddxxbf = (field(EY, xp+ds)
			         -2.0*field(EY, xp)
			         +field(EY, xp-ds))/(ds*ds);
		}	
		else
		{
			bf = field(HY, xp);
			ddxxbf = (field(HY, xp+ds)
			         -2.0*field(HY, xp)
			         +field(HY, xp-ds))/(ds*ds);
		}
		s9 = ddxxbf+(k0*k0*eps-gamma*gamma)*bf;
		if(s9.abs() > s9sum) s9sum = s9.abs();
		}
	} 
	fprintf(stderr, "curl E =  - i om mu0 H:  (%g, %g, %g)\n", 
		s1sum, s2sum, s3sum);
	fprintf(stderr, "curl H = i om ep0 ep E:  (%g, %g, %g)\n", 
		s4sum, s5sum, s6sum);
	fprintf(stderr, "div E = 0:               %g\n", s7sum);
	fprintf(stderr, "div H = 0:               %g\n", s8sum);
	fprintf(stderr, "wave equation:           %g\n", s9sum);
	}
	fprintf(stderr, "\n");
	return;
}

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

/* z component of the Poyntingvector, integrated along the x-axis / the interior layers, 
   normalization constant in case of a purely evanescent mode */
double cMode::power() const
{
	Complex s = CC0;
	for(int l=0; l<=N-1; ++l)
	{
		double x0, x1;
		x0 = piece[l].bif;	
		x1 = piece[l].tif;	
		if(l == 0 && bdt_b == OPEN)
		{
			if(piece[0].kappa.im < 0) x0 = bp_b;
			else                      x0 = piece[0].tif;
		}
		if(l == N-1 && bdt_t == OPEN)
		{
			if(piece[N-1].kappa.im < 0) x1 = bp_t;
			else                        x1 = piece[N-1].bif;
		}
		Complex pp = piece[l].intbfq(x0, x1);
		if(pol == TE) s += pp;
		else          s += pp/piece[l].eps;
	}
	if(fabs(gamma.re) < bqTol(0.0)) 
	{
		s = gamma*s;
		if(pol == TE) return 0.5*invommu0*s.im;
		return 0.5*invomep0*s.im;
	}
	s = gamma*s;
	if(pol == TE) return 0.5*invommu0*s.re;
	return 0.5*invomep0*s.re;
}

/* integration of a fieldproduct cp1^* cp2 over the interval [i.x0, i.x1] */
Complex cMode::integrate(Fcomp cp1, Fcomp cp2, Interval i) const
{
	return integrate(cp1, *this, cp2, i);
}

/* integrate a fieldproduct of two modes over Interval i,
   use the complex conjugate of *this' contribution */
Complex cMode::integrate(Fcomp cp, const cMode& m, Fcomp cpm, Interval i) const
{
	if(((int) cp ) > ((int) HZ)) cplxwgerror("cMode: integrate, cp");
	if(((int) cpm) > ((int) HZ)) cplxwgerror("cMode: integrate, cpm");
	if((pol == TE) && (cp == EX || cp == EZ || cp == HY)) return CC0;
	if((pol == TM) && (cp == HX || cp == HZ || cp == EY)) return CC0;
	if((m.pol == TE) && (cpm == EX || cpm == EZ || cpm == HY)) return CC0;
	if((m.pol == TM) && (cpm == HX || cpm == HZ || cpm == EY)) return CC0;
	
	double x, xt, xt1, xt2;
	Complex s, z;
	int l, lm;
	Interval iv;
	if(i.x0 > i.x1) { x=i.x1; i.x1=i.x0; i.x0=x; }
	if(  bdt_b != OPEN &&   bp_b > i.x0) i.x0 =   bp_b;
	if(  bdt_t != OPEN &&   bp_t < i.x1) i.x1 =   bp_t;
	if(m.bdt_b != OPEN && m.bp_b > i.x0) i.x0 = m.bp_b;
	if(m.bdt_t != OPEN && m.bp_t < i.x1) i.x1 = m.bp_t;
	s = CC0;
	x = i.x0;
	while(fabs(i.x1-x) > S_HDIST)
	{
		l   =   wg.layeridx(x+S_HDIST/2.0);
		lm  = m.wg.layeridx(x+S_HDIST/2.0);
		xt1 =   wg.layer(l).x1;
		xt2 = m.wg.layer(lm).x1;
		xt = xt1; 
		if( xt2<xt) xt = xt2;
		if(i.x1<xt) xt = i.x1;
		iv.x0 = x; iv.x1 = xt;
	
		int fod = 0, fodm = 0;
		Complex f = CC0, fm = CC0;
		Complex e;
		if(pol == TE)
		{
			switch(cp)
			{
			case EY: 
				fod = 0; f = CC1;
				break;
			case HX: 
				fod = 0; f = -gamma*invommu0;
				break;
			case HZ: 
				fod = 1; f = CCI*invommu0;
				break;
			default: cplxwgerror("cMode: integrate, TE, cp (II)");
				break;
			}
		}
		else 
		{
			e = wg.eps(0.5*(iv.x0+iv.x1));
			switch(cp)
			{
			case HY: 
				fod = 0; f = CC1;
				break;
			case EX: 
				fod = 0; f = gamma*invomep0/e;
				break;
			case EZ: 
				fod = 1; f = CCmI*invomep0/e;
				break;
			default: cplxwgerror("cMode: integrate, TM, cp (II)");
				break;
			}
		}
		if(m.pol == TE)
		{
			switch(cpm)
			{
			case EY: 
				fodm = 0; fm = CC1;
				break;
			case HX: 
				fodm = 0; fm = -m.gamma*m.invommu0;
				break;
			case HZ: 
				fodm = 1; fm = CCI*m.invommu0;
				break;
			default: cplxwgerror("cMode: integrate, TE, cpm (II)");
				break;
			}
		}
		else 
		{
			e = m.wg.eps(0.5*(iv.x0+iv.x1));
			switch(cpm)
			{
			case HY: 
				fodm = 0; fm = CC1;
				break;
			case EX: 
				fodm = 0; fm = m.gamma*m.invomep0/e;
				break;
			case EZ: 
				fodm = 1; fm = CCmI*m.invomep0/e;
				break;
			default: cplxwgerror("cMode: integrate, TM, cpm (II)");
				break;
			}
		}
		z = CC1;
		if(fod == 0)
		{
			if(fodm == 0) z = piece[l].intff(m.piece[lm], iv);
			else          z = piece[l].intfd(m.piece[lm], iv);
		}
		else 
		{
			if(fodm == 0) z = piece[l].intdf(m.piece[lm], iv);
			else          z = piece[l].intdd(m.piece[lm], iv);
		}
		s += z*f.conj()*fm;
		x = xt;
	}
	return s;
}

/* determine, whether the absolute squared field can be integrated over region r
   r = 'x': the entire x-axis,
   r = 'b': lowermost segment of the waveguide,
   r = 't': uppermost segment of the waveguide */
bool cMode::isintegrable(char r) const
{
	if(r == 'x') 
	{
		if(isintegrable('b') && isintegrable('t')) return true;
		return false;
	}
	if(r == 'b')
	{
		if(bdt_b == LIMD || bdt_b == LIMN) return true;
		int l = 0;
		if(piece[l].A.sqabs() > MININTCONTAMP) return false;
		if(piece[l].kappa.im < 0.0) return true;
		return false;
	}
	if(r == 't')
	{
		if(bdt_t == LIMD || bdt_t == LIMN) return true;
		int l = wg.nx+1;
		if(piece[l].B.sqabs() > MININTCONTAMP) return false;
		if(piece[l].kappa.im < 0.0) return true;
		return false;
	}
	return false;
}

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

/* principal field component, set maximum absolute level */
void cMode::setmaxF(double &ph)
{
	Interval i;
	double m, p;
	maxF = 0.0;
	for(int l=0; l<=N-1; ++l) 
	{
		i = wg.layer(l);
		if(l == 0       && bdt_b != OPEN) i.x0 = bp_b; 
		if(l == wg.nx+1 && bdt_t != OPEN) i.x1 = bp_t; 
		m = piece[l].maxf(i, p);
		if(m > maxF)
		{
			maxF = m;
			ph = p;
		}
	}
	return;
}

/* normalize mode to power() == 1, if possible */
void cMode::normalize()
{
	normalize(1.0);
}
/* normalize mode to power() == p, if possible */
void cMode::normalize(double p)
{
	double fac = fabs(power());
	if(fac < 1.0e-10)
	{ 
		fprintf(stderr, "cMode::normalize: fac = %g\n", fac);
		fac = 1.0e-10;
	}
	fac = sqrt(p)/sqrt(fac);
	for(int i=0; i<=N-1; ++i) piece[i].A *= fac;
	for(int i=0; i<=N-1; ++i) piece[i].B *= fac;
	double ph;
	setmaxF(ph);
	Complex efac = expmi(ph);
	for(int i=0; i<=N-1; ++i) piece[i].A *= efac;
	for(int i=0; i<=N-1; ++i) piece[i].B *= efac;
	return;
}

/* field profile -> .xyf data file 
	cp: EX - HZ
	foa: MOD, ORG, SQR, REP, IMP  
	disp: output region on the x-axis 
	ext0, ext1: filename id characters
	np: number of output points */
void cMode::writeprofile(Fcomp cp, Afo foa, Interval disp, int np,
                        char ext0, char ext1) const
{
	if(np <= 0) cplxwgerror("cMode: writeprofile: np <= 0");
	char name[13] = "t______";
	name[1] = polchr(pol);
	name[2] = afochr(foa);
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	double x0 = disp.x0;
	double x1 = disp.x1;


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

	Dvector f(np);
	Dvector x(np);
	if(np == 1) f(0) = realize(field(cp, FORW, 0.5*(x0+x1)), foa); 
	else 
	{
		double dx = (x1-x0)/((double) (np-1));
		for(int j=0; j<=np-1; ++j)  
		{
			x(j) = x0+((double) j)*dx;
			f(j) = realize(field(cp, FORW, x(j)), foa);
		}
	}
	toxyf(name, x, f);
	return;
}
/* as before, with default parameters */
void cMode::writeprofile(char ext0, char ext1) const
{
	Interval disp;
	if(bdt_b == OPEN) disp.x0 = wg.hx(0)-2.0;
	else              disp.x0 = bp_b;
	if(bdt_t == OPEN) disp.x1 = wg.hx(wg.nx)+2.0;
	else              disp.x1 = bp_t;
	writeprofile(principalcomp(pol), REP, disp, 1000, ext0, ext1);
	writeprofile(principalcomp(pol), IMP, disp, 1000, ext0, ext1);
	return;
}

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

/* field profile -> MATLAB .m file  
   cp: EX - HZ
   foa: MOD, ORG, SQR, REP, IMP
   disp: output region 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 command, to be included into *L.m */
#define SLHDISTF 1.0e-8
void cMode::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;
	Complex epsold, epsnew;
	double minf, maxf;
	int nsec, nbd;
	char name[13] = "t_______.m";
	name[1] = polchr(pol);
	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;

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

	switch(foa)
	{
	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: case ORG:
		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(wg.nx+2, np+1);
	Dmatrix pos(wg.nx+2, np+1);
	Ivector nump(wg.nx+2);
	Dvector of(np+1);
	Dvector op(np+1);
	Dvector bd(wg.nx+1);
	Dvector ri(wg.nx+2);
	nsec = 0;
	nbd  = 0;

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

	dat = fopen(name, "w+");
	mlout_title(dat, name, "Mode profile");
	name[8] = 0;
	Dmatrix fld = realize(cfld, foa);
	minf = -1.0;
	maxf = 1.0;
	if(pltype == 'L')
	{
		switch(foa)
		{
			case MOD: 
				maxf = cfld.abs().max();
				minf = 0.0;
				if(minf >= maxf) { minf = 0.0; maxf = 1.0; }
				break;
			case SQR:
				maxf = cfld.sqabs().max();
				minf = 0.0;
				if(minf >= maxf) { minf = 0.0; maxf = 1.0; }
				break;
			case ORG: case REP: case IMP:
				maxf = cfld.abs().max();
				minf = -maxf;
				if(minf >= maxf) { minf = -1.0; maxf = 1.0; }
	 			break;
		}
		mlout_geo(dat, wg, 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') 
	{
		mlout_Lsecplot(name, dat, cp, foa, nbd, nsec); 
		mlout_print(dat, name, 'e');
	}
	else mlout_Vsecplot(dat, cp, nsec, ext0, ext1); 
	fclose(dat);
	return;
}

/* as before, with default values */
void cMode::plot(char ext0, char ext1) const
{
	Interval disp;
	if(bdt_b == OPEN) disp.x0 = wg.hx(0)-2.0;
	else              disp.x0 = bp_b;
	if(bdt_t == OPEN) disp.x1 = wg.hx(wg.nx)+2.0;
	else              disp.x1 = bp_t;
	plot(principalcomp(pol), ORG, disp, 1000, ext0, ext1, 'L');
	plot(principalcomp(pol), ORG, disp, 1000, ext0, ext1, 'V');
	return;
}

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

/* compare the present object to m, return 1 if equal, otherwise return 0; 
   both modes are assumed to be normaliz()-ed eigenfunctions, incl. phase adjustment */
#define COMPTOL_BETA    1.0e-7
#define COMPTOL_BP      1.0e-7
#define COMPTOL_FLD     1.0e-7
#define COMPTOL_KAPPA   1.0e-10
int cMode::equal(const cMode& m) const
{
	if(this == &m) return 1;
	cWaveguide wg0 =   wg; wg0.optimize();
	cWaveguide wg1 = m.wg; wg1.optimize();
	if(wg0.equal(wg1) != 1) return 0;
	if(pol != m.pol) return 0;
	if((gamma-m.gamma).abs() > COMPTOL_BETA) return 0; 
	if(fabs(beta-m.beta) > COMPTOL_BETA) return 0; 
	if(fabs(alpha-m.alpha) > COMPTOL_BETA) return 0; 
	if(bdt_b != m.bdt_b) return 0;
	if(bdt_t != m.bdt_t) return 0;
	if(fabs(bp_b-m.bp_b) > COMPTOL_BP) return 0; 
	if(fabs(bp_t-m.bp_t) > COMPTOL_BP) return 0; 
	if(special != m.special) return 0;
	if(piece[0].bdktype != m.piece[0].bdktype)
	{
		if((piece[0].kappa-m.piece[0].kappa).abs() > COMPTOL_KAPPA) return 0;
	} 
	if(piece[wg.nx+1].bdktype != m.piece[m.wg.nx+1].bdktype) 
	{
		if((piece[wg.nx+1].kappa-m.piece[m.wg.nx+1].kappa).abs() > COMPTOL_KAPPA) return 0;
	} 
	if(fabs(maxF-m.maxF)/(0.5*(maxF+m.maxF)) > COMPTOL_FLD) return 0;
	for(int i=0; i<=wg.nx; ++i)
		if(fabs(field(wg.hx(i)).abs()/maxF - m.field(wg.hx(i)).abs()/m.maxF) > COMPTOL_FLD) return 0;
	return 1;
}

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

/* output to file */
void cMode::write(FILE *dat)
{
	comment("cMode", dat);	
	wg.write(dat);
	comment("pol", dat); fputint(pol, dat);
	comment("k0", dat); fputdouble(k0, dat);
	comment("1/(om mu0)", dat); fputdouble(invommu0, dat);
	comment("1/(om ep0)", dat); fputdouble(invomep0, dat);
	comment("gammaq", dat); gammaq.write(dat);
	comment("gamma", dat); gamma.write(dat);
	comment("neff", dat); neff.write(dat);
	comment("beta", dat); fputdouble(beta, dat);
	comment("alpha", dat); fputdouble(alpha, dat);
	comment("ord", dat); fputint(ord, dat);
	comment("bdt_b", dat); fputint(bdt_b, dat);
	comment("bdt_t", dat); fputint(bdt_t, dat);
	comment("bp_b", dat); fputdouble(bp_b, dat);
	comment("bp_t", dat); fputdouble(bp_t, dat);
	comment("maxF", dat); fputdouble(maxF, dat);
	comment("N", dat); fputint(N, dat);
	for(int l=0; l<=N-1; ++l) piece[l].write(dat);
	comment("ids", dat); 
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "%s\n", ids);	
	}
	else 
	{
		for(int i=0; i<=7; ++i) fputchar(ids[i], dat);
	}
	comment("btids", dat); 
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "%s\n", btids);	
	}
	else 
	{
		for(int i=0; i<=7; ++i) fputchar(btids[i], dat);
	}
	comment("special", dat); fputint(special, dat);
	comment("orgM", dat); orgM.write(dat);
	comment("cst", dat); fputint(cst, dat);
	return;
}

/* input from file */
void cMode::read(FILE *dat)
{
	wg.read(dat);
	pol = Polarization(fgetint(dat));
	k0 = fgetdouble(dat);
	invommu0 = fgetdouble(dat);
	invomep0 = fgetdouble(dat);
	gammaq.read(dat);
	gamma.read(dat);
	neff.read(dat);
	beta = fgetdouble(dat);
	alpha = fgetdouble(dat);
	ord = fgetint(dat);
	bdt_b = Boundary_type(fgetint(dat));
	bdt_t = Boundary_type(fgetint(dat));
	bp_b = fgetdouble(dat);
	bp_t = fgetdouble(dat);
	maxF = fgetdouble(dat);
	N = fgetint(dat);
	if(piece != NULL) delete[] piece;
	piece = new cSmPiece[N];
	for(int l=0; l<=N-1; ++l) piece[l].read(dat);
	for(int i=0; i<=7; ++i) ids[i] = fgetchar(dat);
	for(int i=0; i<=7; ++i) btids[i] = fgetchar(dat);
	special = fgetint(dat);
	orgM.read(dat);
	cst = fgetint(dat);
	return;
}

/* output to default file */
void cMode::write(char ext0, char ext1)
{
	char name[10] = "t___.smo";
	FILE *dat;
	name[1] = polchr(pol);
	name[2] = ext0;
	name[3] = ext1;
	fprintf(stderr, "Neff: "); neff.nice(stderr); fprintf(stderr, " >> %s\n", name); 
	dat = fopen(name, "w+");
	write(dat);
	fclose(dat);
	return;
}

/* input from default file */
void cMode::read(Polarization p, char ext0, char ext1)
{
	char name[10] = "t___.smo";
	FILE *dat;
	name[1] = polchr(p);
	name[2] = ext0;
	name[3] = ext1;
	fprintf(stderr, "<< %s : ", name);
	dat = fopen(name, "r");
	read(dat);
	fclose(dat);
	fprintf(stderr, "Neff: "); neff.nice(stderr); 
	return;
}

/* translate mode by dx */
void cMode::translate(double dx)
{ 
	wg.translate(dx);
	for(int l=0; l<=N-1; ++l) piece[l].translate(dx);	
	bp_b += dx;
	bp_t += dx;
	return;
}

/* - ModeArray: arrays of complex modes --------------------------------- */

/* initialize */
cModeArray::cModeArray()
          : num(0)
{
	mvec = NULL;
}

/* destroy */
cModeArray::~cModeArray()
{
	if(mvec != NULL) delete[] mvec; 
	mvec = NULL;
	num = 0; 
}

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

/* delete all Mode entries */
void cModeArray::clear()
{
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0; 
}

/* input from FILE dat */
void cModeArray::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) cplxwgerror("cModeArray: read");
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0; 
	num = fgetint(dat);
	if(num <= 0) { num = 0; return; }
	mvec = new cMode[num];
	for(int i=0; i<=num-1; ++i)
	{
		mvec[i].read(dat); 
	}
	return;
}

/* output to FILE dat */
void cModeArray::write(FILE *dat)
{
	comment("cModeArray", dat);
	comment("num", dat); fputint(num, dat);
	for(int i=0; i<=num-1; ++i)
	{
		comment("*** -------- ***", dat);
		mvec[i].write(dat); 
	}
	return;
}

/* input from default file */
void cModeArray::read(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	int i;
	name[0] = 's';
	name[1] = 'l';
	name[2] = 'a';
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 's';
	name[7] = 'm';
	name[8] = 'a';
	name[9] = 0;

	fprintf(stderr, "<< %s ", name);
	
	dat = fopen(name, "r");
	read(dat);
	fclose(dat);

	fprintf(stderr, "#%d\n", num);
	for(i=0; i<=num-1; ++i)
	{
		fprintf(stderr, "   (%d) %s, Neff: \n", i, mvec[i].ids); mvec[i].neff.nice(stderr); fprintf(stderr, "\n"); 
	}
	return;
}

/* output to default file */
void cModeArray::write(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	int i;
	name[0] = 's';
	name[1] = 'l';
	name[2] = 'a';
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 's';
	name[7] = 'm';
	name[8] = 'a';
	name[9] = 0;

	fprintf(stderr, ">> %s #%d\n", name, num);
	for(i=0; i<=num-1; ++i)
	{
		fprintf(stderr, "   (%d) %s, Neff: \n", i, mvec[i].ids); mvec[i].neff.nice(stderr); fprintf(stderr, "\n"); 
	}
	dat = fopen(name, "w+");
	write(dat);
	fclose(dat);

	return;
}

/* subscripting 
cMode& cModeArray::operator() (int i)  
{ 
	if(i<0 || i>=num) cplxwgerror("cModeArray: () out of range");
	return mvec[i]; 
}
cMode cModeArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) cplxwgerror("cModeArray: () out of range");
	return mvec[i]; 
}
*/

/* add a mode */
void cModeArray::add(cMode m)
{
	cMode *avec;
	avec = new cMode[num+1];
	cMode* ap = avec;
	cMode* 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 cModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		cplxwgerror("cModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	cMode *avec;
	avec = new cMode[num-1];
	cMode* ap = avec;
	cMode* 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 cModeArray ma */ 
void cModeArray::merge(cModeArray& ma)
{
	cMode *avec;
	avec = new cMode[num+ma.num];
	cMode* ap = avec;
	cMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	cMode* 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 cModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	int mina;
	cMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		maxi = j; 
		maxb = mvec[j].beta;
		mina = mvec[j].alpha;
		for(k=j+1; k<=num-1; ++k)
		{
			if(maxb<mvec[k].beta-(1.0e-13) ||
			   (fabs(maxb-mvec[k].beta)<(1.0e-13)
			    && mina > mvec[k].alpha))
			{
				maxb = mvec[k].beta;
				mina = mvec[k].alpha;
				maxi = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t;
	}
	return;
}

/* remove doublets from the array */
void cModeArray::prune()
{
	sort();
	int om = 0;
	while(om <= num-2)
	{
		if(mvec[om].equal(mvec[om+1]) == 1) 
		{
			if(mvec[om].cst < mvec[om+1].cst)
			     remove(om+1);
			else if(mvec[om].cst > mvec[om+1].cst)
			     remove(om);
			else 
			{
				// if(mvec[om].checkcontinuity() < mvec[om+1].checkcontinuity())
				if(fabs(mvec[om].beta-mvec[om].orgM.beta) < fabs(mvec[om+1].beta-mvec[om+1].orgM.beta))
				     remove(om+1);
				else remove(om);
			}
		}
		else ++om;
	}
	sort();	
	return;
}

/* compare the present object to ma, return 1 if equal, otherwise return 0 */
int cModeArray::equal(cModeArray& ma)
{
	if(this == &ma) return 1;
	if(num != ma.num) return 0;
	for(int i=0; i<=num-1; ++i) if(mvec[i].equal(ma(i)) != 1) return 0; 
	return 1;
}

/* ................................................................... */

/* 
 * mode interference evaluation and visualization 
 * all modes are assumed to belong to the same waveguide ! 
 */

/* bidirectional field superposition at point x, z 
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative)
   cp:   EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex cModeArray::field(const Cvector& famp, const Cvector& bamp, 
                         Fcomp cp, double x, double z) const
{
	return field(famp, bamp, cp, x, z, 0.0, 0.0);
}
	
/* bidirectional field superposition at point x, z in [z0, z1],
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1
   cp:   EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex cModeArray::field(const Cvector& famp, const Cvector& bamp, 
                         Fcomp cp, double x, double z,
                         double z0, double z1) const
{
	Complex s, ef, eb, f, b;
	Complex ex, ey, ez, hx, hy, hz;
	if(famp.nel < num) cplxwgerror("cModeArray: field: array length mismatch");
	if(bamp.nel < num) cplxwgerror("cModeArray: field: array length mismatch");
	s = CC0;
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ: 
		for(int i=0; i<=num-1; ++i)
		{
			f = mvec[i].field(cp, FORW, x, z-z0);
			if(fabs(famp(i).abs()) > NEGLECT_FLD_CONTRIB 
			   && std::isfinite(f.re) && std::isfinite(f.im)) 
					s += famp(i)*f;  
			b = mvec[i].field(cp, BACK, x, z-z1);
			if(fabs(bamp(i).abs()) > NEGLECT_FLD_CONTRIB 
			   && std::isfinite(b.re) && std::isfinite(b.im)) 
					s += bamp(i)*b;  
		}
		return s;
	case SX:
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		s = ey*hz.conj()-ez*hy.conj();
		return Complex(0.5*s.re, 0.0);
	case SY:
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		s = ez*hx.conj()-ex*hz.conj();
		return Complex(0.5*s.re, 0.0);
	case SZ:
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		s = ex*hy.conj()-ey*hx.conj();
		return Complex(0.5*s.re, 0.0);
	default: cplxwgerror("cModeArray: field: mysterious cp");
	}
	return CC0;
}

/* bidirectional field superposition, local field strength at x, z
     id: one of mE, mH, qE, qH, mW
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative) */
double cModeArray::lfs(const Cvector& famp, const Cvector& bamp, 
                      FSid id, double x, double z) const
{
	Complex ex, ey, ez, hx, hy, hz;
	double eps_x;
	switch(id)
	{
	case mE: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
	case mH: 
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
	case qE: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
	case qH: 
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
	case mW: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		eps_x = mvec[num-1].wg.eps(x).re;
		return 0.25*(      MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			    +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default: 
		cplxwgerror("cModeArray: lfs(), FSid?");
		break;
	}
	return 0.0;
}

/* bidirectional field superposition, local field strength at x, z
   at point x, z in [z0, z1],   
     id: one of mE, mH, qE, qH, mW
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1 */
double cModeArray::lfs(const Cvector& famp, const Cvector& bamp, 
                      FSid id, double x, double z, double z0, double z1) const
{
	Complex ex, ey, ez, hx, hy, hz;
	double eps_x;
	switch(id)
	{
	case mE: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
	case mH: 
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
	case qE: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
	case qH: 
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
	case mW: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		eps_x = mvec[num-1].wg.eps(x).re;
		return 0.25*(      MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			    +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default: 
		cplxwgerror("cModeArray: lfs(z0, z1), FSid?");
		break;
	}
	return 0.0;
}

/* bidirectional field superposition, local energy density at point x, z
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative) */
double cModeArray::endens(const Cvector& famp, const Cvector& bamp, 
                         double x, double z) const
{
	return lfs(famp, bamp, mW, x, z);
}

/* bidirectional field superposition, 
   local energy density at point x, z in [z0, z1],
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1 */
double cModeArray::endens(const Cvector& famp, const Cvector& bamp, 
                         double x, double z, double z0, double z1) const
{
	return lfs(famp, bamp, mW, x, z, z0, z1);
}

/* the six electric and magnetic components at x, z */
void cModeArray::emfield(const Cvector& famp, const Cvector& bamp,
                        double x, double z, double z0, double z1,
                        Complex& ex, Complex& ey, Complex& ez,
                        Complex& hx, Complex& hy, Complex& hz) const
{
	Complex s, ef, eb, f, b;
	if(famp.nel < num) cplxwgerror("cModeArray: field: array length mismatch");
	if(bamp.nel < num) cplxwgerror("cModeArray: field: array length mismatch");

	ex = ey = ez = hx = hy = hz = CC0;
	Complex af, ab;
	bool icf, icb;
	for(int i=0; i<=num-1; ++i)
	{
		af = famp(i);
		ab = bamp(i);
		ef = mvec[i].efac(FORW, z-z0);
		eb = mvec[i].efac(BACK, z-z1);

		icf = true;
		if(fabs(af.re) < NEGLECT_FLD_CONTRIB && 
		   fabs(af.im) < NEGLECT_FLD_CONTRIB) icf = false;
		if(!std::isfinite(ef.re)) icf = false;  
		if(!std::isfinite(ef.im)) icf = false;  
		icb = true;
		if(fabs(ab.re) < NEGLECT_FLD_CONTRIB && 
		   fabs(ab.im) < NEGLECT_FLD_CONTRIB) icb = false;
		if(!std::isfinite(eb.re)) icb = false;  
		if(!std::isfinite(eb.im)) icb = false;  

		if(icf==false && icb==false) return;

		af *= ef;
		ab *= eb;

		if(icf) ex += af*mvec[i].field(EX, FORW, x);  
		if(icb) ex += ab*mvec[i].field(EX, BACK, x);   

		if(icf) ey += af*mvec[i].field(EY, FORW, x);    
		if(icb) ey += ab*mvec[i].field(EY, BACK, x);     

		if(icf) ez += af*mvec[i].field(EZ, FORW, x);    
		if(icb) ez += ab*mvec[i].field(EZ, BACK, x);     

		if(icf) hx += af*mvec[i].field(HX, FORW, x);    
		if(icb) hx += ab*mvec[i].field(HX, BACK, x);     

		if(icf) hy += af*mvec[i].field(HY, FORW, x);    
		if(icb) hy += ab*mvec[i].field(HY, BACK, x);     

		if(icf) hz += af*mvec[i].field(HZ, FORW, x);    
		if(icb) hz += ab*mvec[i].field(HZ, BACK, x);     
	}
	return;
}

void cModeArray::emfield(const Cvector& famp, const Cvector& bamp,
                        double x, double z,
                        Complex& ex, Complex& ey, Complex& ez,
                        Complex& hx, Complex& hy, Complex& hz) const
{
	emfield(famp, bamp, x, z, 0.0, 0.0, ex, ey, ez, hx, hy, hz);
	return;
}

/* the optical field, discretized on a regular rectangular mesh,  
   npx x npz points on a display window dwx x dwz
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   foa: ORG, MOD, SQR, REP, IMP (ORG = REP) */
Cmatrix cModeArray::field(const Cvector& famp, const Cvector& bamp, 
                          Fcomp cp, Interval dwx, int npx, 
                                    Interval dwz, int npz) const
{
	if(npx <= 1 || npz <= 1) cplxwgerror("cModeArray: discretized field: npx, npz"); 
	Cmatrix fld(npx, npz);
	double dx = dwx.len()/((double)(npx-1));
	double dz = dwz.len()/((double)(npz-1));
	double x, z;

	z = dwz.x0;
	for(int j=0; j<=npz-1; ++j)
	{
		x = dwx.x0;
		for(int i=0; i<=npx-1; ++i) 
		{
			fld(i, j) = field(famp, bamp, cp, x, z);
			x += dx;
		}
		z += dz;
	}
	return fld;
}
Dmatrix cModeArray::field(const Cvector& famp, const Cvector& bamp, 
                          Fcomp cp, Afo foa, Interval dwx, int npx, 
                                             Interval dwz, int npz) const
{
	Cmatrix f = field(famp, bamp, cp, dwx, npx, dwz, npz);
	return realize(f, foa);
}

/* - Output to MATLAB .m files ---------------------------------------- */

/* mode interference plots, image of component cp
   famp: amplitudes of the forward propagating modes at z=0 
   bamp: amplitudes of the backward propagating modes at z=0
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
   foa: MOD, SQR, REP, IMP, ORG = REP
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void cModeArray::plot(const Cvector& famp, const Cvector& bamp, 
                     Fcomp cp, Afo foa,
                     double xbeg, double xend, double zbeg, double zend, 
                     int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "pl_____I.m";
	name[2] = afochr(foa);
        name[3] = fldchr(cp);
        name[4] = cpchr(cp);      
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "plot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	SegWgStruct s(-1);
	s(0) = mvec[0].wg.re();
	mlout_gengeoxz(dat, s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Dmatrix fld = field(famp, bamp, 
	                    cp, foa, Interval(xbeg, xend), npx, 
	                             Interval(zbeg, zend), npz);
	mlout_fld(dat, npx, npz, cp, fld);
	name[8] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lowlevcolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* image of the relative phase
   famp: amplitudes of the forward propagating modes at z=0 
   bamp: amplitudes of the backward propagating modes at z=0
   cp: one of EX, EY, EZ, HX, HY, HZ
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void cModeArray::phasemap(const Cvector& famp, const Cvector& bamp, 
                         Fcomp cp,
                         double xbeg, double xend, double zbeg, double zend,
                         int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "plp____I.m";
        name[3] = fldchr(cp);
        name[4] = cpchr(cp);      
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "phasemap([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Phasemap");
	SegWgStruct s(-1);
	s(0) = mvec[0].wg.re();
	mlout_gengeoxz(dat, s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Dmatrix fld = field(famp, bamp, 
	                    cp, Interval(xbeg, xend), npx, 
	                        Interval(zbeg, zend), npz).arg();
	mlout_fld(dat, npx, npz, cp, fld);
	name[8] = 0; 
	char desc[8];
	desc[0] = 'a'; desc[1] = 'r'; desc[2] = 'g'; desc[3] = ' ';
	desc[4] = fldchr(cp); desc[5] = cpchr(cp); desc[6] = 0;
	mlout_genimage(cp, MOD, name, desc, dat); 
	mlout_lincolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* image plot of the field "strength", id: one of mE, mH, qE, qH, mW  
   famp: amplitudes of the forward propagating modes at z=0 
   bamp: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void cModeArray::plot(FSid id, const Cvector& famp, const Cvector& bamp, 
                     double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	if(npx <= 1 || npz <= 1) 
		cplxwgerror("cModeArray: plot, lfs: npx, npz"); 
	FILE *dat;
	char name[13] = "pl____I.m";
	name[2] = mqchr(id);
	name[3] = idchr(id);
	name[4] = ext0;
	name[5] = ext1;
	fprintf(stderr, "plot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	SegWgStruct s(-1);
	s(0) = mvec[0].wg.re();
	mlout_gengeoxz(dat, s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Dmatrix fld(npx, npz);
	double x, dx, z, dz;
	dx = (xend-xbeg)/(npx-1);
	dz = (zend-zbeg)/(npz-1);
	z = zbeg;
	for(int j=0; j<=npz-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) 
		{
			fld(i, j) = lfs(famp, bamp, id, x, z);
			x += dx;
		}
		z += dz;
	}
	mlout_fld(dat, npx, npz, SZ, fld);
	name[7] = 0; 
	char desc[10];
	iddesc(id, desc);
	mlout_genimage(SZ, MOD, name, desc, dat); 
	mlout_lowlevcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* electromagnetic energy density image
   famp: amplitudes of the forward propagating modes at z=0 
   bamp: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void cModeArray::plot(const Cvector& famp, const Cvector& bamp, 
                     double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	return plot(mW, famp, bamp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
}

/* export full field data ("all") into a viewer m-file
   famp: amplitudes of the forward propagating modes at z=0 
   bamp: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters */
void cModeArray::viewer(const Cvector& famp, const Cvector& bamp,
                       double xbeg, double xend, double zbeg, double zend,
                       int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "fld__A.m";
	name[3] = ext0;
	name[4] = ext1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   xbeg, npx, xend, zbeg, npz, zend, name);
	dat = fopen(name, "w+");
	name[6] = 0; 
	Polarization pol = mvec[0].pol;
	double lambda = mvec[0].wg.lambda;
	for(int i=1; i<=num-1; ++i)
	{
		if(mvec[i].wg.lambda != lambda)
                                cplxwgerror("cModeArray: viewer: wavelengths mixed");
		if(mvec[i].pol != pol)
                                cplxwgerror("cModeArray: viewer: polarizations mixed");
	}
	mlout_viewertop(dat, name, pol, lambda, 0);
	SegWgStruct s(-1);
	s(0) = mvec[0].wg.re();
	mlout_gengeoxz(dat, s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix f;
	Fcomp cp;
	if(pol == TE)
	{
		cp = EX;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EZ;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HX;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HY;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HZ;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
	}
	else
	{
		cp = EX;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EY;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EZ;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HX;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HY;
		f = field(famp, bamp, 
		          cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HZ;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
	}
	Dmatrix n(npx, npz);
	double dx = (xend-xbeg)/((double)(npx-1));
	double x;
	for(int j=0; j<=npz-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) 
		{ 
			n(i, j) = s(0).n(x); 
			x += dx; 
		}
	}
	mlout_n(dat, npx, npz, n);
	mlout_fldviewer(dat, name);
	fclose(dat);
	return;
}

/* - Mode analysis of complex waveguides -------------------------------- */

/* ------------------ "tuning" parameters ------------------------------- */

/* tolerance for fixing the propagation constants */
double CPLXs_GqTol = 1.0e-13;

/* absolute tolerance */ 
double gqTol(Complex gq)
{
	double agq = gq.abs();
	if(agq < 10.0) return CPLXs_GqTol;
	return agq*CPLXs_GqTol;
}

/* minimum distance at which propagation constants are considered different */
double CPLXs_GqSep = 1.0e-8;

/* absolute separation limit */ 
double gqSep(Complex gq)
{
	double agq = gq.abs();
	if(agq < 10.0) return CPLXs_GqSep;
	return agq*CPLXs_GqSep;
}

/* tolerance for the maximum relative error in the continuity of the 
   basic field component for accepting a field as a valid mode */
double CPLXs_DiscErrTol = 1.0e-6;

/* --- mode analysis, helper functions, realized problem ----------- */

/* trap propagation constants */
#define MAXNUMSCF 100   // record this many suspicious configurations for the
                        // further propagation constant search       
int trapCR(ModeArray& ma, Mode& m, 
           double bqmin, double bqmax, double sf, int sm, 
           Dmatrix& susp, int& nsc, int lfsc, int verb)
{
	int newv;
	double step;
	double bq0, bq1, bq2, f0, f1, f2, mbq; 
	double err;

	if(bqmax <= bqmin) return 0;
	if(bqmax-bqmin < bqSep(bqmax)) return 0;

	newv = 0;
	double nst = 10.0*pow(10.0, sf);
	step = (bqmax-bqmin)/nst;
	if(step < bqTol(bqmax)) step = bqTol(bqmax); 
	if(step < bqTol(bqmin)) step = bqTol(bqmin); 

// fprintf(stderr, "trapCR(%.15g <--T%c-- %.15g <%.15g>, %d)\n", 
//                 bqmin, polCHR(m.pol), bqmax, step, sm);

	bq1 = bqmax;
	f1 = m.travres(bq1);
	bq0 = bqmax;
	f0 = f1; 
	for(double bq=bqmax-step; bq>=bqmin; bq-=step)
	{
		bq2 = bq;
		f2 = m.travres(bq2);
		if(f1*f2 <= 0.0) 
		{
			mbq = m.localize(bq2, bq1, bqTol(bq2));	
			m.polish(mbq);
			err = m.checkcontinuity();
// fprintf(stderr, "\n< %s >\n", m.ids);
// fprintf(stderr, "n_eff = %.16g\n", m.neff);
// fprintf(stderr, "betaq = %.16g\n", m.betaq);
// fprintf(stderr, "spec  = %d\n",    m.special);
// fprintf(stderr, "disc  = %.16g\n", m.checkcontinuity());
// fprintf(stderr, "dtol  = %.16g\n", CPLXs_DiscErrTol);
			if(err < CPLXs_DiscErrTol)
			{
				ma.add(m);
				if(verb == 2) fprintf(stderr, "X ");
				if(verb == 1) fprintf(stderr, ".");
				if(sm == 1) return 1;
			}
		}
		if(lfsc != 0)
		{
			if(f1 > 0.0 && f0 > f1 && f1 < f2) 
			{
				if(nsc < MAXNUMSCF-1)
				{
					susp(nsc, 0) = bq2;
					susp(nsc, 1) = bq0;
					++nsc;
				}
			}
			if(f1 < 0.0 && f0 < f1 && f1 > f2) 
			{
				if(nsc < MAXNUMSCF-1)
				{
					susp(nsc, 0) = bq2;
					susp(nsc, 1) = bq0;
					++nsc;
				}
			}
		}
		bq0 = bq1; f0 = f1;
		bq1 = bq2; f1 = f2;
	}
	return newv;
}

int trapCR(ModeArray& ma, Mode& m, 
          double bqmin, double bqmax, double sf, int sm, int lfsc, 
          int verb)
{
	int newv, nsc = 0;
	Dmatrix susp(MAXNUMSCF, 2);
	newv = trapCR(ma, m, bqmin, bqmax, sf, sm, susp, nsc, 1, verb);
	if(newv > 0) return newv;
	if(lfsc == 0) return newv;
	double minbq, maxbq;
	while(nsc >= 1)
	{
		minbq = susp(nsc-1, 0);
		maxbq = susp(nsc-1, 1);
		--nsc;
		newv += trapCR(ma, m, minbq, maxbq, sf, sm, susp, nsc, 0, verb);
		if(sm == 1 && newv > 0) return newv;
	}
	return newv;
}


/* --- mode analysis --------------------------------------------------- */

/* determine modes of the complex waveguide, 
   starting with real squared propagation constants between qqmin and qqmax, 
   for the problem specified by wg, pol, 
   boundaries of type bdtb at bpb, of type bdtt at bpt; 
   verb: 0: no control output, 1: progress, 2: all */
void findmodesC(cWaveguide wg, Polarization pol, 
                Boundary_type bdtb, double bpb,
                Boundary_type bdtt, double bpt,
                double gqmin, double gqmax, cModeArray& ma, int verb) 
{
	double k0 = val_k0(wg.lambda);
	double tgq, elim; 
	int cl, split, nsm;

// fprintf(stderr, "\nFM_I( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
// fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", bdtb, bpb, bdtt, bpt);
// fprintf(stderr, "      gqmin = %g, gqmax = %g )\n", gqmin, gqmax);
// verb = 2; 
// wg.write(stderr);
	
	if(gqmax < gqmin+bqSep(gqmax)) return;	
//	if(gqmax < gqmin+bqSep(gqmax)) cplxwgerror("findmodesC: gqmin, gqmax");	
	cl = wg.checksymmetry();
	split = 1;
	if(cl <= 0 || cl >= wg.nx+1) split = 0;
	if(bdtb != bdtt) split = 0;
	if(bdtb != OPEN && fabs((wg.hx(0)-bpb)-(bpt-wg.hx(wg.nx)))>COMPTOL_HX) 
		split = 0;
	if(split == 1)
	{
		cWaveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		wg0.optimize();
		wg1.optimize();
		cModeArray sm; 
		double cp = (wg.hx(cl-1)+wg.hx(cl))/2.0;
		if(verb == 1) fprintf(stderr, "s");
		if(verb == 2) fprintf(stderr, "sym ");
		sm.clear();
		findmodesC(wg0, pol, bdtb, bpb, LIMN, cp, 
		            gqmin, gqmax, sm, verb); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				;
			else sm(j).expand(); 
		}
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "a");
		if(verb == 2) fprintf(stderr, "asy ");
		sm.clear();
		findmodesC(wg0, pol, bdtb, bpb, LIMD, cp, 
		            gqmin, gqmax, sm, verb); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				sm(j).mirror(cp);
			else sm(j).expand(); 
		}
		ma.merge(sm);
		ma.sort();
		return;
	}	
	Waveguide wr = wg.re();
	elim = MAXREFIND*MAXREFIND; 
	cl = -1;
	for(int l=1; l<=wg.nx; ++l)
	{ 
		double te = wr.decoupepseff(l);
		if(te < elim) 
		{
			elim = te;
			cl = l;
		}
	}
	tgq = k0*k0*elim;
	if(tgq < gqmin) tgq = gqmin;
	if(gqmax > tgq+bqSep(tgq))
	{
		cWaveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		wg0.optimize();
		wg1.optimize();
		cModeArray sm;
		if(verb == 1) fprintf(stderr, "d-");
		if(verb == 2) fprintf(stderr, "d^(0, %d) ", cl);
		sm.clear();
		findmodesC(wg0, pol, bdtb, bpb, OPEN, AWAY, 
		           tgq, gqmax, sm, verb); 
		for(int j=0; j<=sm.num-1; ++j) { sm(j).special = 2; sm(j).setids(); }
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "d+");
		if(verb == 2) fprintf(stderr, "d^(%d, %d) ", cl, wg.nx+1);
		sm.clear();
		findmodesC(wg1, pol, OPEN, -AWAY, bdtt, bpt, 
		           tgq, gqmax, sm, verb); 
		for(int j=0; j<=sm.num-1; ++j) { sm(j).special = 2; sm(j).setids(); }
		ma.merge(sm);
		if(tgq > gqmin+bqSep(gqmin))
		{
			if(verb == 1) fprintf(stderr, "dr");
			if(verb == 2) fprintf(stderr, "d_(0, %d) ", wg.nx+1);
			sm.clear();
			findmodesC(wg, pol, bdtb, bpb, bdtt, bpt, 
				   gqmin, tgq, sm, verb); 
			ma.merge(sm);
		}
		ma.sort();
		return;
	}
	if(bdtb != OPEN)
	{
		elim = wr.lbdecepseff(bpb); 
		tgq = k0*k0*elim;
		if(tgq < gqmin) tgq = gqmin;
		if(gqmax > tgq+bqSep(tgq))
		{
			cModeArray sm;
			if(verb == 1) fprintf(stderr, "l^");
			if(verb == 2) fprintf(stderr, "l^ ");
			sm.clear();
			findmodesC(wg, pol, OPEN, -AWAY, bdtt, bpt, 
				   tgq, gqmax, sm, verb); 
			for(int j=0; j<=sm.num-1; ++j) { sm(j).special = 2; sm(j).setids(); }
			ma.merge(sm);
			if(gqmin < tgq-bqSep(tgq))
			{
				if(verb == 1) fprintf(stderr, "l_");
				if(verb == 2) fprintf(stderr, "l_ ");
				sm.clear();
				findmodesC(wg, pol, bdtb, bpb, bdtt, bpt, 
					    gqmin, tgq, sm, verb); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt != OPEN)
	{
		elim = wr.ubdecepseff(bpt); 
		tgq = k0*k0*elim;
		if(tgq < gqmin) tgq = gqmin;
		if(gqmax > tgq+bqSep(tgq))
		{
			cModeArray sm;
			if(verb == 1) fprintf(stderr, "u^");
			if(verb == 2) fprintf(stderr, "u^ ");
			sm.clear();
			findmodesC(wg, pol, bdtb, bpb, OPEN, AWAY, 
				   tgq, gqmax, sm, verb); 
			for(int j=0; j<=sm.num-1; ++j) { sm(j).special = 2; sm(j).setids(); }
			ma.merge(sm);
			if(gqmin < tgq-bqSep(tgq))
			{
				if(verb == 1) fprintf(stderr, "u_");
				if(verb == 2) fprintf(stderr, "u_ ");
				sm.clear();
				findmodesC(wg, pol, bdtb, bpb, bdtt, bpt, 
					   gqmin, tgq, sm, verb); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt==OPEN && bdtb==OPEN)
	{
		double emin = wg.defaultepseffmin();	
		tgq = k0*k0*emin;
		if(gqmin < tgq) gqmin = tgq;
		double emax = wg.defaultepseffmax();	
		tgq = k0*k0*emax;
		if(gqmax > tgq) gqmax = tgq;
		if(gqmax < gqmin+bqSep(gqmin)) return;
	}
	if(wg.isreal() && wg.defaultepseffmin() > bqSep(0.0)/k0/k0)
	{
		if(verb == 1) fprintf(stderr, ">RR");
		if(verb == 2) fprintf(stderr, ">RR ");
		ModeArray mar; 
		findmodes(wr, pol, bdtb, bpb, bdtt, bpt, gqmin, gqmax, mar, verb);
		for(int i=0; i<=mar.num-1; ++i)
		{
			cMode mm(mar(i));
			ma.add(mm);
		}
		ma.sort();
		return;
	}
	double remin = wr.eps(0);
	for(int l=1; l<=wr.nx+1; ++l)
	{
		if(wr.eps(l) < remin) remin = wr.eps(l);
	}
	if(remin > bqSep(0.0)/k0/k0)
	{
		if(verb == 1) fprintf(stderr, ">R");
		if(verb == 2) fprintf(stderr, ">R ");
		ModeArray mar; 
		findmodes(wr, pol, bdtb, bpb, bdtt, bpt, gqmin, gqmax, mar, verb);
		for(int i=0; i<=mar.num-1; ++i)
		{
			cMode mm(mar(i));
			cModeArray mat;
			mat = complexify(mm, wg, 0);
			if(mat.num >= 1) ma.merge(mat);
		}
		ma.sort();
		return;
	}
	Mode mr(wr, bdtb, bpb, bdtt, bpt, pol);
	if(gqmin >= gqmax) return;
//	fprintf(stderr, "\nFM_FCR( wg.nx = %d, pol = T%c\n", wr.nx, polCHR(pol));
//	fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//	                 bdtb, bpb, bdtt, bpt);
//	fprintf(stderr, "      gqmin = %.15g, gqmax = %g)\n", 
//			 gqmin, gqmax);
//	mr.travresinspect(gqmin+bqTol(gqmin), gqmax-bqTol(gqmax));

	ModeArray mar; 
	trapCR(mar, mr, gqmin, gqmax, 3.0, 0, 1, verb);

	mar.sort();
	int om = 0;
	while(om <= mar.num-2)
	{
		double ov = overlap(mar(om), FORW, mar(om+1), FORW).abs();
		if(ov>0.001)
		{ 
	// fprintf(stderr, " <%g|%g> = %g\n", mar(om).betaq, mar(om+1).betaq, ov);
			mar.remove(om+1);
			if(verb == 2) fprintf(stderr, "(r) ");
		}
		else 
		{
	// fprintf(stderr, " <%g|%g> = %g\n", mar(om).betaq, mar(om+1).betaq, ov);
			++om;
		}
	}
	mar.sort();

	for(int i=0; i<=mar.num-1; ++i)
	{
		cMode mm(mar(i));
		cModeArray mat;
		mat = complexify(mm, wg, 0);
		if(mat.num >= 1) ma.merge(mat);
	}
	ma.sort();
	if(verb >= 1) fprintf(stderr, "Ok.\n");
	return;
}

/* --- mode analysis, simpler driver routines -------------------------- */

/* guided mode analysis, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   ma: the modes found (output)
   quiet == 1: suppress log output */
int modeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, int quiet)         
{
	wg.consistency();
	wg.optimize();
	double k0 = val_k0(wg.lambda);

	if(quiet != 1)
	{
fprintf(stderr, "\n------------- cSlams - MA ------------------------------------------ '21 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
fprintf(stderr, "Tol: %g  ", CPLXs_GqTol);
fprintf(stderr, "Sep: %g\n", CPLXs_GqSep);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", wg.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=wg.nx; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    n:");
for(int j=0; j<=wg.nx+1; ++j) { fprintf(stderr, "  "); wg.n(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "  eps:");
for(int j=0; j<=wg.nx+1; ++j) { fprintf(stderr, "  "); wg.eps(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "lambda: %.10g  k0: %g\n", wg.lambda, k0);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	ma.clear();

	double gqmin = wg.defaultepseffmin()*k0*k0;
	double gqmax = wg.defaultepseffmax()*k0*k0;
	int verb = 0; if(quiet != 1) verb = 2;
	findmodesC(wg, pol, OPEN, -AWAY, OPEN, AWAY, gqmin, gqmax, ma, verb);
	ma.sort();
	ma.prune();
	if(quiet != 1)
	{
		if(ma.num<=0) fprintf(stderr, "No modes found.\n");
		else
fprintf(stderr, "\n----------------------------------------------------------------------------\n");
		for(int j=0; j<=ma.num-1; ++j)
		{
fprintf(stderr, "[%d] %s", j, ma(j).ids);
fprintf(stderr, ": n_eff = "); ma(j).neff.nice(stderr);
fprintf(stderr, ", gamma = "); ma(j).gamma.nice(stderr);
fprintf(stderr, "\n");
		}
fprintf(stderr, "\n");
	}
	return ma.num;
}

int modeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma)         
{
	return modeanalysis(wg, pol, ma, 0);
}


/* leaky mode analysis, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   ma: the modes found (output)
   epsout: permittivity applied to outer layers for initial bound mode computation
   quiet == 1: suppress log output */
int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, double epsout, int quiet)         
{
	wg.consistency();
	wg.optimize();
	double k0 = val_k0(wg.lambda);
	ma.clear();

	cWaveguide wg0 = wg;
	Waveguide w0 = wg0.re();
	wg0.unleak(epsout);

	if(quiet != 1)
	{
fprintf(stderr, "\n------------- cSlams - LMA ----------------------------------------- '21 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
fprintf(stderr, "Tol: %g  ", CPLXs_GqTol);
fprintf(stderr, "Sep: %g\n", CPLXs_GqSep);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "lambda: %.10g  k0: %g\n", wg.lambda, k0);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", wg.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=wg.nx; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    n:");
for(int j=0; j<=wg.nx+1; ++j) { fprintf(stderr, "  "); wg.n(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "  eps:");
for(int j=0; j<=wg.nx+1; ++j) { fprintf(stderr, "  "); wg.eps(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}

	double gqmin = wg.defaultepseffmin()*k0*k0;
	double gqmax = wg.defaultepseffmax()*k0*k0;
	int verb = 0; if(quiet != 1) verb = 2;
	findmodesC(wg, pol, OPEN, -AWAY, OPEN, AWAY, gqmin, gqmax, ma, verb);
	ma.sort();

	if(quiet != 1)
	{
		if(ma.num<=0) fprintf(stderr, "No modes found.\n");
		else
fprintf(stderr, "\n-- GMA ---------------------------------------------------------------------\n");
		for(int j=0; j<=ma.num-1; ++j)
		{
fprintf(stderr, "%s", ma(j).ids);
fprintf(stderr, ": n_eff = "); ma(j).neff.nice(stderr);
fprintf(stderr, ", gamma = "); ma(j).gamma.nice(stderr);
fprintf(stderr, "\n");
		}
fprintf(stderr, "\n");
	}

	if(wg.equal(wg0) == 0)
	{
		if(quiet != 1)
		{
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", wg0.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=wg0.nx; ++j) fprintf(stderr, "%6.10g  ", wg0.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    n:");
for(int j=0; j<=wg0.nx+1; ++j) { fprintf(stderr, "  "); wg0.n(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "  eps:");
for(int j=0; j<=wg0.nx+1; ++j) { fprintf(stderr, "  "); wg0.eps(j).nice(stderr); }
fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
		}
		cModeArray gma;
		gma.clear();

		gqmin = wg0.defaultepseffmin()*k0*k0;
		gqmax = wg0.defaultepseffmax()*k0*k0;
		verb = 0; if(quiet != 1) verb = 2;
		findmodesC(wg0, pol, OPEN, -AWAY, OPEN, AWAY, gqmin, gqmax, gma, verb);
		gma.sort();

		if(quiet != 1)
		{
			if(gma.num<=0) fprintf(stderr, "No modes found.\n");
			else
fprintf(stderr, "\n-- _GMA_ -------------------------------------------------------------------\n");
			for(int j=0; j<=gma.num-1; ++j)
			{
fprintf(stderr, "%s", gma(j).ids);
fprintf(stderr, ": n_eff = "); gma(j).neff.nice(stderr);
fprintf(stderr, ", gamma = "); gma(j).gamma.nice(stderr);
fprintf(stderr, "\n");
			}
fprintf(stderr, "\n");
		}
		
		for(int i=0; i<=gma.num-1; ++i)
		{
			cMode mm = gma(i);
			cModeArray tma = leakify(mm, wg, quiet);
			ma.merge(tma);
		}
	}
	ma.sort();
	ma.prune();
	if(quiet != 1)
	{
		if(ma.num<=0) fprintf(stderr, "No modes found.\n");
		else
fprintf(stderr, "\n-- LMA ---------------------------------------------------------------------\n");
		for(int j=0; j<=ma.num-1; ++j)
		{
fprintf(stderr, "[%d] %s", j, ma(j).ids);
fprintf(stderr, ": n_eff = "); ma(j).neff.nice(stderr);
fprintf(stderr, ", gamma = "); ma(j).gamma.nice(stderr);
fprintf(stderr, "\n");
		}
fprintf(stderr, "\n");
	}

	return ma.num;
}

int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, double epsout)         
{
	return leakymodeanalysis(wg, pol, ma, epsout, 0);
}

int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, int quiet)         
{
	wg.consistency();
	wg.optimize();

	Waveguide w0 = wg.re();
	double emin = AWAY;
	for(int l=0; l<=w0.nx+1; ++l)
	{
		if(w0.eps(l) < emin) emin = w0.eps(l);
	}
	if(w0.nx >= 2)
	{
		double emini = AWAY;
		for(int l=1; l<=w0.nx; ++l)
		{
			if(w0.eps(l) < emini) emini = w0.eps(l);
		}
		if(emini > emin) emin = emini;
	}
	return leakymodeanalysis(wg, pol, ma, emin, 0);
}

int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma)         
{
	return leakymodeanalysis(wg, pol, ma, 0);
}