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

/*
 * slamode.cpp
 * Modes of multilayer slab waveguides, mode overlaps.
 */

#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"

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

/* field evaluation: 
   neglect contributions with amplitudes below NEGLECT_FLD_CONTRIB */
double NEGLECT_FLD_CONTRIB = 1.0e-300; 

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

/* initiate: set polarization, interfaces, permittivity */
void SmPiece::init(Polarization p, Interval i, double ev, double lambda)
{
	pol = p;
	bif = i.x0; 
	tif = i.x1; 
	k0 = val_k0(lambda);
	k0nq = k0*k0*ev;
	inve = 1.0/ev;
	return;
}
/* ... VEIMS nonstandard effective TM profile */
void SmPiece::init(Interval i, double ev, double ecv, double lambda)
{
	pol = TM;
	bif = i.x0; 
	tif = i.x1; 
	k0 = val_k0(lambda);
	k0nq = k0*k0*ev;
	inve = 1.0/ecv;
	return;
}

#define GAMMAQSHIFT 1.0e-40

/* mode solver, transverse resonance evaluation: start */
void SmPiece::start_tre(double bq, double& phi, double &phis, 
                        Boundary_type bdt)
{
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = tif;
			A = 1.0; 
			B = 0.0;
			phi = A;
			phis = g*A;
			break;
		case OSC:
			slabmodeerror("start_tre: OPEN + OSC");
			break;
		}
		break;
	case LIMD:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			A = 1.0/(1.0-exp(-g*(tif-xb)));
			B = -A;
			phi = 1.0;
			phis = g*A*(1.0+exp(-g*(tif-xb)));
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			A = 1.0;
			B = 0;
			phi = sin(g*(tif-xa));
			phis = g*cos(g*(tif-xa));
			break;
		}
		break;
	case LIMN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			A = 1.0/(1.0+exp(-g*(tif-xb)));
			B = A; 
			phi = 1.0; 
			phis = g*A*(1.0-exp(-g*(tif-xb)));
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			B = 1.0;
			A = 0;
			phi = cos(g*(tif-xb));
			phis = -g*sin(g*(tif-xb));
			break;
		}
		break;
	} 
	if(pol == TM) phis *= inve;
	return;
	
}

/* mode solver, transverse resonance evaluation: continue */
void SmPiece::cont_tre(double bq, double& phi, double &phis)
{
	if(pol == TM) phis /= inve;
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
 	double dg, expdg, sindg, cosdg;
	switch(bhv)
	{
	case HYP:
		xa = tif; xb = bif;
		dg = g*(tif-bif);
		expdg = exp(-dg); 
		A = 0.5*(phi+phis/g)/expdg; 
		B = 0.5*(phi-phis/g); 
		phi  =   A+expdg*B;
		phis = g*A-g*expdg*B;
		break;
	case OSC:
		xa = xb = 0.5*(tif+bif);
		dg = g*(tif-xa);
		sindg = sin(dg); cosdg = cos(dg);
		A = -sindg*phi+cosdg*phis/g; 
		B =  cosdg*phi+sindg*phis/g; 
		phi  = sindg*A+cosdg*B;
		phis = g*cosdg*A-g*sindg*B;
		break;
	}
	if(pol == TM) phis *= inve;
	return;
}	

/* mode solver, transverse resonance evaluation: finish */
double SmPiece::finish_tre(double bq, double& phi, double &phis, 
                           Boundary_type bdt)
{
	if(pol == TM) phis /= inve;
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
	double trs = 1.0;
	double dg, sindg, cosdg;
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			xa = bif; xb = bif;	
			A = 0.5*(phi+phis/g); 
			B = 0.5*(phi-phis/g); 
			trs = A;
			break;
		case OSC:
			slabmodeerror("finish_tre: OPEN + OSC");
			break;
		}
		break;
	case LIMD: case LIMN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			dg = g*(tif-bif);
			A = 0.5*(phi+phis/g)/exp(-dg); 
			B = 0.5*(phi-phis/g); 
			if(bdt == LIMD) trs = A+B;
			else            trs = A-B;
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			dg = g*(tif-xa);
			sindg = sin(dg); cosdg = cos(dg);
			A = -sindg*phi+cosdg*phis/g; 
			B =  cosdg*phi+sindg*phis/g; 
			if(bdt == LIMD) trs = B;
			else            trs = A;
			break;
		}
		break;
	} 
	return trs;
}

/* set parameters in an upper segment for a modal solution */
void SmPiece::polish_ub(Boundary_type bdt)
{
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			A = 0.0;
			break;
		case OSC:
			slabmodeerror("polish_ub: OPEN + OSC");
			break;
		}
		break;
	case LIMD: case LIMN:
		switch(bhv)
		{
		case HYP:
			if(bdt == LIMD) A = -B;
			else            A =  B;
			break;
		case OSC:
			if(bdt == LIMD) B = 0.0;
			else            A = 0.0;
			break;
		}
		break;
	} 
	return;
}

/* mode solver, reverse transverse resonance evaluation: start */
void SmPiece::start_revtre(double bq, double& phi, double &phis, 
                           Boundary_type bdt)
{
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			xa = bif; xb = bif;
			A = 0.0; 
			B = 1.0;
			phi = B;
			phis = -g*B;
			break;
		case OSC:
			slabmodeerror("start_revtre: OPEN + OSC");
			break;
		}
		break;
	case LIMD:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			B = 1.0/(1.0-exp(-g*(tif-xb)));
			A = -B;
			phi = 1.0;
			phis = -g*B*(1.0+exp(-g*(tif-xb)));
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			A = 1.0;
			B = 0;
			phi = sin(g*(bif-xa));
			phis = g*cos(g*(bif-xa));
			break;
		}
		break;
	case LIMN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			B = 1.0/(1.0+exp(-g*(tif-xb)));
			A = B; 
			phi = 1.0; 
			phis = -g*B*(1.0-exp(-g*(tif-xb)));
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			B = 1.0;
			A = 0;
			phi = cos(g*(bif-xb));
			phis = -g*sin(g*(bif-xb));
			break;
		}
		break;
	} 
	if(pol == TM) phis *= inve;
	return;
	
}

/* mode solver, reverse transverse resonance evaluation: continue */
void SmPiece::cont_revtre(double bq, double& phi, double &phis)
{
	if(pol == TM) phis /= inve;
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
 	double dg, expdg, sindg, cosdg;
	switch(bhv)
	{
	case HYP:
		xa = tif; xb = bif;
		dg = g*(tif-bif);
		expdg = exp(-dg); 
		A = 0.5*(phi+phis/g); 
		B = 0.5*(phi-phis/g)/expdg; 
		phi  =   A*expdg+B;
		phis = g*A*expdg-g*B;
		break;
	case OSC:
		xa = xb = 0.5*(tif+bif);
		dg = g*(tif-xa);
		sindg = sin(dg); cosdg = cos(dg);
		A =  sindg*phi+cosdg*phis/g; 
		B =  cosdg*phi-sindg*phis/g; 
		phi  =  -sindg*A+cosdg*B;
		phis = g*cosdg*A+g*sindg*B;
		break;
	}
	if(pol == TM) phis *= inve;
	return;
}	

/* mode solver, reverse transverse resonance evaluation: finish */
double SmPiece::finish_revtre(double bq, double& phi, double &phis, 
                              Boundary_type bdt)
{
	if(pol == TM) phis /= inve;
	double gq;
	gq = bq - k0nq;
	if(fabs(gq) < GAMMAQSHIFT) gq = -GAMMAQSHIFT;
	if(gq > 0.0) bhv = HYP; 
	else         bhv = OSC; 
	g = sqrt(fabs(gq));
	double trs = 1.0;
	double dg, sindg, cosdg;
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = tif;	
			A = 0.5*(phi+phis/g); 
			B = 0.5*(phi-phis/g); 
			trs = B;
			break;
		case OSC:
			slabmodeerror("finish_revtre: OPEN + OSC");
			break;
		}
		break;
	case LIMD: case LIMN:
		switch(bhv)
		{
		case HYP:
			xa = tif; xb = bif;
			dg = g*(tif-bif);
			A = 0.5*(phi+phis/g); 
			B = 0.5*(phi-phis/g)/exp(-dg); 
			if(bdt == LIMD) trs = A+B;
			else            trs = A-B;
			break;
		case OSC:
			xa = xb = 0.5*(tif+bif);
			dg = g*(tif-xa);
			sindg = sin(dg); cosdg = cos(dg);
			A = sindg*phi+cosdg*phis/g; 
			B = cosdg*phi-sindg*phis/g; 
			if(bdt == LIMD) trs = B;
			else            trs = A;
			break;
		}
		break;
	} 
	return trs;
}

/* set parameters in a lower segment for a modal solution */
void SmPiece::polish_lb(Boundary_type bdt)
{
	switch(bdt)
	{
	case OPEN:
		switch(bhv)
		{
		case HYP:
			B = 0.0;
			break;
		case OSC:
			slabmodeerror("polish_lb: OPEN + OSC");
			break;
		}
		break;
	case LIMD: case LIMN:
		switch(bhv)
		{
		case HYP:
			if(bdt == LIMD) B = -A;
			else            B =  A;
			break;
		case OSC:
			if(bdt == LIMD) B = 0.0;
			else            A = 0.0;
			break;
		}
		break;
	} 
	return;
}

/* field evaluation within a specific layer: principal component */
double SmPiece::bfld(double x) const
{
//	return (bhv == HYP) ? A*exp(g*(x-xa))+B*exp(-g*(x-xb))
//                          : A*sin(g*(x-xa))+B*cos( g*(x-xb));
	double f = 0.0;
	if(bhv == HYP)
	{
		if(fabs(A)>NEGLECT_FLD_CONTRIB) f += A*exp( g*(x-xa));
		if(fabs(B)>NEGLECT_FLD_CONTRIB) f += B*exp(-g*(x-xb));
		return f;
	}
	if(fabs(A)>NEGLECT_FLD_CONTRIB) f += A*sin(g*(x-xa));
	if(fabs(B)>NEGLECT_FLD_CONTRIB) f += B*cos(g*(x-xb));
	return f;
}

/* field evaluation within a specific layer: 
   the derivative of the principal component with respect to x */
double SmPiece::bder(double x) const
{
//	return (bhv == HYP) ? g*(A*exp(g*(x-xa))-B*exp(-g*(x-xb)))
//	                    : g*(A*cos(g*(x-xa))-B*sin( g*(x-xb)));
	double f = 0.0;
	if(bhv == HYP)
	{
		if(fabs(A)>NEGLECT_FLD_CONTRIB) f += g*A*exp( g*(x-xa));
		if(fabs(B)>NEGLECT_FLD_CONTRIB) f -= g*B*exp(-g*(x-xb));
		return f;
	}
	if(fabs(A)>NEGLECT_FLD_CONTRIB) f += g*A*cos(g*(x-xa));
	if(fabs(B)>NEGLECT_FLD_CONTRIB) f -= g*B*sin(g*(x-xb));
	return f;
}

/* number of nodes in the principal field component in interval x0, x1 */
int SmPiece::nodes(double x0, double x1) const
{
	return (bhv == HYP) ? exp_nodes(A, B, g, xa, xb, x0, x1) 
	                    : harmon_nodes(A, B, g, xa, x0, x1); 
}

/* number of nodes in the derivative of the b. f. c. in Interval x0, x1 */
int SmPiece::dnodes(double x0, double x1) const
{
	return (bhv == HYP) ? exp_nodes(g*A, -g*B, g, xa, xb, x0, x1) 
	                    : harmon_nodes(-g*B, g*A, g, xa, x0, x1); 
}

/* maximum absolute level of the principal field component within a layer */
double SmPiece::maxf(Interval i) const
{
	return (bhv == HYP) ? exp_max(A, B, g, xa, xb, i.x0, i.x1) 
	                    : harmon_max(A, B, g, xa, i.x0, i.x1); 
}


/* integrate the squared principal field component over the entire layer */
double SmPiece::intbfq() const
{
	return (bhv == HYP) ? intfqHYP(A, B, g, xa, xb, bif, tif)
			    : intfqOSC(A, B, g, xa, bif, tif); 
}

/* integrate the squared principal field component over Interval iv */
double SmPiece::intbfq(double x0, double x1) const
{
	return (bhv == HYP) ? intfqHYP(A, B, g, xa, xb, x0, x1)
			    : intfqOSC(A, B, g, xa, x0, x1); 
}

/* integrate the present principal field component
   times the principal component of modepiece mp over interval iv */
double SmPiece::intff(const SmPiece& mp, const Interval& iv) const
{
	double s = 0.0;
	switch(bhv)
	{
	case HYP:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfEE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfEE(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEE(-g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return s;
			break;
		case OSC:
			s=A*(mp.A*intfES( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEC( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfES(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEC(-g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return s;
			break;
		}
		break;
	case OSC:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfSE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfSE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfCE( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfCE( g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return s;
			break;
		case OSC:
			s=A*(mp.A*intfSS( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfSC( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfCS( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfCC( g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return s;
			break;
		}
		break;
	}
	return 0.0; 
}

/* integrate the present principal field component times the derivative
   of the principal component of modepiece mp over interval iv */
double SmPiece::intfd(const SmPiece& mp, const Interval& iv) const
{
	double s = 0.0;
	switch(bhv)
	{
	case HYP:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfEE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfEE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfEE(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfEE(-g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return mp.g*s;
			break;
		case OSC:
			s=A*(mp.A*intfEC( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfES( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfEC(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfES(-g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return mp.g*s;
			break;
		}
		break;
	case OSC:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfSE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfSE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfCE( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfCE( g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return mp.g*s;
			break;
		case OSC:
			s=A*(mp.A*intfSC( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfSS( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 +B*(mp.A*intfCC( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfCS( g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return mp.g*s;
			break;
		}
		break;
	}
	return 0.0; 
}

/* integrate the derivative of the present principal field component
   times the principal component of modepiece mp over interval iv */
double SmPiece::intdf(const SmPiece& mp, const Interval& iv) const
{
	double s = 0.0;
	switch(bhv)
	{
	case HYP:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfEE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfEE(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEE(-g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return g*s;
			break;
		case OSC:
			s=A*(mp.A*intfES( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEC( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfES(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfEC(-g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return g*s;
			break;
		}
		break;
	case OSC:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfCE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfCE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfSE( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfSE( g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return g*s;
			break;
		case OSC:
			s=A*(mp.A*intfCS( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfCC( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfSS( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    +mp.B*intfSC( g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return g*s;
			break;
		}
		break;
	}
	return 0.0; 
}

/* integrate the derivative of the present principal field component times 
   the derivative of the principal component of modepiece mp over interval iv */
double SmPiece::intdd(const SmPiece& mp, const Interval& iv) const
{
	double s = 0.0;
	switch(bhv)
	{
	case HYP:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfEE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfEE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfEE(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfEE(-g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return g*mp.g*s;
			break;
		case OSC:
			s=A*(mp.A*intfEC( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfES( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfEC(-g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfES(-g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return g*mp.g*s;
			break;
		}
		break;
	case OSC:
		switch(mp.bhv)
		{
		case HYP:
			s=A*(mp.A*intfCE( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfCE( g, xa, -mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfSE( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfSE( g, xb, -mp.g, mp.xb, iv.x0, iv.x1));
			return g*mp.g*s;
			break;
		case OSC:
			s=A*(mp.A*intfCC( g, xa,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfCS( g, xa,  mp.g, mp.xb, iv.x0, iv.x1))
			 -B*(mp.A*intfSC( g, xb,  mp.g, mp.xa, iv.x0, iv.x1)
			    -mp.B*intfSS( g, xb,  mp.g, mp.xb, iv.x0, iv.x1));
			return g*mp.g*s;
			break;
		}
		break;
	}
	return 0.0; 
}

/* integrate the principal field component times a phase term over interval iv */
Complex SmPiece::crintfCiS(double cr, double xr, const Interval& iv) const
{
	double sr, si;
	switch(bhv)
	{
	case HYP:
		sr  = A*intfEC( g, xa, cr, xr, iv.x0, iv.x1);
		sr += B*intfEC(-g, xb, cr, xr, iv.x0, iv.x1);
		si  = A*intfES( g, xa, cr, xr, iv.x0, iv.x1);
		si += B*intfES(-g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr, -si);
		break;
	case OSC:
		sr  = A*intfSC( g, xa, cr, xr, iv.x0, iv.x1);
		sr += B*intfCC( g, xb, cr, xr, iv.x0, iv.x1);
		si  = A*intfSS( g, xa, cr, xr, iv.x0, iv.x1);
		si += B*intfCS( g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr, -si);
		break;
	}
	return CC0; 
}

/* integrate the derivative of the b. f. c. times a phase term over iv */
Complex SmPiece::crintdCiS(double cr, double xr, const Interval& iv) const
{
	double sr, si;
	switch(bhv)
	{
	case HYP:
		sr  = A*intfEC( g, xa, cr, xr, iv.x0, iv.x1);
		sr -= B*intfEC(-g, xb, cr, xr, iv.x0, iv.x1);
		si  = A*intfES( g, xa, cr, xr, iv.x0, iv.x1);
		si -= B*intfES(-g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr*g, -si*g);
		break;
	case OSC:
		sr  = A*intfCC( g, xa, cr, xr, iv.x0, iv.x1);
		sr -= B*intfSC( g, xb, cr, xr, iv.x0, iv.x1);
		si  = A*intfCS( g, xa, cr, xr, iv.x0, iv.x1);
		si -= B*intfSS( g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr*g, -si*g);
		break;
	}
	return CC0; 
}

/* integrate the principal field component times an exponential over interval iv */
Complex SmPiece::crintfExp(double cr, double xr, const Interval& iv) const
{
	double sr; 
	switch(bhv)
	{
	case HYP:
		sr  = A*intfEE( g, xa, cr, xr, iv.x0, iv.x1);
		sr += B*intfEE(-g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr, 0.0);
		break;
	case OSC:
		sr  = A*intfSE( g, xa, cr, xr, iv.x0, iv.x1);
		sr += B*intfCE( g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr, 0.0);
		break;
	}
	return CC0; 
}

/* integrate the derivative of the b. f. c. times an exponential over iv */
Complex SmPiece::crintdExp(double cr, double xr, const Interval& iv) const
{
	double sr; 
	switch(bhv)
	{
	case HYP:
		sr  = A*intfEE( g, xa, cr, xr, iv.x0, iv.x1);
		sr -= B*intfEE(-g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr*g, 0.0);
		break;
	case OSC:
		sr  = A*intfCE( g, xa, cr, xr, iv.x0, iv.x1);
		sr -= B*intfSE( g, xb, cr, xr, iv.x0, iv.x1);
		return Complex(sr*g, 0.0);
		break;
	}
	return CC0; 
}

/* integrate the principal field component times a complex exponential 
   (exp -i kr (x-xr) over interval iv */
Complex SmPiece::crintfcE(Complex kr, double xr, const Interval& iv) const
{
	Complex kf, k;
	double xf; 
	Complex s = CC0; 
	Complex c; 
	switch(bhv)
	{
	case HYP:
		c = Complex(A, 0.0); kf = Complex(0.0,  g); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(B, 0.0); kf = Complex(0.0, -g); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		return s;
		break;
	case OSC:
		c = Complex(0.0, -0.5*A); kf = Complex(-g, 0.0); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(0.0,  0.5*A); kf = Complex( g, 0.0); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex( 0.5*B, 0.0); kf = Complex(-g, 0.0); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex( 0.5*B, 0.0); kf = Complex( g, 0.0); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		return s;
		break;
	}
	return CC0; 
}

/* integrate the derivative of the pr. f. component times a complex exponential 
   (exp -i kr (x-xr) over interval iv */
Complex SmPiece::crintdcE(Complex kr, double xr, const Interval& iv) const
{
	Complex kf, k;
	double xf; 
	Complex s = CC0; 
	Complex c; 
	switch(bhv)
	{
	case HYP:
		c = Complex( g*A, 0.0); kf = Complex(0.0,  g); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(-g*B, 0.0); kf = Complex(0.0, -g); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		return s;
		break;
	case OSC:
		c = Complex(0.5*g*A,  0.0); kf = Complex(-g, 0.0); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(0.5*g*A,  0.0); kf = Complex( g, 0.0); xf = xa; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(0.0,  0.5*g*B); kf = Complex(-g, 0.0); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		c = Complex(0.0, -0.5*g*B); kf = Complex( g, 0.0); xf = xb; k = kf+kr; s += c*expi(kf*(xf-xr))*intcexp(k, xr, iv.x0, iv.x1);
		return s;
		break;
	}
	return CC0; 
}

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

/* output to file */
void SmPiece::write(FILE *dat) const
{
	comment("(( SmPiece))", 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); fputdouble(A, dat);
	comment("B", dat); fputdouble(B, dat);
	comment("g", dat); fputdouble(g, dat);
	comment("xa", dat); fputdouble(xa, dat);
	comment("xb", dat); fputdouble(xb, dat);
	comment("bhv", dat); fputint(bhv, dat);
	comment("k0nq", dat); fputdouble(k0nq, dat);
	comment("1/(n*n)", dat); fputdouble(inve, dat);
	return;
}

/* input from file */
void SmPiece::read(FILE *dat)
{
	pol = Polarization(fgetint(dat));
	bif = fgetdouble(dat);
	tif = fgetdouble(dat);
	k0 = fgetdouble(dat);
	A = fgetdouble(dat);
	B = fgetdouble(dat);
	g = fgetdouble(dat);
	xa = fgetdouble(dat);
	xb = fgetdouble(dat);
	bhv = FldBhv(fgetint(dat));
	k0nq = fgetdouble(dat);
	inve = fgetdouble(dat);
	return;
}

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

/* initialize */
Mode::Mode()
{
	N = 0;
	piece = NULL;
	special = 0;
	vM = 0;
}
Mode::Mode(const Waveguide& 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 SmPiece[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;
	vM = 0;
}

/* VEIMS: nonstandard TM effective horizontal profile */
Mode::Mode(const Waveguide& g, const Dvector& ec, 
           Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt)
{
	wg = g;
	pol = TM;
	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 SmPiece[N];
	switch(bdtb)
	{
	case OPEN:
		l.x0 = -AWAY; l.x1 = wg.hx(0);
		piece[0].init(l, wg.eps(0), ec(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(l, wg.eps(0), ec(0), wl);
		break;
	}
	for(int i=1; i<=N-2; ++i) piece[i].init(wg.layer(i), 
	                                        wg.eps(i), ec(i), wl);	
	switch(bdtt)
	{
	case OPEN:
		l.x0 = wg.hx(wg.nx); l.x1 = AWAY;
		piece[N-1].init(l, wg.eps(wg.nx+1), ec(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(l, wg.eps(wg.nx+1), ec(wg.nx+1), wl);
		break;
	}
	special = 0;
	vM = 0;
}

/* destroy */
Mode::~Mode()
{
	
	if(piece != NULL) delete[] piece;
	piece = NULL;
	N = 0;
}

/* copy constructor */
Mode::Mode(const Mode& m)
{
	wg = m.wg;	
	pol = m.pol;	
	k0 = m.k0;	
	invommu0 = m.invommu0;	
	invomep0 = m.invomep0;	
	beta = m.beta;
	betaq = m.betaq;
	cbeta = m.cbeta;
	neff = m.neff;
	npcB = m.npcB;
	ord = m.ord;
	ppt = m.ppt;
	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 SmPiece[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];
	zero = m.zero;
	fac = m.fac;
	sfd = m.sfd;
	dnq = m.dnq;
	img = m.img;
	brv = m.brv;
	evi = m.evi;
	special = m.special;
	vM = m.vM;
	ky = m.ky;
	kz = m.kz;
	ckz = m.ckz;
	cVzero = m.cVzero;
	cVsfd = m.cVsfd;
	cVdnq = m.cVdnq;
	cVfac = m.cVfac;
	cVbrv = m.cVbrv;
	cVppt = m.cVppt;
}

/* assignment */
Mode& Mode::operator=(const Mode& m)
{
	if(this != &m)
	{
		wg = m.wg;	
		pol = m.pol;	
		k0 = m.k0;	
		invommu0 = m.invommu0;	
		invomep0 = m.invomep0;	
		beta = m.beta;
		betaq = m.betaq;
		cbeta = m.cbeta;
		neff = m.neff;
		npcB = m.npcB;
		ord = m.ord;
		ppt = m.ppt;
		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 SmPiece[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];
		zero = m.zero;
		fac = m.fac;
		sfd = m.sfd;
		dnq = m.dnq;
		img = m.img;
		brv = m.brv;
		evi = m.evi;
		special = m.special;
		vM = m.vM;
		ky = m.ky;
		kz = m.kz;
		ckz = m.ckz;
		cVzero = m.cVzero;
		cVsfd = m.cVsfd;
		cVdnq = m.cVdnq;
		cVfac = m.cVfac;
		cVbrv = m.cVbrv;
		cVppt = m.cVppt;
	}
	return *this;
}
/* ------------------------------------------------------------------------ */

/* mode solver: transverse resonance evaluation */
double Mode::travres(double tbq)
{
	double phi, phis;
	piece[0].start_tre(tbq, phi, phis, bdt_b); 
	for(int l=1; l<=N-2; ++l) piece[l].cont_tre(tbq, phi, phis); 
	return piece[N-1].finish_tre(tbq, phi, phis, bdt_t); 
}

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

/* node counting: determine the number of modes (returned) with squared 
   propagation constants larger than the present betaq, the return value is
   equal to the order+1 of the highest order mode above the present betaq; 
   meant to be applied to an unfinished mode object */
int Mode::nummodesabove(double tbq)
{
	travres(tbq);
	int z = nodes();
	if(bdt_t != LIMN) return z; 
	if(piece[N-1].bfld(bp_t)*piece[N-1].bder(bp_t)<0.0) ++z;
	return z;
}

/* write data files for travres, nummodesabove & exit; 
   for diagnosis purposes */
void Mode::travresinspect(double bqmin, double bqmax)
{
	int numtrip = 20000;
	Dvector b(numtrip);
	Dvector tr(numtrip);
	Dvector nd(numtrip);
	for(int j=0; j<=numtrip-1; ++j)
	{
		b(j) = bqmin+j/((double) numtrip)*(bqmax-bqmin);
		tr(j) = travres(b(j));
		nd(j) = nummodesabove(b(j));
	}
	toxyf("tr", b, tr);
	toxyf("nd", b, nd);
//	slabmodeerror("travresinspect: stop");
	return;
}


/* localize a zero in the mode transverse resonance function */
double Mode::localize(double bql, double bqu, double lim)
{
	double fu, fl, f;
	double bq;
	if(bqu < bql) slabmodeerror("localize (0)");
	fl = travres(bql);
	fu = travres(bqu);
	if(fl*fu > 0.0) slabmodeerror("localize (1)");
	while(bqu-bql > lim)
	{
		bq = (bql+bqu)*0.5;
		f  = travres(bq);
		if(fl*f < 0.0) { bqu = bq; fu  = f; }
		else           { bql = bq; fl  = f; }
	} 
	bq = (bqu+bql)*0.5;
	f  = travres(bq);
	return bq;
}

/* mode solver: set the relevant variables for an actual modal solution */ 
void Mode::polish(double bq)
{
	betaq = bq;
	if(betaq >= 0.0)
	{
		beta = sqrt(betaq);
		ppt = PROPAG;
	}
	else
	{
		beta = sqrt(-betaq);
		ppt = EVANESC;
	}
	neff = beta/k0;
	double emin = wg.defaultepseffmin();
	double emax = wg.defaultepseffmax();
	if(emin < emax) npcB = (betaq/k0/k0-emin)/(emax-emin); 
	else            npcB = -1.0;
	travres(betaq);
	piece[N-1].polish_ub(bdt_t);
	setfieldmap();
	normalize();
	classify();
	setids();

	Mode m = (*this);
	m.revtravres(betaq);
	m.piece[0].polish_lb(bdt_b);
	m.setfieldmap();
	m.normalize();
	m.classify();
	m.setids();

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

/* 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 Mode::expand()
{
	if(bdt_t == OPEN) slabmodeerror("expand: bdt");
	Waveguide 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;
	SmPiece *npiece;
	npiece = new SmPiece[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].g     = piece[l].g;
		npiece[tl].bhv   = piece[l].bhv;
		npiece[tl].k0nq  = piece[l].k0nq;
		npiece[tl].inve = piece[l].inve;
		if(sym == 1)
		{
			if(piece[l].bhv == HYP)
			{
				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].A;
				npiece[tl].xa = p+(p-piece[l].xa);
				npiece[tl].B  = piece[l].B;
				npiece[tl].xb = p+(p-piece[l].xb);
			}
		}
		else
		{
			if(piece[l].bhv == HYP)
			{
				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].A;
				npiece[tl].xa = p+(p-piece[l].xa);
				npiece[tl].B  = -piece[l].B;
				npiece[tl].xb = p+(p-piece[l].xb);
			}
		}
	}
	delete[] piece;
	N = nN;
	piece = npiece;

	double emin = wg.defaultepseffmin();
	double emax = wg.defaultepseffmax();
	if(emin < emax) npcB = (betaq/k0/k0-emin)/(emax-emin); 
	else            npcB = -1.0;

	normalize();
	setids();
	return;
}

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

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

	SmPiece *npiece;
	npiece = new SmPiece[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].g     = piece[l].g;
		npiece[tl].bhv   = piece[l].bhv;
		npiece[tl].k0nq  = piece[l].k0nq;
		npiece[tl].inve = piece[l].inve;
		if(piece[l].bhv == HYP)
		{
			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].A;
			npiece[tl].xa = p+(p-piece[l].xa);
			npiece[tl].B  = piece[l].B;
			npiece[tl].xb = p+(p-piece[l].xb);
		}
	}
	delete[] piece;
	piece = npiece;
	return;
}

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

/* number of nodes in the principal field */
int Mode::nodes() const
{
	Interval i;
	double dx = 1.0e-10;
	double b, t, p;
	int n = 0;
	for(int l=0; l<=N-1; ++l) 
	{
		i = wg.layer(l);
		if(l == 0       && bdt_b != OPEN) i.x0 = bp_b+dx; 
		if(l == wg.nx+1 && bdt_t != OPEN) i.x1 = bp_t-dx; 
		b = i.x0+dx;
		t = i.x1-dx;
		n += piece[l].nodes(b, t);
		if(l <= N-2)
                {
			p = i.x1;
			double fl = piece[l].bfld(p);
			double fr = piece[l+1].bfld(p);
			double fb = piece[l].bfld(t);
			double ft = piece[l+1].bfld(p+dx);
			if(fl*fr >= 0.0) 
			{
				if(fb*ft < 0.0) ++n;
			}
			else
			{	
				int done = 0;
				double fls = piece[l].bder(p);
				double frs = piece[l+1].bder(p);
				if(   fls > 0.0 && frs > 0.0
				   && fl  < 0.0 && fr  > 0.0)
				{
					++n; 
					done = 1;
				}
				if(   fls > 0.0 && frs > 0.0
				   && fl  > 0.0 && fr  < 0.0)
				{
					if(fb > 0.0) --n;
					if(ft < 0.0) --n;
					++n;
					done = 1;
				}
				if(   fls < 0.0 && frs < 0.0
				   && fl  > 0.0 && fr  < 0.0)
				{
					++n; 
					done = 1;
				}
				if(   fls < 0.0 && frs < 0.0
				   && fl  < 0.0 && fr  > 0.0)
				{
					if(fb < 0.0) --n;
					if(ft > 0.0) --n;
					++n;
					done = 1;
				}
				if(done == 0)
				{
					fprintf(stderr, "nodes@interface:\n");
					fprintf(stderr, " p  = %g\n", p);
					fprintf(stderr, "dx  = %g\n", dx);
					fprintf(stderr, "fb  = %g\n", fb);
					fprintf(stderr, "fl  = %g\n", fl);
					fprintf(stderr, "fr  = %g\n", fr);
					fprintf(stderr, "ft  = %g\n", ft);
					fprintf(stderr, "fl' = %g\n", fls);
					fprintf(stderr, "fr' = %g\n", frs);
					piece[l].write(stderr);
					piece[l+1].write(stderr);
					slabmodeerror("nodes: confused");
				}
			}
                }
	}
	return n;
}

/* number of nodes in the derivative of the principal field */
int Mode::dnodes() const
{
	Interval i;
	int n = 0;
	double dx = 1.0e-10;
	double b, t;
	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; 
		b = i.x0+dx;
		t = i.x1-dx;
		n += piece[l].dnodes(b, t);
		if(l <= N-2) 
		{
			if(piece[l].bder(t)
			   *piece[l+1].bder(i.x1+dx)<0.0) ++n;
		}
	}
	return n;
}

/* classify: determine mode order, if possible,
             otherwise: set the special switch */
#define NONODESLIM 1.0e-6
void Mode::classify()
{
	special = 0;
	ord = nodes();
	Interval i;
	double lim = maxF*NONODESLIM;
	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; 
		if(piece[l].maxf(i) < lim) 
		{
			special = 1;
			return;
		}
	}
	return;
}

/* set the mode id string */
void Mode::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;
	return;
}
/* ... for known mode order */
void Mode::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 */
double Mode::field(double x) const
{
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return 0.0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return 0.0;
	return piece[wg.layeridx(x)].bfld(x); 
}
/* and its derivative */
double Mode::d_field(double x) const
{
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return 0.0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return 0.0;
	return piece[wg.layeridx(x)].bder(x); 
}

/* field values at point x,
   cp: EX - HZ, SX, SY, SZ 
   note the convention for real mode field components */
/* the forward field */
double Mode::field(Fcomp cp, double x) const
{
	if((cp == SX) | (cp == SY)) return 0.0;
	if(bdt_b != OPEN && x<bp_b-S_HDIST) return 0.0;
	if(bdt_t != OPEN && x>bp_t+S_HDIST) return 0.0;
	int l = wg.layeridx(x);
	double e, f;
	if(cp == SZ) 
	{
		if(pol == TE)  
		{
			f = piece[l].bfld(x);
			return 0.5*beta*invommu0*f*f;
		}
		e = wg.eps(l);
		f = piece[l].bfld(x);
		return 0.5*beta*invomep0*f*f/e;
	}
	if(zero(cp)) return 0.0;
	f = (sfd(cp)) ? piece[l].bfld(x) : piece[l].bder(x); 
	if(dnq(cp)) f *= piece[l].inve;
	if(evi(cp)) f *= -1.0;
	return fac(cp)*f;
}

/* field values for forward or backward propagating modes,
   Fcomp: EX - HZ */
double Mode::field(Fcomp cp, Propdir d, double x) const
{
	if(((int) cp) > ((int) HZ)) slabmodeerror("field, Propdir: cp");
	double f = field(cp, x);
	if(d == FORW) return f;
	return (brv(cp)) ? -f : f; 
}

/* factor of phase/amplitude change along a propagation distance z
   d: FORW, BACK selects the forward / backward traveling field */
Complex Mode::efac(Propdir d, double z) const
{
	if(vM) return expmi(ckz(d)*z);
	if(d == FORW)
		return (ppt == PROPAG) ? expi(-beta*z) : exp(-beta*z);
	return (ppt == PROPAG) ? expi(beta*z) : exp(beta*z);
}
/* forward and backward values, for distances dzf and dzb */
void Mode::efac(double dzf, Complex& ef, double dzb, Complex& eb) const
{
	if(vM)
	{	
		ef = expmi(kz*dzf);
		eb =  expi(kz*dzb);
		return;
	}	
	ef = (ppt == PROPAG) ? expi(-beta*dzf) : exp(-beta*dzf);
	eb = (ppt == PROPAG) ? expi( beta*dzb) : exp( beta*dzb);
	return;
}

/* factor of phase/amplitude change along a propagation distance z
   d: FORW, BACK selects the forward / backward traveling field 
   pert: a propagation constant shift */
Complex Mode::efac(Propdir d, Complex pert, double z) const
{
	if(vM) slabmodeerror("efac, pert: vM not implemented");
	Complex cbet;
	if(ppt == PROPAG) cbet = Complex(beta,  0.0)+pert;
	else              cbet = Complex(0.0, -beta)+pert;
	if(d == FORW) return exp(CCI*cbet*(-z));
	else          return exp(CCI*cbet*( z));
}

/* complex mode profiles 
   Fcomp: EX - HZ */
Complex Mode::cfield(Fcomp cp, Propdir d, double x) const
{

	if(vM) return cVfield(cp, d, x);
	double f = field(cp, d, x);
	return (img(cp)) ? Complex(0.0, f) : Complex(f);
}
/* forward and backward field */
void Mode::cfield(Fcomp cp, double x, Complex& ff, Complex& bf) const
{
	if(vM)
	{	
		ff = cVfield(cp, FORW, x);
		bf = ff;
		if(cVbrv(cp)) bf.neg();
		return;
	} 
	double f = field(cp, FORW, x);
	ff = (img(cp)) ? Complex(0.0, f) : Complex(f);
	bf = ff;
	if(brv(cp)) bf.neg();
	return;
}

/* if viewed as a vertically propagating mode: complex mode profiles 
   Fcomp: EX - HZ */
Complex Mode::vertfield(Fcomp cp, Propdir d, double z) const
{
	if(vM) return vertVfield(cp, d, z);
	Complex f;
	switch(cp)
	{
	case EX: f = cfield(EZ, d, z); f.neg(); return f;
	         break;
	case EY: return cfield(EY, d, z); 
	         break;
	case EZ: f = cfield(EX, d, z); f.neg(); return f;
	         break;
	case HX: f = cfield(HZ, d, z); f.neg(); return f;
	         break;
	case HY: return cfield(HY, d, z); 
	         break;
	case HZ: f = cfield(HX, d, z); f.neg(); return f;
	         break;
	default: slabmodeerror("vertfield: cp");
	         break;
	}
	return CC0;
}

/* 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 Mode::cfield(Fcomp cp, Propdir d, double x, double z) const
{
	if(vM) return cVfield(cp, d, x)*expmi(ckz(d)*z);
	return cfield(cp, d, x)*efac(d, z);
}

/* check the continuity of the principal field components,
   returns a maximum error quantity */
double Mode::checkcontinuity() const
{
	double xm, xp;
	double 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 = fabs(p-m)/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 Mode::checkmode() const
{
	if(vM)
	{
		cVcheckmode();
		return;
	}
	Complex ex, ey, ez, hx, hy, hz, cbet;
	Complex s1, s2, s3, s4, s5, s6, s7, s8;
	double s1sum, s2sum, s3sum, s4sum, s5sum, s6sum, s7sum, s8sum;
	Complex dxex, dxey, dxez, dxhx, dxhy, dxhz;
	double x0, x1, xp, s, ds;
	double ommu, omep, eps;
	double ddxxbf, s9, s9sum, bf, dxbf;
	Complex ext, eyt, ezt, hxt, hyt, hzt;
	Complex extsum, eytsum, eztsum, hxtsum, hytsum, hztsum;

	ommu = 1.0/val_invommu0(wg.lambda);
	omep = 1.0/val_invomep0(wg.lambda);
	fprintf(stderr, "Modecheck:\n");
	fprintf(stderr, "beta: %g,  neff: %g,  npcB: %g\n", beta, neff, npcB);
	fprintf(stderr, "pol: T%c,  ", polCHR(pol) );
	if(ppt == PROPAG) fprintf(stderr, "ppt: PROPAG,  ");
	else              fprintf(stderr, "ppt: EVANESC,  ");

	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 = cfield(EX, dir, x0);
	ey = cfield(EY, dir, x0);
	ez = cfield(EZ, dir, x0);
	hx = cfield(HX, dir, x0);
	hy = cfield(HY, dir, x0);
	hz = cfield(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 = cfield(EX, dir, x1);
	ey = cfield(EY, dir, x1);
	ez = cfield(EZ, dir, x1);
	hx = cfield(HX, dir, x1);
	hy = cfield(HY, dir, x1);
	hz = cfield(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");
	if(ppt == PROPAG) cbet = Complex( 0.0, -beta);
	else              cbet = Complex(-beta,  0.0);
	if(dir == BACK) cbet = cbet*(-1.0);
	if(bdt_b == OPEN) x0 = wg.hx(0)-2.0; 
	if(bdt_t == OPEN) x0 = wg.hx(0)+2.0; 
	s = (x1-x0)/200.0;	
	ds = DERIVSTEP/2.0;
	s1sum=s2sum=s3sum=s4sum=s5sum=s6sum=s7sum=s8sum=s9sum=0.0;
	extsum = CC0;
	eytsum = CC0;
	eztsum = CC0;
	hxtsum = CC0;
	hytsum = CC0;
	hztsum = CC0;
	for(xp=x0+s; xp<=x1-s; xp+=s)
	{
		if(fabs(wg.eps(xp+ds)-wg.eps(xp-ds))<1.0e-12)
		{
		eps = wg.eps(xp); 
		ex = cfield(EX, dir, xp);
		ey = cfield(EY, dir, xp);
		ez = cfield(EZ, dir, xp);
		hx = cfield(HX, dir, xp);
		hy = cfield(HY, dir, xp);
		hz = cfield(HZ, dir, xp);

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

	if(ppt == PROPAG)
	{ 	
		if(pol == TE)
		{
			bf = field(EY, xp);
			dxbf = (field(EY, xp+ds)-field(EY, xp-ds))/(2.0*ds);
			ext = CC0;
			eyt = Complex(bf, 0.0);
			ezt = CC0; 
			hxt = Complex(bf, 0.0)*(-1.0)*beta/ommu; 
			hyt = CC0;
			hzt = Complex(0.0, dxbf)/ommu; 
		}
		else
		{
			bf = field(HY, xp);
			dxbf = (field(HY, xp+ds)-field(HY, xp-ds))/(2.0*ds);
			ext = Complex(bf, 0.0)*beta/(omep*eps); 
			eyt = CC0;
			ezt = CCI*Complex(dxbf, 0.0)*(-1.0)/(omep*eps); 
			hxt = CC0;
			hyt = Complex(bf, 0.0);
			hzt = CC0; 
		}
	}
	else
	{ 	
		if(pol == TE)
		{
			bf = field(EY, xp);
			dxbf = (field(EY, xp+ds)-field(EY, xp-ds))/(2.0*ds);
			ext = CC0;
			eyt = Complex(bf, 0.0);
			ezt = CC0; 
			hxt = Complex(0.0, bf)*beta/ommu; 
			hyt = CC0;
			hzt = Complex(0.0, dxbf)/ommu; 
		}
		else
		{
			bf = field(HY, xp);
			dxbf = (field(HY, xp+ds)-field(HY, xp-ds))/(2.0*ds);
			ext = Complex(bf, 0.0)*beta/(omep*eps); 
			eyt = CC0;
			ezt = Complex(dxbf, 0.0)/(omep*eps); 
			hxt = CC0;
			hyt = Complex(0.0, bf);
			hzt = CC0; 
		}
	}
	if(dir == BACK)
	{
		ezt = ezt*(-1.0); 
		hxt = hxt*(-1.0); 
		hyt = hyt*(-1.0); 
	}
	if(extsum.abs() < (ext-ex).abs()) extsum = ext-ex;
	if(eytsum.abs() < (eyt-ey).abs()) eytsum = eyt-ey;
	if(eztsum.abs() < (ezt-ez).abs()) eztsum = ezt-ez;
	if(hxtsum.abs() < (hxt-hx).abs()) hxtsum = hxt-hx;
	if(hytsum.abs() < (hyt-hy).abs()) hytsum = hyt-hy;
	if(hztsum.abs() < (hzt-hz).abs()) hztsum = hzt-hz;
		
		s1 = cbet*ey-CCI*hx*ommu; 
		s2 = dxez-cbet*ex-CCI*hy*ommu; 
		s3 = (dxey*(-1.0))-CCI*hz*ommu; 
		s4 = cbet*hy+CCI*ex*omep*eps;
		s5 = dxhz-cbet*hx+CCI*ey*omep*eps;
		s6 = (dxhy*(-1.0))+CCI*ez*omep*eps;
		s7 = dxex+cbet*ez;
		s8 = dxhx+cbet*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);
		}
		if(ppt == PROPAG) s9 = ddxxbf+(k0*k0*eps-beta*beta)*bf;
		else              s9 = ddxxbf+(k0*k0*eps+beta*beta)*bf;
		if(s9 > s9sum) s9sum = s9;
		}
	} 
	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, "bf -> EX: %g + i %g\n", extsum.re, extsum.im);
	fprintf(stderr, "bf -> EY: %g + i %g\n", eytsum.re, eytsum.im);
	fprintf(stderr, "bf -> EZ: %g + i %g\n", eztsum.re, eztsum.im);
	fprintf(stderr, "bf -> HX: %g + i %g\n", hxtsum.re, hxtsum.im);
	fprintf(stderr, "bf -> HY: %g + i %g\n", hytsum.re, hytsum.im);
	fprintf(stderr, "bf -> HZ: %g + i %g\n", hztsum.re, hztsum.im);
	}
	fprintf(stderr, "\n");
	return;
}
void Mode::cVcheckmode() const
{
	Complex ex, ey, ez, hx, hy, hz, cbet;
	Complex s1, s2, s3, s4, s5, s6, s7, s8;
	double s1sum, s2sum, s3sum, s4sum, s5sum, s6sum, s7sum, s8sum;
	Complex dxex, dxey, dxez, dxhx, dxhy, dxhz;
	double x0, x1, xp, s, ds;
	double ommu, omep, eps;

	ommu = 1.0/val_invommu0(wg.lambda);
	omep = 1.0/val_invomep0(wg.lambda);
	fprintf(stderr, "Modecheck (vM): %s  ky = %g+i%g, kz = %g+i%g\n", ids, ky.re, ky.im, kz.re, kz.im);
	fprintf(stderr, "betaq = %g, beta: %g,  neff: %g\n", betaq, beta, neff);
	fprintf(stderr, "pol: T%c,  ", polCHR(pol) );
	if(ppt == PROPAG) fprintf(stderr, "ppt: PROPAG,  ");
	else              fprintf(stderr, "ppt: EVANESC,  ");

	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 = cfield(EX, dir, x0);
	ey = cfield(EY, dir, x0);
	ez = cfield(EZ, dir, x0);
	hx = cfield(HX, dir, x0);
	hy = cfield(HY, dir, x0);
	hz = cfield(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 = cfield(EX, dir, x1);
	ey = cfield(EY, dir, x1);
	ez = cfield(EZ, dir, x1);
	hx = cfield(HX, dir, x1);
	hy = cfield(HY, dir, x1);
	hz = cfield(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");
	if(bdt_b == OPEN) x0 = wg.hx(0)-2.0; 
	if(bdt_t == OPEN) x0 = wg.hx(0)+2.0; 
	s = (x1-x0)/200.0;	
	ds = DERIVSTEP/2.0;
	s1sum=s2sum=s3sum=s4sum=s5sum=s6sum=s7sum=s8sum=0.0;
	for(xp=x0+s; xp<=x1-s; xp+=s)
	{
		if(fabs(wg.eps(xp+ds)-wg.eps(xp-ds))<1.0e-12)
		{
		eps = wg.eps(xp); 
		ex = cfield(EX, dir, xp);
		ey = cfield(EY, dir, xp);
		ez = cfield(EZ, dir, xp);
		hx = cfield(HX, dir, xp);
		hy = cfield(HY, dir, xp);
		hz = cfield(HZ, dir, xp);

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

		s1 = CCmI*ky*ez      +CCI*ckz(dir)*ey + CCI*ommu*hx;
		s2 = CCmI*ckz(dir)*ex-dxez            + CCI*ommu*hy;
		s3 = dxey            +CCI*ky*ex       + CCI*ommu*hz;

		s4 = CCmI*ky*hz      +CCI*ckz(dir)*hy - CCI*omep*eps*ex;
		s5 = CCmI*ckz(dir)*hx-dxhz            - CCI*omep*eps*ey;
		s6 = dxhy            +CCI*ky*hx       - CCI*omep*eps*ez;

		s7 = dxex-CCI*ky*ey-CCI*ckz(dir)*ez;
		s8 = dxhx-CCI*ky*hy-CCI*ckz(dir)*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();
		}
	} 
	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, "\n");
	return;
}

/* principal component, integrate a product of field or
   derivatives over interval i */
double Mode::integrate(FldorDer fod1, FldorDer fod2, Interval i) const
{
	return integrate(fod1, *this, fod2, i);
}

/* integration of a fieldproduct over the interval [i.x0, i.x1] 
   cp1, cp2 != SZ */
double Mode::integrate(Fcomp cp1, Fcomp cp2, Interval i) const 
{
	return integrate(cp1, *this, cp2, i);
}
	
/* integration of a fieldproduct over the interval [i.x0, i.x1], 
   for bidirectional mode versions, 
   d1, d2 = FORW, BACK,  
   cc1, cc2 == 1: take the complex conjugate of the component, 
   cp1, cp2 != SX, SY, SZ */
Complex Mode::integrate(Fcomp cp1, Propdir d1, int cc1,  
                        Fcomp cp2, Propdir d2, int cc2, Interval i) const 
{
	return integrate(d1, cp1, cc1, *this, d2, cp2, cc2, i);
}
	
/* ------------------------------------------------------------------------ */

/* z component of the Poyntingvector, integrated along the x-axis, 
   the modal power; ppt EVANESC: normalization */
#define SELECT_K_COMP 1.0e-14
double Mode::power() const
{
	if(vM) 
	{
		if(cVppt == PROPAG) return cVoverlap(FORW, (*this), FORW).abs();
		else                return cVoverlap(FORW, (*this), BACK).abs();
	}
	double s = 0.0;
	double f;

	f = (bdt_b == OPEN) ? piece[0].intbfq(-AWAY, wg.hx(0))
	                    : piece[0].intbfq(bp_b, wg.hx(0));
	s += (pol == TE) ? f : f*piece[0].inve;
	for(int l=1; l<=N-2; ++l)
	{
		f = piece[l].intbfq();
		s += (pol == TE) ? f : f*piece[l].inve;
	} 
	f = (bdt_t == OPEN) ? piece[N-1].intbfq(wg.hx(wg.nx), AWAY)
	                    : piece[N-1].intbfq(wg.hx(wg.nx), bp_t);
	s += (pol == TE) ? f : f*piece[N-1].inve;
	return (pol == TE) ? s*0.5*beta*invommu0 : s*0.5*beta*invomep0;	
}

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


/* principal field component, set maximum absolute level */
void Mode::setmaxF()
{
	Interval i;
	double m;
	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);
		if(m > maxF) maxF = m;
	}
	return;
}

/* for guided modes only: extent of the mode profile;
   outside the returned interval the field strength
   (absolute value, principal component) is lower than fac * maxF;
   typically fac << 1.0 */
Interval Mode::extent(double fac) const
{
	Interval e;
	if(bdt_b != OPEN) return Interval(bp_b, bp_t);
	if(bdt_t != OPEN) return Interval(bp_b, bp_t);
	double xl = wg.hx(0);
	double xr = wg.hx(wg.nx);
	if(piece[0].bhv != HYP)  return Interval(xl, xr);
	if(piece[N-1].bhv != HYP)  return Interval(xl, xr);
	double g, f;
	g = piece[0].g;
	f = fabs(piece[0].bfld(xl));
	e.x0 = xl+1.0/g*log(fabs(fac)*maxF/f);
	g = piece[N-1].g;
	f = fabs(piece[N-1].bfld(xr));
	e.x1 = xr-1.0/g*log(fabs(fac)*maxF/f);
	return e;
}

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

/* normalize mode to power() == 1, if possible */
void Mode::normalize()
{
	normalize(1.0);
	return;
}
/* normalize mode to power() == p, if possible */
void Mode::normalize(double p)
{
	double fac = fabs(power());
	if(fac < 1.0e-10)
	{ 
		fprintf(stderr, "Mode::normalize: vM=%d (?)\n", vM);
		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;
	setmaxF();
	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 Mode::writeprofile(Fcomp cp, Afo foa, Interval disp, int np,
                        char ext0, char ext1) const
{
	if(np <= 0) slabmodeerror("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 ORG:
		fprintf(stderr, "   %c%c [%g, %g] ", 
		                fldchr(cp), cpchr(cp), x0, x1);
		break;
	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:
		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(cfield(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(cfield(cp, FORW, x(j)), foa);
		}
	}
	toxyf(name, x, f);
	return;
}
/* as before, with default parameters */
void Mode::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), ORG, 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 Mode::plot(Fcomp cp, Afo foa, Interval disp, int np, 
                char ext0, char ext1, char pltype) const
{
	FILE *dat;
	double x, dx, t;
	double x0, x1;
	double xbd;
	double epsold, epsnew;
	double minf, maxf, f;
	int nsec, nbd;
	char name[13] = "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 ORG:
		fprintf(stderr, "   %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case MOD:
		fprintf(stderr, "   |%c%c| [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case SQR:
		fprintf(stderr, "   |%c%c|^2 [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case REP:
		fprintf(stderr, "   Re %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	case IMP:
		fprintf(stderr, "   Im %c%c [%g, %g] (%d) >> %s\n", 
				fldchr(cp), cpchr(cp), x0, x1, np, name);
		break;
	}

	Cmatrix cfld(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);
	pos(nsec, nump(nsec)) = x;  
	cfld(nsec, nump(nsec)) = cfield(cp, FORW, x);  
	++nump(nsec);
	for(int i=1; i<=np; ++i)
	{
		x = x0+((double) i)*dx/((double) np);
		epsnew = wg.eps(x);
		if(fabs(epsnew-epsold)<1.0e-10)
		{
			pos(nsec, nump(nsec)) = x;  
			cfld(nsec, nump(nsec)) = cfield(cp, FORW, x);  
			++nump(nsec);
		}
		else
		{
			xbd = (wg.layer(x)).x0; 
			pos(nsec, nump(nsec)) = xbd;  
			cfld(nsec, nump(nsec)) = 	
				cfield(cp, FORW, xbd-dx*SLHDISTF);  
			++nump(nsec);
			bd(nbd) = xbd; 
			++nbd;
			epsold = epsnew;
			ri(nbd) = sgn_sqrt(epsold);
			++nsec;
			nump(nsec) = 0;
			pos(nsec, nump(nsec)) = xbd;  
			cfld(nsec, nump(nsec)) = 	
				cfield(cp, FORW, xbd+dx*SLHDISTF);  
			++nump(nsec);
			pos(nsec, nump(nsec)) = x;  
			cfld(nsec, nump(nsec)) = cfield(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 = fld.min();
	maxf = fld.max();
	if(minf >= maxf) { minf -= 1.0;	maxf = minf+2.0; }

	if(pltype == 'L')
	{
		switch(foa)
		{
			case MOD: case SQR:
				minf = 0.0;
	 			break;
			case ORG: case REP: case IMP:
				f = fabs(minf); maxf = fabs(maxf);
				if(f > maxf) maxf = f;
				minf = -maxf;
	 			break;
		}
		mlout_geo(dat, 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 Mode::plot(char ext0, char ext1, char pltype) 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, pltype);
	return;
}
void Mode::plot(char ext0, char ext1) const
{
	plot(ext0, ext1, 'L');
	return;
}

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

/* compare the present object to m, return 1 if equal, otherwise return 0; 
   both modes are assumed to be properly normalized eigenfunctions */
#define COMPTOL_BETA    1.0e-7
#define COMPTOL_BP      1.0e-7
int Mode::equal(const Mode& m) const
{
	if(this == &m) return 1;
	if(wg.equal(m.wg) != 1) return 0;
	if(pol != m.pol) return 0;
	if(fabs(beta-m.beta) > COMPTOL_BETA) return 0; 
	if(fabs(betaq-m.betaq) > COMPTOL_BETA) return 0; 
	if(ppt != m.ppt) 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(vM != m.vM) return 0;
	if((ky-m.ky).sqabs() > COMPTOL_BETA) return 0; 
	if((kz-m.kz).sqabs() > COMPTOL_BETA) return 0; 
	return 1;
}

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

/* output to file */
void Mode::write(FILE *dat) const
{
	comment("Mode", 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("beta", dat); fputdouble(beta, dat);
	comment("betaq", dat); fputdouble(betaq, dat);
	comment("cbeta", dat); cbeta.write(dat);
	comment("neff", dat); fputdouble(neff, dat);
	comment("npcB", dat); fputdouble(npcB, dat);
	comment("ord", dat); fputint(ord, dat);
	comment("ppt", dat); fputint(ppt, 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("zero", dat); zero.write(dat);
	comment("fac", dat); fac.write(dat);
	comment("sfd", dat); sfd.write(dat);
	comment("dnq", dat); dnq.write(dat);
	comment("img", dat); img.write(dat);
	comment("brv", dat); brv.write(dat);
	comment("evi", dat); evi.write(dat);
	comment("special", dat); fputint(special, dat);
	comment("vM", dat); fputint(vM, dat);
	comment("ky", dat); ky.write(dat);
	comment("kz", dat); kz.write(dat);
	comment("ckz", dat); ckz.write(dat);
	comment("cVzero", dat); cVzero.write(dat);
	comment("cVsfd", dat); cVsfd.write(dat);
	comment("cVdnq", dat); cVdnq.write(dat);
	comment("cVfac", dat); cVfac.write(dat);
	comment("cVbrv", dat); cVbrv.write(dat);
	comment("cVppt", dat); fputint(cVppt, dat);
	return;
}

/* input from file */
void Mode::read(FILE *dat)
{
	wg.read(dat);
	pol = Polarization(fgetint(dat));
	k0 = fgetdouble(dat);
	invommu0 = fgetdouble(dat);
	invomep0 = fgetdouble(dat);
	beta = fgetdouble(dat);
	betaq = fgetdouble(dat);
	cbeta.read(dat);
	neff = fgetdouble(dat);
	npcB = fgetdouble(dat);
	ord = fgetint(dat);
	ppt = Propagation_type(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 SmPiece[N];
	for(int l=0; l<=N-1; ++l) piece[l].read(dat);
	for(int i=0; i<=7; ++i) ids[i] = fgetchar(dat);
	zero.read(dat);
	fac.read(dat);
	sfd.read(dat);
	dnq.read(dat);
	img.read(dat);
	brv.read(dat);
	evi.read(dat);
	special = fgetint(dat);
	vM = fgetint(dat);
	ky.read(dat);
	kz.read(dat);
	ckz.read(dat);
	cVzero.read(dat);
	cVsfd.read(dat);
	cVdnq.read(dat);
	cVfac.read(dat);
	cVbrv.read(dat);
	cVppt = Propagation_type(fgetint(dat));
	return;
}

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

/* input from default file */
void Mode::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, "beta: %.10g,  neff: %.10g,  npcB: %.10g\n",
			 beta, neff, npcB); 
	return;
}

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

/* small permittivity perturbation: first order propagation constant
   shift due to perturbation p, for propagation direction dir */
Complex Mode::phaseshift(Perturbation p, Propdir dir)
{
	Complex pt, nrm, db = CC0;
	switch(dir)
	{
	case FORW:
		// xx
		pt = p.e(0, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, BACK, 0, EX, FORW, 0, p.pr); 	
		// xy
		pt = p.e(0, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, BACK, 0, EY, FORW, 0, p.pr); 	
		// xz
		pt = p.e(0, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, BACK, 0, EZ, FORW, 0, p.pr); 	
		// yx
		pt = p.e(1, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, BACK, 0, EX, FORW, 0, p.pr); 	
		// yy
		pt = p.e(1, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, BACK, 0, EY, FORW, 0, p.pr); 	
		// yz
		pt = p.e(1, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, BACK, 0, EZ, FORW, 0, p.pr); 	
		// zx
		pt = p.e(2, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, BACK, 0, EX, FORW, 0, p.pr); 	
		// zy
		pt = p.e(2, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, BACK, 0, EY, FORW, 0, p.pr); 	
		// zz
		pt = p.e(2, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, BACK, 0, EZ, FORW, 0, p.pr); 	
		
		nrm = integrate(EX, FORW, 0, HY, FORW, 0, Interval(-AWAY, AWAY))
		    - integrate(EY, FORW, 0, HX, FORW, 0, Interval(-AWAY, AWAY));
		nrm = nrm*2.0;
		break;
	case BACK:
		// xx
		pt = p.e(0, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, FORW, 0, EX, BACK, 0, p.pr); 	
		// xy
		pt = p.e(0, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, FORW, 0, EY, BACK, 0, p.pr); 	
		// xz
		pt = p.e(0, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EX, FORW, 0, EZ, BACK, 0, p.pr); 	
		// yx
		pt = p.e(1, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, FORW, 0, EX, BACK, 0, p.pr); 	
		// yy
		pt = p.e(1, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, FORW, 0, EY, BACK, 0, p.pr); 	
		// yz
		pt = p.e(1, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EY, FORW, 0, EZ, BACK, 0, p.pr); 	
		// zx
		pt = p.e(2, 0); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, FORW, 0, EX, BACK, 0, p.pr); 	
		// zy
		pt = p.e(2, 1); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, FORW, 0, EY, BACK, 0, p.pr); 	
		// zz
		pt = p.e(2, 2); 
		if(pt.re != 0.0 || pt.im != 0.0) 
			db += pt*integrate(EZ, FORW, 0, EZ, BACK, 0, p.pr); 	
		
		nrm = integrate(EX, BACK, 0, HY, BACK, 0, Interval(-AWAY, AWAY))
		    - integrate(EY, BACK, 0, HX, BACK, 0, Interval(-AWAY, AWAY));
		nrm = nrm*2.0;
		break;
	}
	db = db/nrm*(k0*INVSQRTMU0/INVSQRTEP0);
	return db; 
}

/* small interface displacement: first order propagation constant shift
   due to a shift of the dielectric interface hx(b) in the modes
   waveguide by a distance dx, for propagation direction dir */
Complex Mode::phaseshift(int b, double dx, Propdir dir)
{
	Complex f, nrm, db;
	if(b<0 || b>wg.nx) slabmodeerror("phaseshift: b");
	double x = wg.hx(b);
	db = cfield(EX, dir, x+S_HDIST)*cfield(EX, dir, x-S_HDIST);
	f = cfield(EY, dir, x); db += f*f;
	f = cfield(EZ, dir, x); db -= f*f;
	db *= (wg.eps(b)-wg.eps(b+1));
	nrm = integrate(EX, dir, 0, HY, dir, 0, Interval(-AWAY, AWAY))
	    - integrate(EY, dir, 0, HX, dir, 0, Interval(-AWAY, AWAY));
	nrm *= 2.0;
	db = db/nrm*(k0*INVSQRTMU0/INVSQRTEP0);
	return db*dx; 
}

/* small displacement of a boundary position: 
   first order propagation constant shift due to a shift dx of the 
   boundary position at the bottom (p == 'b') or at the top (p == 't')
   of the computational window, for propagation direction dir, 
   so far implemented as a finite difference expression only ... */
Complex Mode::phaseshift_bdp(char p, double dx, Propdir dir)
{
	if(p != 'b' &&  p != 't') slabmodeerror("phaseshift_bdp: p");
	double fd = SLAMS_BQTOL*1.0e+5;
	ModeArray map;
	double bpb_p = bp_b, bpt_p = bp_t;	
	if(p == 'b') bpb_p = bp_b+fd;
	else         bpt_p = bp_t+fd;
	map.clear();
	double dbq = fabs(betaq);
	if(dbq < 10.0) dbq = SLAMS_BQSEP*100.0;
	else dbq = dbq*SLAMS_BQSEP*100.0; 
	while(map.num <= 0 && dbq < k0*k0*4.0)
	{
		findmodes(wg, pol, bdt_b, bpb_p, bdt_t, bpt_p, 
		          betaq-dbq, betaq+dbq, map, 0);
		dbq *= 10.0;
	}
	if(map.num <= 0) slabmodeerror("phaseshift_bdp: findmodes");
	int nmi = 0; 
	double dst = (map(0).cbeta(FORW) - cbeta(FORW)).abs();
	int j = 1;
	while(j <= map.num-1)
	{
		double d = (map(j).cbeta(FORW) - cbeta(FORW)).abs();
		if(d < dst) { nmi = j; dst = d; }
		++j;
	}
	Complex db = (map(nmi).cbeta(dir)-cbeta(dir))/fd;
	return db*dx; 
}

/* small alteration of the vacuum wavelength: first order propagation
   constant shift due to a shift of the vacuum wavelength by dl,
   for the version of the mode with propagation direction dir */
Complex Mode::phaseshift(double dl, Propdir dir)
{
	Complex nrm, db, fqs;
	db = CC0;
	for(int l=0; l<=wg.nx+1; ++l)
	{
		fqs  = integrate(EX, dir, 0, EX, dir, 0, wg.layer(l));
		fqs += integrate(EY, dir, 0, EY, dir, 0, wg.layer(l));
		fqs -= integrate(EZ, dir, 0, EZ, dir, 0, wg.layer(l));
		db += wg.eps(l)*fqs;
	}
	db *= k0*INVSQRTMU0/INVSQRTEP0;
	fqs  = integrate(HX, dir, 0, HX, dir, 0, Interval(-AWAY, AWAY));
	fqs += integrate(HY, dir, 0, HY, dir, 0, Interval(-AWAY, AWAY));
	fqs -= integrate(HZ, dir, 0, HZ, dir, 0, Interval(-AWAY, AWAY));
	db += fqs*k0*INVSQRTEP0/INVSQRTMU0;
	nrm = integrate(EX, dir, 0, HY, dir, 0, Interval(-AWAY, AWAY))
	    - integrate(EY, dir, 0, HX, dir, 0, Interval(-AWAY, AWAY));
	db = db/(nrm*(-2.0)*wg.lambda);
	return db*dl; 
/* alternative evaluation as interface- and boundary-shifts, 
   via homogeneity arguments */
//	Complex db = cbeta(dir);
//	for(int b=0; b<=wg.nx; ++b) db += phaseshift(b, 1.0, dir)*wg.hx(b);
//	if(bdt_b != OPEN) db += phaseshift_bdp('b', 1.0, dir)*bp_b;
//	if(bdt_t != OPEN) db += phaseshift_bdp('t', 1.0, dir)*bp_t;
//	db *= (-1.0/wg.lambda);	
//	return db*dl; 
}
/* ... the corresponding expression for the alteration of the effective mode
   index due to a vacuum wavelength shift */
Complex Mode::neffshift(double dl, Propdir dir)
{
	Complex db = phaseshift(dl, dir); 
	return (db*wg.lambda+cbeta(dir)*dl)/(2.0*PI);
}

/* translate mode by dx */
void Mode::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;
}

/* --- mode overlaps ----------------------------------------------------- */

/* product 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, alternatively forward or backward propagating */
Complex overlap(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2)
{
	if(m1.vM) return m1.cVoverlap(d1, m2, d2);
	return m1.overlap(d1, m2, d2);
} 

/* product 
   0.25\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dx 
   (non complex conjugate fields)
   of two modes, alternatively forward or backward propagating */
Complex overlap_nc(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2)
{
	if(m1.vM) slabmodeerror("overlap_nc: vM not implemented");
	return m1.overlap_nc(d1, m2, d2);
} 



/* product 
   0.5\int(E_1x H_2y - E_1y H_2x)dx 
   (non complex conjugate fields)
   of two modes, alternatively forward or backward propagating */
Complex overlap_ncns(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2)
{
	if(m1.vM) slabmodeerror("overlap_ncns: vM not implemented");
	return m1.overlap_ncns(d1, m2, d2);
} 

/* for the BEP implementation: products 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, combinations of forward and backward propagating versions */
void bidir_overlaps(const Mode& m1, const Mode& m2, 
                    Complex& off, Complex& ofb, Complex& obf, Complex& obb)
{
	if(m1.vM) 
	{	
		cVbidir_overlaps(m1, m2, off, ofb, obf, obb);
		return;
	}
	double exhymeyhx, hyexmhxey;
	Interval i(-AWAY, AWAY);
	exhymeyhx = m1.integrate(EX, m2, HY, i)-m1.integrate(EY, m2, HX, i);
	hyexmhxey = m1.integrate(HY, m2, EX, i)-m1.integrate(HX, m2, EY, i);

	off = CC0;
	if(m2.ppt == PROPAG)  off += exhymeyhx; 
	else                  off += Complex(0.0, exhymeyhx); 
	if(m1.ppt == PROPAG)  off += hyexmhxey; 
	else                  off -= Complex(0.0, hyexmhxey); 
	off *= 0.25;

	ofb = CC0;
	if(m2.ppt == PROPAG)  ofb -= exhymeyhx; 
	else                  ofb -= Complex(0.0, exhymeyhx); 
	if(m1.ppt == PROPAG)  ofb += hyexmhxey; 
	else                  ofb -= Complex(0.0, hyexmhxey); 
	ofb *= 0.25;

	obf = CC0;
	if(m2.ppt == PROPAG)  obf += exhymeyhx; 
	else                  obf += Complex(0.0, exhymeyhx); 
	if(m1.ppt == PROPAG)  obf -= hyexmhxey; 
	else                  obf += Complex(0.0, hyexmhxey); 
	obf *= 0.25;

	obb = CC0;
	if(m2.ppt == PROPAG)  obb -= exhymeyhx; 
	else                  obb -= Complex(0.0, exhymeyhx); 
	if(m1.ppt == PROPAG)  obb -= hyexmhxey; 
	else                  obb += Complex(0.0, hyexmhxey); 
	obb *= 0.25;

	return;
} 

/* integrals of products of profile components belonging to different modes */

/* principal components, integrate a product of fields or derivatives 
   over interval i */
double Mode::integrate(FldorDer fod1, const Mode& m, FldorDer fod2, Interval i) const
{
	double x, xt, xt1, xt2;
	double 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 = 0.0;
	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;
		if(fod1 == FLD)
		{
			z = (fod2 == FLD) ? piece[l].intff(m.piece[lm], iv)
			                  : piece[l].intfd(m.piece[lm], iv);
		}
		else
		{
			z = (fod2 == FLD) ? piece[l].intdf(m.piece[lm], iv)
			                  : piece[l].intdd(m.piece[lm], iv);
		}
		s += z;
		x = xt;
	}
	return s;
}


/* integrate a fieldproduct of two modes over Interval i, real version */
double Mode::integrate(Fcomp cp, const Mode& m, Fcomp cpm, Interval i) const
{
	if(((int) cp) > ((int) HZ)) slabmodeerror("integrate: cp");
	if(zero(cp)) return 0.0;
	if(m.zero(cpm)) return 0.0;

	double x, xt, xt1, xt2;
	double 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 = 0.0;
	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;
		if(sfd(cp))
		{
			z = (m.sfd(cpm)) ? piece[l].intff(m.piece[lm], iv)
			                 : piece[l].intfd(m.piece[lm], iv);
		}
		else
		{
			z = (m.sfd(cpm)) ? piece[l].intdf(m.piece[lm], iv)
			                 : piece[l].intdd(m.piece[lm], iv);
		}
		if(  dnq(cp )) z *= piece[l].inve;
		if(m.dnq(cpm)) z *= m.piece[lm].inve;
		s += z;
		x = xt;
	}
	if(  evi(cp )) s *= -1.0;
	if(m.evi(cpm)) s *= -1.0;
	return fac(cp)*m.fac(cpm)*s;
}

/* product 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   with another mode, alternatively forward or backward propagating */
Complex Mode::overlap(Propdir d1, const Mode& m2, Propdir d2) const
{
	if(vM) return cVoverlap(d1, m2, d2);
	Complex s = CC0;
	double st;
	Interval i(-AWAY, AWAY);
	st = integrate(EX, m2, HY, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s += st; 
	else                  s += Complex(0.0, st); 
	st = integrate(EY, m2, HX, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s -= st; 
	else                  s -= Complex(0.0, st); 
	st = integrate(HY, m2, EX, i);
	if(d1 == BACK) st *= -1.0;
	if(ppt == PROPAG)  s += st; 
	else                  s -= Complex(0.0, st); 
	st = integrate(HX, m2, EY, i);
	if(d1 == BACK) st *= -1.0;
	if(ppt == PROPAG)  s -= st; 
	else                  s += Complex(0.0, st); 
	return s*0.25;
} 
/* ... for vTE / vTM modes */
Complex Mode::cVoverlap(Propdir d1, const Mode& m2, Propdir d2) const
{
	Complex ex1hy2, ey1hx2, hy1ex2, hx1ey2;
	Interval i(-AWAY, AWAY);
	ex1hy2 = cVintegrate(EX, m2, HY, i);
	ey1hx2 = cVintegrate(EY, m2, HX, i);
	hy1ex2 = cVintegrate(HY, m2, EX, i);
	hx1ey2 = cVintegrate(HX, m2, EY, i);

	if(d1 == BACK) if(   cVbrv(EX)) ex1hy2.neg();
	if(d2 == BACK) if(m2.cVbrv(HY)) ex1hy2.neg();

	if(d1 == BACK) if(   cVbrv(EY)) ey1hx2.neg();
	if(d2 == BACK) if(m2.cVbrv(HX)) ey1hx2.neg();

	if(d1 == BACK) if(   cVbrv(HY)) hy1ex2.neg();
	if(d2 == BACK) if(m2.cVbrv(EX)) hy1ex2.neg();

	if(d1 == BACK) if(   cVbrv(HX)) hx1ey2.neg();
	if(d2 == BACK) if(m2.cVbrv(EY)) hx1ey2.neg();

	return (ex1hy2-ey1hx2+hy1ex2-hx1ey2)*0.25;
}


/* product 
   0.25\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dx 
   (non complex conjugate fields)
   with another mode, alternatively forward or backward propagating */
Complex Mode::overlap_nc(Propdir d1, const Mode& m2, Propdir d2) const
{
	if(vM) slabmodeerror("overlap_nc: vM not implemented");
	Complex s = CC0;
	double st;
	Interval i(-AWAY, AWAY);
	st = integrate(EX, m2, HY, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s += st; 
	else                  s += Complex(0.0, st); 
	st = integrate(EY, m2, HX, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s -= st; 
	else                  s -= Complex(0.0, st); 
	st = integrate(HY, m2, EX, i);
	if(d1 == BACK) st *= -1.0;
	if(ppt == PROPAG)  s += st; 
	else                  s += Complex(0.0, st); 
	st = integrate(HX, m2, EY, i);
	if(d1 == BACK) st *= -1.0;
	if(ppt == PROPAG)  s -= st; 
	else                  s -= Complex(0.0, st); 
	return s*0.25;
} 

/* product 
   0.5\int(E_1x H_2y - E_1y H_2x)dx 
   (non complex conjugate fields)
   with another mode, alternatively forward or backward propagating */
Complex Mode::overlap_ncns(Propdir d1, const Mode& m2, Propdir d2) const
{
	if(vM) slabmodeerror("overlap_ncns: vM not implemented");
	Complex s = CC0;
	double st;
	Interval i(-AWAY, AWAY);
	st = integrate(EX, m2, HY, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s += st; 
	else                  s += Complex(0.0, st); 
	st = integrate(EY, m2, HX, i);
	if(d2 == BACK) st *= -1.0;
	if(m2.ppt == PROPAG)  s -= st; 
	else                  s -= Complex(0.0, st); 
	return s*0.5;
} 

/* field evaluation, determine map: Fcomp -> field */
void Mode::setfieldmap() 
{
	if((int) EX != 0) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) EY != 1) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) EZ != 2) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) HX != 3) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) HY != 4) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) HZ != 5) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) SX != 6) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) SY != 7) slabmodeerror("Fcomp <-> int conv. corrupted");
        if((int) SZ != 8) slabmodeerror("Fcomp <-> int conv. corrupted");
	zero = Ivector(6);
	fac  = Dvector(6);
	sfd  = Ivector(6);
	dnq  = Ivector(6);
	img  = Ivector(6);
	brv  = Ivector(6);
	evi  = Ivector(6);
	switch(pol)
	{
	case TE:	
		zero((int) EX) = 1; 
		zero((int) EY) = 0; 
		zero((int) EZ) = 1; 
		zero((int) HX) = 0; 
		zero((int) HY) = 1; 
		zero((int) HZ) = 0; 

		fac((int) EX) = 0.0;
		fac((int) EY) = 1.0; 
		fac((int) EZ) = 0.0;
		fac((int) HX) = -beta*invommu0; 
		fac((int) HY) = 0.0;
		fac((int) HZ) = invommu0; 

		sfd((int) EX) = 1;
		sfd((int) EY) = 1;
		sfd((int) EZ) = 1;
		sfd((int) HX) = 1;
		sfd((int) HY) = 1;
		sfd((int) HZ) = 0;

		dnq((int) EX) = 0;
		dnq((int) EY) = 0;
		dnq((int) EZ) = 0;
		dnq((int) HX) = 0;
		dnq((int) HY) = 0;
		dnq((int) HZ) = 0;

		switch(ppt)
		{
		case PROPAG:  
			evi((int) EX) = 0;
			evi((int) EY) = 0;
			evi((int) EZ) = 0;
			evi((int) HX) = 0;
			evi((int) HY) = 0;
			evi((int) HZ) = 0;
			break;
		case EVANESC:
			evi((int) EX) = 0;
			evi((int) EY) = 0;
			evi((int) EZ) = 0;
			evi((int) HX) = 1;
			evi((int) HY) = 0;
			evi((int) HZ) = 0;
			break;
		}
		break;

	case TM:
		zero((int) EX) = 0; 
		zero((int) EY) = 1; 
		zero((int) EZ) = 0; 
		zero((int) HX) = 1; 
		zero((int) HY) = 0; 
		zero((int) HZ) = 1; 

		fac((int) EX) = beta*invomep0; 
		fac((int) EY) = 0.0; 
		fac((int) EZ) = -invomep0; 
		fac((int) HX) = 0.0; 
		fac((int) HY) = 1.0;
		fac((int) HZ) = 0.0; 

		sfd((int) EX) = 1;
		sfd((int) EY) = 1;
		sfd((int) EZ) = 0;
		sfd((int) HX) = 1;
		sfd((int) HY) = 1;
		sfd((int) HZ) = 1;

		dnq((int) EX) = 1;
		dnq((int) EY) = 0;
		dnq((int) EZ) = 1;
		dnq((int) HX) = 0;
		dnq((int) HY) = 0;
		dnq((int) HZ) = 0;

		switch(ppt)
		{
		case PROPAG:  
			evi((int) EX) = 0;
			evi((int) EY) = 0;
			evi((int) EZ) = 0;
			evi((int) HX) = 0;
			evi((int) HY) = 0;
			evi((int) HZ) = 0;
			break;
		case EVANESC:
			evi((int) EX) = 0;
			evi((int) EY) = 0;
			evi((int) EZ) = 1;
			evi((int) HX) = 0;
			evi((int) HY) = 0;
			evi((int) HZ) = 0;
			break;
		}
		break;
	}
	if((int) FORW != 0) slabmodeerror("PropDir <-> int conv. corrupted");
        if((int) BACK != 1) slabmodeerror("PropDir <-> int conv. corrupted");
	cbeta = Cvector(2);
	switch(ppt)
	{
	case PROPAG:  
		cbeta(FORW) = Complex( beta,  0.0); 
		cbeta(BACK) = Complex(-beta,  0.0); 
		img((int) EX) = 0;
		img((int) EY) = 0;
		img((int) EZ) = 1;
		img((int) HX) = 0;
		img((int) HY) = 0;
		img((int) HZ) = 1;
		break;
	case EVANESC: 
		cbeta(FORW) = Complex(0.0, -beta); 
		cbeta(BACK) = Complex(0.0,  beta); 
		img((int) EX) = 0;
		img((int) EY) = 0;
		img((int) EZ) = 0;
		img((int) HX) = 1;
		img((int) HY) = 1;
		img((int) HZ) = 1;
		break;
	}
	brv((int) EX) = 0;
	brv((int) EY) = 0;
	brv((int) EZ) = 1;
	brv((int) HX) = 1;
	brv((int) HY) = 1;
	brv((int) HZ) = 0;
	return;
}

/* integrate a fieldproduct of two modes over Interval i, real version */
double integrate(const Mode& sm1, Fcomp cp1,
                 const Mode& sm2, Fcomp cp2, Interval i)
{
	return sm1.integrate(cp1, sm2, cp2, i);
}

/* ... integrate directional versions of a product of two mode profiles; 
   cc1, cc2 == 1: take the complex conjugate of the component */ 
Complex Mode::integrate(Propdir d1, Fcomp c1, int cc1,
        const Mode& m2, Propdir d2, Fcomp c2, int cc2, Interval i) const
{
	if(vM) slabmodeerror("Mode::integrate (x): vM not implemented");
	if(i.x1 <= i.x0) return 0.0;

	Complex s;
	s = Complex(integrate(c1, m2, c2, i));

	if(d1 == BACK && brv(c1) != 0) s.neg();
	if(d2 == BACK && brv(c2) != 0) s.neg();
	if(img(c1))
	{
		if(cc1) s.mulmi();
		else    s.muli();
	}
	if(m2.img(c2))
	{
		if(cc2) s.mulmi();
		else    s.muli();
	}
	return s;
}

/* ... integrate directional versions of a product of two mode profiles; 
   cc1, cc2 == 1: take the complex conjugate of the component */ 
Complex integrate(const Mode& m1, Propdir d1, Fcomp c1, int cc1,
                  const Mode& m2, Propdir d2, Fcomp c2, int cc2,
                  Interval i)
{
	if(m1.vM) slabmodeerror("Mode::integrate (y): vM not implemented");
	return m1.integrate(d1, c1, cc1, m2, d2, c2, cc2, i); 
}

/* for the QUEP implementation: 
   overlaps of perpendiculary propagating modes & helper functions */
Complex Mode::cross_int(Fcomp cp, int l,
                        Complex ceexp, double xoffs, Interval i) const
{
	if(i.x1<i.x0) return CC0;	
	Complex s;
	double cr, ci;
	cr = ceexp.re;
	ci = ceexp.im;
	s = CC0;
	if(fabs(cr)>1.0e-10) 
	{
		s = (sfd(cp)) ? piece[l].crintfCiS(cr, xoffs, i)
		              : piece[l].crintdCiS(cr, xoffs, i);
		if(dnq(cp)) s *= piece[l].inve; 
	}	
	else 
	{
		s = (sfd(cp)) ? piece[l].crintfExp(ci, xoffs, i)
		              : piece[l].crintdExp(ci, xoffs, i);
		if(dnq(cp)) s *= piece[l].inve; 
	}	
	return s;
}

Complex Mode::cross_int(Fcomp cp, 
                        Complex ceexp, double xoffs, Interval i) const
{
	Complex s;
	double x, xt;
	int l;
	Interval ip;
	if(i.x0 > i.x1) slabmodeerror("cross_int: i.x0 > i.x1"); 
	s = CC0;
	x = i.x0;
	while(fabs(i.x1-x) > S_HDIST)
	{
		l = wg.layeridx(x+S_HDIST/2.0);
		xt = wg.layer(l).x1;
		if(i.x1<xt) xt = i.x1;
		ip.x0 = x; ip.x1 = xt;
		s += cross_int(cp, l, ceexp, xoffs, ip);
		x = xt;
	}
	return s;
}

/* for test purposes: numerical evaluation of integrals */
Complex crossoverlap_num(const Mode& mh, Propdir hdir,
                         const Mode& mv, Propdir vdir,
                         double zpos, double xoffs, Interval i)
{
	if(mh.bdt_b != OPEN && i.x0 < mh.bp_b)
			slabmodeerror("co_num: bounds, i.x0");
	if(mh.bdt_t != OPEN && i.x1 > mh.bp_t)
			slabmodeerror("co_num: bounds, i.x1");
	
	double dx, x, xa, xb;	
	int nx = 10000;
	dx = i.len()/((double) nx);
	Complex sum = CC0;
	Complex f, fh, fv;
	for(int xi=0; xi<=nx-1; ++xi)
	{
		xa = i.x0+((double) xi)*dx;
		xb = i.x0+((double) (xi+1))*dx;
		x = 0.5*(xa+xb);
		f = CC0;	
		fh = (mh.cfield(EX, hdir, x)).conj();
		fv = mv.vertfield(HY, vdir, zpos)*mv.efac(vdir, x-xoffs);
		f = f + fh*fv;
		fh = (mh.cfield(EY, hdir, x)).conj();
		fv = mv.vertfield(HX, vdir, zpos)*mv.efac(vdir, x-xoffs);
		f = f - fh*fv;
		fh = (mh.cfield(HY, hdir, x)).conj();
		fv = mv.vertfield(EX, vdir, zpos)*mv.efac(vdir, x-xoffs);
		f = f + fh*fv;
		fh = (mh.cfield(HX, hdir, x)).conj();
		fv = mv.vertfield(EY, vdir, zpos)*mv.efac(vdir, x-xoffs);
		f = f - fh*fv;
		sum = sum + f*dx;
	}
	return 0.25*sum;
}

Complex Mode::crossoverlap(Propdir hdir, const Mode& mv, Propdir vdir,
                           double zpos, Interval i)
{
	if(vM) return cVcrossoverlap(hdir, mv, vdir, zpos, i);
	Complex sum = CC0;
	Complex ceexp, fcf, cint;
	double xoffs; 
	if(bdt_b != OPEN && i.x0 < bp_b) i.x0 = bp_b;
	if(bdt_t != OPEN && i.x1 > bp_t) i.x1 = bp_t;
	if(vdir == FORW) xoffs = i.x0;
	else             xoffs = i.x1;
	ceexp = mv.cbeta(vdir); 
	if(zero(EX)==0 && mv.zero(HY)==0)
	{
		fcf = mv.vertfield(HY, vdir, zpos);
		cint = fac(EX)*cross_int(EX, ceexp, xoffs, i);
		if(img(EX)) cint.mulmi();
		if(hdir == BACK && brv(EX) != 0) cint.neg();
		if(evi(EX)) cint.neg();
		sum += fcf*cint;
	}
	if(zero(EY)==0 && mv.zero(HZ)==0)
	{
		fcf = mv.vertfield(HX, vdir, zpos);
		cint = fac(EY)*cross_int(EY, ceexp, xoffs, i);
		if(img(EY)) cint.mulmi();
		if(hdir == BACK && brv(EY) != 0) cint.neg();
		if(evi(EY)) cint.neg();
		sum -= fcf*cint;
	}
	if(zero(HY)==0 && mv.zero(EZ)==0)
	{
		fcf = mv.vertfield(EX, vdir, zpos);
		cint = fac(HY)*cross_int(HY, ceexp, xoffs, i);
		if(img(HY)) cint.mulmi();
		if(hdir == BACK && brv(HY) != 0) cint.neg();
		if(evi(HY)) cint.neg();
		sum += fcf*cint;
	}
	if(zero(HX)==0 && mv.zero(EY)==0)
	{
		fcf = mv.vertfield(EY, vdir, zpos);
		cint = fac(HX)*cross_int(HX, ceexp, xoffs, i);
		if(img(HX)) cint.mulmi();
		if(hdir == BACK && brv(HX) != 0) cint.neg();
		if(evi(HX)) cint.neg();
		sum -= fcf*cint;
	}
	sum *= 0.25;
	return sum;
}   

void Mode::crossoverlap(const Mode& mv, double zpos, Interval i,
                        Complex& cff, Complex& cfb, 
                        Complex& cbf, Complex& cbb)
{
	if(vM)
	{
		cVcrossoverlap(mv, zpos, i, cff, cfb, cbf, cbb);
		return;
	}	
	cff = cfb = cbf = cbb = CC0;
	Complex ceexp_f, fcf_f, cint_f;
	Complex ceexp_b, fcf_b, cint_b;
	if(bdt_b != OPEN && i.x0 < bp_b) i.x0 = bp_b;
	if(bdt_t != OPEN && i.x1 > bp_t) i.x1 = bp_t;
	ceexp_f = mv.cbeta(FORW); 
	ceexp_b = mv.cbeta(BACK); 
	if(zero(EX)==0 && mv.zero(HY)==0)
	{
		  fcf_f = mv.vertfield(HY, FORW, zpos);
		 cint_f = cross_int(EX, ceexp_f, i.x0, i);
		  fcf_b = mv.vertfield(HY, BACK, zpos);
		 cint_b = cross_int(EX, ceexp_b, i.x1, i);
		if(img(EX)) { cint_f.mulmi(); cint_b.mulmi(); }
		if(evi(EX)) { cint_f.neg();   cint_b.neg(); }
		cff += fac(EX)*fcf_f*cint_f;
		cfb += fac(EX)*fcf_b*cint_b;
		if(brv(EX)) { cint_f.neg(); cint_b.neg(); }
		cbf += fac(EX)*fcf_f*cint_f;
		cbb += fac(EX)*fcf_b*cint_b;
	}
	if(zero(EY)==0 && mv.zero(HZ)==0)
	{
		  fcf_f = mv.vertfield(HX, FORW, zpos);
		 cint_f = cross_int(EY, ceexp_f, i.x0, i);
		  fcf_b = mv.vertfield(HX, BACK, zpos);
		 cint_b = cross_int(EY, ceexp_b, i.x1, i);
		if(img(EY)) { cint_f.mulmi(); cint_b.mulmi(); }
		if(evi(EY)) { cint_f.neg();   cint_b.neg(); }
		cff -= fac(EY)*fcf_f*cint_f;
		cfb -= fac(EY)*fcf_b*cint_b;
		if(brv(EY)) { cint_f.neg(); cint_b.neg(); }
		cbf -= fac(EY)*fcf_f*cint_f;
		cbb -= fac(EY)*fcf_b*cint_b;
	}
	if(zero(HY)==0 && mv.zero(EZ)==0)
	{
		  fcf_f = mv.vertfield(EX, FORW, zpos);
		 cint_f = cross_int(HY, ceexp_f, i.x0, i);
		  fcf_b = mv.vertfield(EX, BACK, zpos);
		 cint_b = cross_int(HY, ceexp_b, i.x1, i);
		if(img(HY)) { cint_f.mulmi(); cint_b.mulmi(); }
		if(evi(HY)) { cint_f.neg();   cint_b.neg(); }
		cff += fac(HY)*fcf_f*cint_f;
		cfb += fac(HY)*fcf_b*cint_b;
		if(brv(HY)) { cint_f.neg(); cint_b.neg(); }
		cbf += fac(HY)*fcf_f*cint_f;
		cbb += fac(HY)*fcf_b*cint_b;
	}
	if(zero(HX)==0 && mv.zero(EY)==0)
	{
		  fcf_f = mv.vertfield(EY, FORW, zpos);
		 cint_f = cross_int(HX, ceexp_f, i.x0, i);
		  fcf_b = mv.vertfield(EY, BACK, zpos);
		 cint_b = cross_int(HX, ceexp_b, i.x1, i);
		if(img(HX)) { cint_f.mulmi(); cint_b.mulmi(); }
		if(evi(HX)) { cint_f.neg();   cint_b.neg(); }
		cff -= fac(HX)*fcf_f*cint_f;
		cfb -= fac(HX)*fcf_b*cint_b;
		if(brv(HX)) { cint_f.neg(); cint_b.neg(); }
		cbf -= fac(HX)*fcf_f*cint_f;
		cbb -= fac(HX)*fcf_b*cint_b;
	}
	cff *= 0.25;
	cfb *= 0.25;
	cbf *= 0.25;
	cbb *= 0.25;
	return; 
}   

/* product 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of a mode with a Gaussian beam gb, numerical evaluation of the integrals,
   alternatively forward or backward propagating fields are considered */
#define CNINT_VAN 1.0e-7
#define CNINT_NUM 10 
Complex Mode::overlap(const Gaussianbeam& gb, Propdir pd) const
{
	if(pol != gb.pol) return CC0; 
	Complex s = CC0;
	Interval lyr;
	double bp;
	for(int l=0; l<=wg.nx+1; ++l)
	{
		lyr.x0 = piece[l].bif;
		lyr.x1 = piece[l].tif;
		if(l == 0 && bdt_b == OPEN)
		{ 
			bp = lyr.x1+log(CNINT_VAN)/(piece[l]).g;
			if(bp > lyr.x0) lyr.x0 = bp;
			bp = gb.x0-sqrt(-log(CNINT_VAN)*gb.w*gb.w);
			if(bp > lyr.x0) lyr.x0 = bp;
		}
		if(l == wg.nx+1 && bdt_t == OPEN)
		{ 
			bp = lyr.x0-log(CNINT_VAN)/(piece[l]).g;
			if(bp < lyr.x1) lyr.x1 = bp;
			bp = gb.x0+sqrt(-log(CNINT_VAN)*gb.w*gb.w);
			if(bp < lyr.x1) lyr.x1 = bp;
		}
		if(lyr.x0 < lyr.x1)
		{		
			double e = wg.eps(l); if(e < 1.0) e = 1.0;
			int ndp = (int) ceil(((double) CNINT_NUM)*lyr.len()
			                     /wg.lambda*sqrt(e));
			int nwl = 1;
			if((piece[l]).bhv == OSC) 
				nwl = (int) ceil((piece[l]).g*lyr.len()/2.0/PI);
			if(CNINT_NUM*nwl > ndp) ndp = CNINT_NUM*nwl;
			if(pol == TE)
			{
				s -= cnintegrate(EY, gb, HX, lyr, ndp);
				s -= cnintegrate(HX, gb, EY, lyr, ndp);
			}
			else 
			{
				s += cnintegrate(EX, gb, HY, lyr, ndp);
				s += cnintegrate(HY, gb, EX, lyr, ndp);
			}
		}
	}
	if(pd == BACK) s.neg(); 
	return s*0.25;
} 

/* ... helper function: numerical integration by Gaussian quadrature rules
   based on: Numerical Recipes in C --- The Art of Scientific Computing
   Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */ 
Complex Mode::cnintegrate(Fcomp cpm, const Gaussianbeam& gb, Fcomp cpb,
                         Interval lyr, int numx) const
{
	double xr, xm, dx, x;
	static double xc[] = {0.0, 0.1488743389, 0.4333953941,
		              0.6794095682, 0.8650633666, 0.9739065285};
	static double wc[] = {0.0, 0.2955242247, 0.2692667193,
		              0.2190863625, 0.1494513491, 0.0666713443};
	Complex sum = CC0; 
	double a, b, step;
	step = lyr.len()/((double) numx);
	Complex s, mf, bf, fval; 
	for(int i=0; i<=numx-1; ++i)
	{
		a  = lyr.x0+((double) i)*step;
		b  = a+step;
		xm = 0.5*(b+a);
		xr = 0.5*(b-a);
		s  = CC0;
		for(int j=1; j<=5; j++) 
		{
			dx = xr*xc[j];
			   x  = xm+dx;
			fval  = (cfield(cpm, FORW, x).conj())*gb.field(cpb, x);
			   x  = xm-dx;
			fval += (cfield(cpm, FORW, x).conj())*gb.field(cpb, x);
			s += wc[j]*fval;
		}
		sum += s*xr;
	}
	return sum;
}

/* for the vQUEP routines: extension to vTE, vTM modes */ 

/* prepare vTE / vTM mode representation, for given ky */
void Mode::vectorify(Complex val_ky)
{
// fprintf(stderr, "%s: %g", ids, power());
	vM = 1;
	ky = val_ky;
	Complex kzq = betaq-ky*ky;
        double r, phi;
        r = sqrt(kzq.abs());
        phi = kzq.args()/2.0;
        kz = expi(phi)*r;
	if(fabs(kz.re) > r*(SELECT_K_COMP)) 
	{
		if(kz.re < 0.0) kz.neg();
	}
	else 
	{
		if(kz.im > 0.0) kz.neg();
	}
	if(fabs(kz.re/r) < SELECT_K_COMP) kz.re = 0.0;
	if(fabs(kz.im/r) < SELECT_K_COMP) kz.im = 0.0;
	if((int) FORW != 0) slabmodeerror("PropDir <-> int conv. corrupted");
        if((int) BACK != 1) slabmodeerror("PropDir <-> int conv. corrupted");
	ckz = Cvector(2);
	ckz(FORW) =  kz;
	ckz(BACK) = -kz;
	if(fabs(kz.re) >= fabs(kz.im)) cVppt = PROPAG;
	else                           cVppt = EVANESC;
	
	if((int) EX != 0) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
        if((int) EY != 1) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
        if((int) EZ != 2) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
        if((int) HX != 3) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
        if((int) HY != 4) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
        if((int) HZ != 5) slabmodeerror("vectorify: Fcomp <-> int conv. corrupted");
	cVzero = Ivector(6);
	cVsfd  = Ivector(6);
	cVdnq  = Ivector(6);
	cVfac  = Cvector(6);
	cVbrv  = Ivector(6);
	switch(pol)
	{
	case TE:	
		cVzero((int) EX) = 1; 
		cVzero((int) EY) = 0; 
		cVzero((int) EZ) = 0; 
		cVzero((int) HX) = 0; 
		cVzero((int) HY) = 0; 
		cVzero((int) HZ) = 0; 

		cVsfd((int) EX) = 1;
		cVsfd((int) EY) = 1;
		cVsfd((int) EZ) = 1;
		cVsfd((int) HX) = 1;
		cVsfd((int) HY) = 0;
		cVsfd((int) HZ) = 0;

		cVdnq((int) EX) = 0;
		cVdnq((int) EY) = 0;
		cVdnq((int) EZ) = 0;
		cVdnq((int) HX) = 0;
		cVdnq((int) HY) = 0;
		cVdnq((int) HZ) = 0;

		cVfac((int) EX) = 0.0;
		cVfac((int) EY) = ckz(FORW)/betaq; 
		cVfac((int) EZ) = -ky/betaq;
		cVfac((int) HX) = -invommu0; 
		cVfac((int) HY) = CCI*invommu0*ky/betaq; 
		cVfac((int) HZ) = CCI*invommu0*ckz(FORW)/betaq; 

		cVbrv((int) EX) = 0;
		cVbrv((int) EY) = 1;
		cVbrv((int) EZ) = 0;
		cVbrv((int) HX) = 0;
		cVbrv((int) HY) = 0;
		cVbrv((int) HZ) = 1;

		break;

	case TM:
		cVzero((int) EX) = 0; 
		cVzero((int) EY) = 0; 
		cVzero((int) EZ) = 0; 
		cVzero((int) HX) = 1; 
		cVzero((int) HY) = 0; 
		cVzero((int) HZ) = 0; 

		cVsfd((int) EX) = 1;
		cVsfd((int) EY) = 0;
		cVsfd((int) EZ) = 0;
		cVsfd((int) HX) = 1;
		cVsfd((int) HY) = 1;
		cVsfd((int) HZ) = 1;

		cVdnq((int) EX) = 1;
		cVdnq((int) EY) = 1;
		cVdnq((int) EZ) = 1;
		cVdnq((int) HX) = 0;
		cVdnq((int) HY) = 0;
		cVdnq((int) HZ) = 0;

		cVfac((int) EX) = invomep0; 
		cVfac((int) EY) = CCmI*invomep0*ky/betaq; 
		cVfac((int) EZ) = CCmI*invomep0*ckz(FORW)/betaq; 
		cVfac((int) HX) = 0.0; 
		cVfac((int) HY) = ckz(FORW)/betaq;
		cVfac((int) HZ) = -ky/betaq; 

		cVbrv((int) EX) = 0;
		cVbrv((int) EY) = 0;
		cVbrv((int) EZ) = 1;
		cVbrv((int) HX) = 0;
		cVbrv((int) HY) = 1;
		cVbrv((int) HZ) = 0;

		break;
	}
// fprintf(stderr, "-> %g", power());
	normalize();
// fprintf(stderr, " -> %g\n", power());
	return;
}

/* vTE/vTM modes, complex mode profile
   values of component cp: EX - HZ at transverse coordinate x,
   for propagation direction d: FORW, BACK */
Complex Mode::cVfield(Fcomp cp, Propdir d, double x) const
{
	if(cp > HZ) slabmodeerror("cVfield: cp");
	if(cVzero(cp)) return CC0;
	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);
	double f = (cVsfd(cp)) ? piece[l].bfld(x) : piece[l].bder(x); 
	f *= beta;
	if(cVdnq(cp)) f *= piece[l].inve;
	if(d == FORW) return cVfac(cp)*f;
	if(cVbrv(cp)) f = -f;
	return cVfac(cp)*f;
}

/* if viewed as a mode propagating along x: complex vert. mode profiles */
Complex Mode::vertVfield(Fcomp cp, Propdir d, double z) const
{
	Complex f;
	switch(cp)
	{
	case EX: f = cVfield(EZ, d, z); f.neg(); return f;
	         break;
	case EY: f = cVfield(EY, d, z); f.neg(); return f;
	         break;
	case EZ: f = cVfield(EX, d, z); f.neg(); return f;
	         break;
	case HX: return cVfield(HZ, d, z); 
	         break;
	case HY: return cVfield(HY, d, z); 
	         break;
	case HZ: return cVfield(HX, d, z); 
	         break;
	default: slabmodeerror("vertVfield: cp");
	         break;
	}
	return CC0;
}

/* for the vBEP / vQUEP implementation: products 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, combinations of forward and backward propagating versions */
void cVbidir_overlaps(const Mode& m1, const Mode& m2, 
                      Complex& off, Complex& ofb, 
                      Complex& obf, Complex& obb)
{
	Complex ex1hy2, ey1hx2, hy1ex2, hx1ey2, eh;
	Interval i(-AWAY, AWAY);
	ex1hy2 = m1.cVintegrate(EX, m2, HY, i);
	ey1hx2 = m1.cVintegrate(EY, m2, HX, i);
	hy1ex2 = m1.cVintegrate(HY, m2, EX, i);
	hx1ey2 = m1.cVintegrate(HX, m2, EY, i);

	off = CC0;
	ofb = CC0;
	obf = CC0;
	obb = CC0;

	eh = ex1hy2;                                                       off += eh;
	eh = ex1hy2;                            if(m2.cVbrv(HY)) eh.neg(); ofb += eh;
	eh = ex1hy2; if(m1.cVbrv(EX)) eh.neg();                            obf += eh;
	eh = ex1hy2; if(m1.cVbrv(EX)) eh.neg(); if(m2.cVbrv(HY)) eh.neg(); obb += eh;
	
	eh = ey1hx2;                                                       off -= eh;
	eh = ey1hx2;                            if(m2.cVbrv(HX)) eh.neg(); ofb -= eh;
	eh = ey1hx2; if(m1.cVbrv(EY)) eh.neg();                            obf -= eh;
	eh = ey1hx2; if(m1.cVbrv(EY)) eh.neg(); if(m2.cVbrv(HX)) eh.neg(); obb -= eh;
	
	eh = hy1ex2;                                                       off += eh;
	eh = hy1ex2;                            if(m2.cVbrv(EX)) eh.neg(); ofb += eh;
	eh = hy1ex2; if(m1.cVbrv(HY)) eh.neg();                            obf += eh;
	eh = hy1ex2; if(m1.cVbrv(HY)) eh.neg(); if(m2.cVbrv(EX)) eh.neg(); obb += eh;
	
	eh = hx1ey2;                                                       off -= eh;
	eh = hx1ey2;                            if(m2.cVbrv(EY)) eh.neg(); ofb -= eh;
	eh = hx1ey2; if(m1.cVbrv(HX)) eh.neg();                            obf -= eh;
	eh = hx1ey2; if(m1.cVbrv(HX)) eh.neg(); if(m2.cVbrv(EY)) eh.neg(); obb -= eh;

	off *= 0.25;
	ofb *= 0.25;
	obf *= 0.25;
	obb *= 0.25;
	
	return;
} 

/* for the vBEP / vQUEP implementation: products 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, combinations of forward and backward propagating versions */
Complex cVoverlap(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2)
{
	return m1.cVoverlap(d1, m2, d2);
} 

/* integrate a fieldproduct of two vTE/ vTM modes (FORW) over Interval i, 
   complex conjugation of the first factors */
Complex Mode::cVintegrate(Fcomp cp, const Mode& m, Fcomp cpm, Interval i) const
{
	if(((int) cp) > ((int) HZ)) slabmodeerror("cVintegrate: cp");
	if(  cVzero(cp) ) return CC0;
	if(m.cVzero(cpm)) return CC0;

	double x, xt, xt1, xt2;
	double 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 = 0.0;
	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;
		if(cVsfd(cp))
		{
			z = (m.cVsfd(cpm)) ? piece[l].intff(m.piece[lm], iv)
			                   : piece[l].intfd(m.piece[lm], iv);
		}
		else
		{
			z = (m.cVsfd(cpm)) ? piece[l].intdf(m.piece[lm], iv)
			                   : piece[l].intdd(m.piece[lm], iv);
		}
		if(  cVdnq(cp )) z *=    piece[l].inve;
		if(m.cVdnq(cpm)) z *= m.piece[lm].inve;
		s += z;
		x = xt;
	}
	s *= beta*m.beta;
	return (cVfac(cp).conj())*(m.cVfac(cpm))*s;
/*
	Complex anr = (cVfac(cp).conj())*(m.cVfac(cpm))*s;
if(anr.abs() > 1.0e-6)
{
	Complex nmr = num_cVintegrate(cp, m, cpm, i);
fprintf(stderr, "\n int %g -> %g\n", i.x0, i.x1);
fprintf(stderr, "%s betaq = %g %c%c\n", ids, betaq, fldchr(cp), cpchr(cp));
fprintf(stderr, "%s betaq = %g %c%c\n", m.ids, m.betaq, fldchr(cpm), cpchr(cpm));
fprintf(stderr, "anr = %g +i %g\n", anr.re, anr.im);
fprintf(stderr, "nmr = %g +i %g\n", nmr.re, nmr.im);
}
	return anr;
*/
}

/* check: numerical evaluation of the integral */
Complex Mode::num_cVintegrate(Fcomp cp, const Mode& m, Fcomp cpm, Interval i) const
{
	double xl, xr, dx;
	Complex fl, fr;
	int Ns = 10000;
	Interval ni = i;
	if(ni.x0 < -10.0) ni.x0 = -10.0;
	if(ni.x1 >  10.0) ni.x1 =  10.0;
	dx = (ni.x1-ni.x0)/((double) Ns);
	Complex sum = CC0;
	for(int s=0; s<=Ns-1; ++s)
	{
		xl=ni.x0+((double) s)*dx;
		xr=xl+dx;
		fl = (cVfield(cp, FORW, xl).conj())*m.cVfield(cpm, FORW, xl);
		fr = (cVfield(cp, FORW, xr).conj())*m.cVfield(cpm, FORW, xr);
		sum += (fl+fr)*0.5;
	}
	sum = sum*dx;
	return sum;
}

/* for the vQUEP implementation: 
   overlaps of perpendiculary propagating vTE and vTM modes 
   & helper functions */
Complex Mode::cVcross_int(Fcomp cp, 
                          Complex ceexp, double xoffs, Interval i) const
{
	Complex s;
	double x, xt;
	int l;
	Interval ip;
	if(i.x0 > i.x1) slabmodeerror("cVcross_int: i.x0 > i.x1"); 
	s = CC0;
	x = i.x0;
	while(fabs(i.x1-x) > S_HDIST)
	{
		l = wg.layeridx(x+S_HDIST/2.0);
		xt = wg.layer(l).x1;
		if(i.x1<xt) xt = i.x1;
		ip.x0 = x; ip.x1 = xt;
		if(ip.x1 > ip.x0) 
		{
			Complex sl;
			sl = (cVsfd(cp)) ? piece[l].crintfcE(ceexp, xoffs, ip)
			                 : piece[l].crintdcE(ceexp, xoffs, ip);
			if(cVdnq(cp)) sl *= piece[l].inve; 
			s += sl;
		}
		x = xt;
	}
	s *= beta;
	return ((cVfac(cp)).conj())*s;
}
Complex Mode::cVcrossoverlap(Propdir hdir, const Mode& mv, Propdir vdir,
                             double zpos, Interval i) const
{
	Complex sum = CC0;
	Complex ceexp, cint;
	double xoffs; 
	if(bdt_b != OPEN && i.x0 < bp_b) i.x0 = bp_b;
	if(bdt_t != OPEN && i.x1 > bp_t) i.x1 = bp_t;
	if(vdir == FORW) xoffs = i.x0;
	else             xoffs = i.x1;
	ceexp = mv.ckz(vdir); 

	if(cVzero(EX)==0 && mv.cVzero(HY)==0)
	{
		cint = cVcross_int(EX, ceexp, xoffs, i);
		if(hdir == BACK && cVbrv(EX) != 0) cint.neg();
		sum += cint*mv.vertVfield(HY, vdir, zpos);
	}
	if(cVzero(EY)==0 && mv.cVzero(HZ)==0)
	{
		cint = cVcross_int(EY, ceexp, xoffs, i);
		if(hdir == BACK && cVbrv(EY) != 0) cint.neg();
		sum -= cint*mv.vertVfield(HX, vdir, zpos);
	}
	if(cVzero(HY)==0 && mv.cVzero(EZ)==0)
	{
		cint = cVcross_int(HY, ceexp, xoffs, i);
		if(hdir == BACK && cVbrv(HY) != 0) cint.neg();
		sum += cint*mv.vertVfield(EX, vdir, zpos);
	}
	if(cVzero(HX)==0 && mv.cVzero(EY)==0)
	{
		cint = cVcross_int(HX, ceexp, xoffs, i);
		if(hdir == BACK && cVbrv(HX) != 0) cint.neg();
		sum -= cint*mv.vertVfield(EY, vdir, zpos);
	}
	sum *= 0.25;

	return sum;
}   

void Mode::cVcrossoverlap(const Mode& mv, double zpos, Interval i,
                        Complex& cff, Complex& cfb, 
                        Complex& cbf, Complex& cbb) const
{
	cff = cfb = cbf = cbb = CC0;
	Complex ceexp_f, ceexp_b, fcf_f, fcf_b, cint_f, cint_b;
	if(bdt_b != OPEN && i.x0 < bp_b) i.x0 = bp_b;
	if(bdt_t != OPEN && i.x1 > bp_t) i.x1 = bp_t;
	ceexp_f = mv.ckz(FORW); 
	ceexp_b = mv.ckz(BACK); 
	if(cVzero(EX)==0 && mv.cVzero(HY)==0)
	{
		cint_f = cVcross_int(EX, ceexp_f, i.x0, i);
		cint_b = cVcross_int(EX, ceexp_b, i.x1, i);
		fcf_f = mv.vertVfield(HY, FORW, zpos);
		fcf_b = mv.vertVfield(HY, BACK, zpos);
		cff += cint_f*fcf_f;
		cfb += cint_b*fcf_b;
		if(cVbrv(EX) != 0) { cint_f.neg(); cint_b.neg(); }
		cbf += cint_f*fcf_f;
		cbb += cint_b*fcf_b;
	}
	if(cVzero(EY)==0 && mv.cVzero(HZ)==0)
	{
		cint_f = cVcross_int(EY, ceexp_f, i.x0, i);
		cint_b = cVcross_int(EY, ceexp_b, i.x1, i);
		fcf_f = mv.vertVfield(HX, FORW, zpos);
		fcf_b = mv.vertVfield(HX, BACK, zpos);
		cff -= cint_f*fcf_f;
		cfb -= cint_b*fcf_b;
		if(cVbrv(EY) != 0) { cint_f.neg(); cint_b.neg(); }
		cbf -= cint_f*fcf_f;
		cbb -= cint_b*fcf_b;
	}
	if(cVzero(HY)==0 && mv.cVzero(EZ)==0)
	{
		cint_f = cVcross_int(HY, ceexp_f, i.x0, i);
		cint_b = cVcross_int(HY, ceexp_b, i.x1, i);
		fcf_f = mv.vertVfield(EX, FORW, zpos);
		fcf_b = mv.vertVfield(EX, BACK, zpos);
		cff += cint_f*fcf_f;
		cfb += cint_b*fcf_b;
		if(cVbrv(HY) != 0) { cint_f.neg(); cint_b.neg(); }
		cbf += cint_f*fcf_f;
		cbb += cint_b*fcf_b;
	}
	if(cVzero(HX)==0 && mv.cVzero(EY)==0)
	{
		cint_f = cVcross_int(HX, ceexp_f, i.x0, i);
		cint_b = cVcross_int(HX, ceexp_b, i.x1, i);
		fcf_f = mv.vertVfield(EY, FORW, zpos);
		fcf_b = mv.vertVfield(EY, BACK, zpos);
		cff -= cint_f*fcf_f;
		cfb -= cint_b*fcf_b;
		if(cVbrv(HX) != 0) { cint_f.neg(); cint_b.neg(); }
		cbf -= cint_f*fcf_f;
		cbb -= cint_b*fcf_b;
	}
	cff *= 0.25;
	cfb *= 0.25;
	cbf *= 0.25;
	cbb *= 0.25;
	return;
}