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

/*
 * veims.cpp
 * Variational Effective Index Mode Solver VEIMS
 * mode analysis for waveguides with 2D rectangular cross section
 * --- the simplest version, separabale mode profile approximations
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"matlvis.h"
#include"slamode.h"
#include"integral.h"
#include"slamarr.h"
#include"slams.h"
#include"veims.h"

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

/* VEIMS mode profile representation: separable fields */
void EIMode::fldrep(Fcomp cp, double &fac, FldorDer &vfod, FldorDer &hfod, double x, double y) const
{
	double iomu0, ioep0, epsr, epseff, epsc, q;
	int rly = st(rsegi).layeridx(x);
	int seg = st.segidx(y);
	if(pol == TE)
	{
		switch(cp)
		{
		case EX: 
			vfod = FLD;
			hfod = FLD;
			fac = 0.0;
			break;
		case EY: 
			epseff = hwg.eps(seg);
			vfod = FLD;
			hfod = FLD;
			fac = beta*vmp.beta/k0/k0/epseff;
			break;
		case EZ:
			epseff = hwg.eps(seg);
			vfod = FLD;
			hfod = DER;
			fac = -vmp.beta/k0/k0/epseff;
			break;
		case HX:	
			iomu0 = vmp.invommu0; 
			vfod = FLD;
			hfod = FLD;
			fac = -vmp.beta*iomu0;
			break;
		case HY:	
			iomu0 = vmp.invommu0; 
			epseff = hwg.eps(seg);
			vfod = DER;
			hfod = DER;
			fac = -vmp.beta*iomu0/k0/k0/epseff;
			break;
		case HZ:	
			iomu0 = vmp.invommu0; 
			epseff = hwg.eps(seg);
			vfod = DER;
			hfod = FLD;
			fac = vmp.beta*beta*iomu0/k0/k0/epseff;
			break;
		default: 
			eimodeerror("setfldrep: illegal cp");	
			break;
		}
		return;
	}
	switch(cp)
	{
	case EX:	
		ioep0 = vmp.invomep0; 
		epsr = st(rsegi).eps(rly);
		vfod = FLD;
		hfod = FLD;
		fac = vmp.beta*ioep0/epsr;
		break;
	case EY:	
		ioep0 = vmp.invomep0; 
		epsr = st(rsegi).eps(rly);
		epsc = hwg_ec(seg);
		q = hwg_q(seg);
		vfod = DER;
		hfod = DER;
		fac = vmp.beta*ioep0/epsr/k0/k0/epsc*q;
		break;
	case EZ:	
		ioep0 = vmp.invomep0; 
		epsr = st(rsegi).eps(rly);
		epsc = hwg_ec(seg);
		q = hwg_q(seg);
		vfod = DER;
		hfod = FLD;
		fac = -vmp.beta*beta*ioep0/epsr/k0/k0/epsc*q;
		break;
	case HX: 
		vfod = FLD;
		hfod = FLD;
		fac = 0.0;
		break;
	case HY: 
		epsc = hwg_ec(seg);
		vfod = FLD;
		hfod = FLD;
		fac = beta*vmp.beta/k0/k0/epsc;
		break;
	case HZ:
		epsc = hwg_ec(seg);
		vfod = FLD;
		hfod = DER;
		fac = -vmp.beta/k0/k0/epsc;
		break;
	default: 
		eimodeerror("setfldrep: illegal cp");	
		break;
	}
	return;
}
void EIMode::fldrep(Fcomp cp, double &fac, FldorDer &vfod, FldorDer &hfod, int l, int m) const
{
	Rect r = wg.rectbounds(l, m);
	fldrep(cp, fac, vfod, hfod, 0.5*(r.x0+r.x1), 0.5*(r.y0+r.y1));
	return;
}

/* mode profile, field values at position (x, y) 
   cp: EX - HZ, SZ */
double EIMode::field(Fcomp cp, double x, double y) const
{
	double vf, hf, fac;
	FldorDer hfod, vfod;
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ:	
		fldrep(cp, fac, vfod, hfod, x, y);
		vf = ((vfod == FLD) ? vmp.field(x) : vmp.d_field(x));
		hf = ((hfod == FLD) ? hmp.field(y) : hmp.d_field(y));
		return fac*vf*hf; 
		break;
	case SZ:
		return 0.5*(field(EX, x, y)*field(HY, x, y)
			   -field(HX, x, y)*field(EY, x, y)); 
		break;
	default: 
		return 0.0;
		break;
	}
	return 0.0;
}

/* store num values of component cp between (x0, y0) and (x1, y1) in a vector 
   cp: EX - HZ, SZ */
Dvector EIMode::field(Fcomp cp, int num, 
                      double x0, double y0, double x1, double y1) const
{
	double x, y;
	double dx, dy;
	int j;
	if(num <= 0) eimodeerror("field: num <= 0");
	Dvector f(num);
	if(num == 1) 
	{
		x = 0.5*(x0+x1);
		y = 0.5*(x0+x1);
		f(0) = field(cp, x, y);
		return f;
	}
	dx = (x1-x0)/(num-1);
	dy = (y1-y0)/(num-1);
	for(j=0; j<=num-1; ++j)
	{
		x = x0+j*dx;	
		y = y0+j*dy;	
		f(j) = field(cp, x, y);
	}
	return f;
}

/* evaluate component cp on a rectangular npx x npy mesh on area disp 
   cp: EX - HZ, SZ 
   foa: ORG, MOD, SQR */
Dmatrix EIMode::field(Rect disp, int npx, int npy, Fcomp cp, Afo foa) const
{
	double x, y;
	double dx, dy;
	int i, j;
	if(npx <= 1) eimodeerror("field: npx <= 1");
	if(npy <= 1) eimodeerror("field: npy <= 1");
	Dmatrix f(npx, npy);
	double ft;
	dx = (disp.x1-disp.x0)/(npx-1);
	dy = (disp.y1-disp.y0)/(npy-1);
	for(i=0; i<=npx-1; ++i)
	{
		x = disp.x0+i*dx;	
		for(j=0; j<=npy-1; ++j)
		{
			y = disp.y0+j*dy;	
			ft = field(cp, x, y);
			switch(foa)
			{
			case MOD: 
				ft = fabs(ft); 
				break;
			case SQR: 
				ft = ft*ft; 
				break;
			case ORG: 
				break;
			case REP: 
				if(cp == EZ || cp == HZ) ft = 0.0;
				break;
			case IMP: 
				if(cp != EZ && cp != HZ) ft = 0.0;
				break;
			}
			f(i, j) = ft;
		}
	}
	return f;
} 

/* --- integrals of mode profiles ------------------------------ */

/* helper function */
double EIMode::prodint(int l, int m, Fcomp cp1, Fcomp cp2, Rect r) const
{
	FldorDer vfod1, vfod2, hfod1, hfod2;
	double fac1, fac2;
	fldrep(cp1, fac1, vfod1, hfod1, l, m);
	fldrep(cp2, fac2, vfod2, hfod2, l, m);
	return fac1*fac2*vmp.integrate(vfod1, vfod2, Interval(r.x0, r.x1))
	                *hmp.integrate(hfod1, hfod2, Interval(r.y0, r.y1));
}

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

/* z component of the Poyntingvector, integrated over the entire x-y-plane */
/* TE part */
double EIMode::tepower() const
{
	return -0.5*recint(HX, EY, XYplane);
}
/* TM part */
double EIMode::tmpower() const
{
	return 0.5*recint(EX, HY, XYplane);
}
/* total power */
double EIMode::power() const
{
	return tepower()+tmpower();
}

/* --- integrals of mode profile product along lines ------------- */

/* integration of field products along horizontal lines */
double EIMode::horlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
		  double x, double ya, double yb) const
{
	FldorDer vfod1, vfod2, hfod1, hfod2;
	double fac1, fac2;
	fldrep(cp1, fac1, vfod1, hfod1, l, m);
	fldrep(cp2, fac2, vfod2, hfod2, l, m);
	return fac1*fac2*((vfod1 == FLD) ? vmp.field(x) : vmp.d_field(x))
	                *((vfod2 == FLD) ? vmp.field(x) : vmp.d_field(x))
	                *hmp.integrate(hfod1, hfod2, Interval(ya, yb));
}

/* integration of field products along vertical lines */
double EIMode::verlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
		  double xa, double xb, double y) const
{
	FldorDer vfod1, vfod2, hfod1, hfod2;
	double fac1, fac2;
	fldrep(cp1, fac1, vfod1, hfod1, l, m);
	fldrep(cp2, fac2, vfod2, hfod2, l, m);
	return fac1*fac2*vmp.integrate(vfod1, vfod2, Interval(xa, xb))
	                *((hfod1 == FLD) ? hmp.field(y) : hmp.d_field(y))
	                *((hfod2 == FLD) ? hmp.field(y) : hmp.d_field(y));
}
 
/* -------------------------------------------------------------------- */

/* set maximum values for E, H, Sz (only a coarse approximation !) */
void EIMode::setfieldmax()
{
	int l, m;
	double x, y, w, s;
	int np = 10;
	Rect r;
	minE = 1.0e+300;
	maxE = -1.0e+300;
	minH = 1.0e+300;
	maxH = -1.0e+300;
	minS = 1.0e+300;
	maxS = -1.0e+300;
	for(l=1; l<=wg.nx; ++l)
	{
		for(m=1; m<=wg.ny;++m)
		{
			r = wg.rectbounds(l, m);
			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(EX, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EY, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EZ, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(EX, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EY, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EZ, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
			}

			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(HX, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HY, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HZ, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(HX, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HY, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HZ, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
			}

			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(SZ, x, y); 
				if(w > maxS) maxS = w;
				if(w < minS) minS = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(SZ, x, y); 
				if(w > maxS) maxS = w;
				if(w < minS) minS = w;
			}
		}
	}
	return;
}

/* coordinates of the maximum mode intensity on rectangle r, 
   average 1/e^2-radius in x- and y-direction 
   only a coarse approximation ! */
double EIMode::maxintensity(Rect r, double& xm, double& ym, 
				  double& rx, double& ry) const
{
	double x, y, w, sx, sy;
	int np = 50;
	double max;
	double a, b;
	max = -1.0e+300;
	xm = 0.0;
	ym = 0.0;

	sx = (r.x1-r.x0)/np;
	sy = (r.y1-r.y0)/np;
	for(x=r.x0; x<=r.x1; x+=sx)
	{
		for(y=r.y0; y<=r.y1; y+=sy)
		{
			w = fabs(field(SZ, x, y)); 
			if(w > max)
			{
				max = w;
				xm = x;
				ym = y;
			}
		}
	}

	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		x -= sx;
		w = fabs(field(SZ, x, y)); 
	}
	a = x;
	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		x += sx;
		w = fabs(field(SZ, x, y)); 
	}
	b = x;
	rx = 0.5*(b-a);

	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		y -= sy;
		w = fabs(field(SZ, x, y)); 
	}
	a = y;
	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		y += sy;
		w = fabs(field(SZ, x, y)); 
	}
	b = y;
	ry = 0.5*(b-a);
	 
	return max;
}

/* normalize mode to power() == nrm */
void EIMode::normalize(double nrm)
{
	double p = power();
	p = 1.0/p;
	hmp.normalize(p);
	setfieldmax();
	return;
}

/* translate mode by (dx, dy) */
void EIMode::translate(double dx, double dy)
{
	wg.translate(dx, dy);
	st.xtranslate(dx);
	st.ytranslate(dy);
	return;
}

/* write profile section along line (x0,y0) -> (x1,y1) 
	cp: EX - HZ, SZ
	ext0, ext1: filename id characters
	np: number of output points
*/
void EIMode::writesec(Fcomp cp, char ext0, char ext1, int np,
	              double x0, double y0, double x1, double y1) const
{
	FILE *dat;
	int i;
	double x, y, dx, dy;
	char name[13] = "t_____.xyf";
	switch(pol)
	{
		case TE: name[1] = 'e'; break;
		case TM: name[1] = 'm'; break;
	}
	name[2] = fldchr(cp);
	name[3] = cpchr(cp);
	name[4] = ext0;
	name[5] = ext1;
	fprintf(stderr, ">> %s\n", name);
	dat = fopen(name, "w+");
	dx = x1-x0; dy = y1-y0;
	for(i=0; i<=np; ++i)
	{
		x = x0+i*dx/np;
		y = y0+i*dy/np;
		if(dx == 0.0) fprintf(dat, "%g  %g\n", y, field(cp, x, y));
		else
		{
			if(dy == 0.0) fprintf(dat, "%g  %g\n", 
						    x, field(cp, x, y));
			else
			{
			   fprintf(dat, "%g  %g\n", i*sqrt(dx*dx+dy*dy)/np,
			   				field(cp, x, y));
			}
		}
	}
	fclose(dat);
	return;
}

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

/* permittivity perturbation: propagation constant shift */
Complex EIMode::phaseshift(ChWgPert p) const
{
	Complex pt;
	double nrm;

	Complex db;
	db = CC0;

	// xx
	pt = p.e(0, 0); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EX, p.r); 	

	// yy
	pt = p.e(1, 1); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EY, p.r); 	

	// zz
	pt = p.e(2, 2); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EZ, EZ, p.r); 	

	// xy, yx
	pt = p.e(0, 1) + p.e(1, 0); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EY, p.r); 	

	// xz, zx
	pt = p.e(0, 2) - p.e(2, 0); 
	pt = CCI*pt;
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EZ, p.r); 	

	// yz, zy
	pt = p.e(1, 2) - p.e(2, 1); 
	pt = CCI*pt;
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EZ, p.r); 	
	
	nrm  = recint(HY, EX, XYplane);
	nrm -= recint(HX, EY, XYplane);
	nrm *= 2.0;
	db = db*(val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm);
	return db; 
}

/* geometry variations: propagation constant shift due to
   moving the horizontal boundary between rectangles [l,m] and [l+1,m]
   ( moving hx(l) for hy(m-1) < y < hy(m) ) */
double EIMode::horgeovar(int l, int m) const
{
	double pt;
	double nrm;
	double y0, y1;
	if(l<0) return 0.0;
	if(m<0) return 0.0;
	if(l>wg.nx+1) return 0.0;
	if(m>wg.ny) return 0.0;
	if(m <= 0) y0 = XYplane.y0;
	else       y0 = wg.hy(m-1);
	if(m >= wg.ny+1) y1 = XYplane.y1;
	else             y1 = wg.hy(m);

	pt  = 0.5*horlineint(l+1, m, EX, EX, wg.hx(l), y0, y1)
	         *wg.eps(l+1,m)/wg.eps(l,m); 
	pt += 0.5*horlineint(l, m,   EX, EX, wg.hx(l), y0, y1)
	         *wg.eps(l,m)/wg.eps(l+1,m); 

	pt += 0.5*(horlineint(l+1, m, EY, EY, wg.hx(l), y0, y1) 
	          +horlineint(l, m,   EY, EY, wg.hx(l), y0, y1)); 

	pt += 0.5*(horlineint(l+1, m, EZ, EZ, wg.hx(l), y0, y1) 
	          +horlineint(l, m,   EZ, EZ, wg.hx(l), y0, y1)); 

	pt *= (wg.eps(l,m)-wg.eps(l+1,m));

	nrm  = recint(HY, EX, XYplane);
	nrm -= recint(HX, EY, XYplane);
	nrm *= 2.0;

	pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm;

	return pt;
}

/* geometry variations: propagation constant shift due to
   moving the vertical boundary between rectangles [l,m] and [l,m+1]
   ( moving hy(m) for hx(l-1) < x < hx(l) ) */
double EIMode::vergeovar(int l, int m) const
{
	double pt;
	double nrm;
	double x0, x1;
	if(m<0) return 0.0;
	if(l<0) return 0.0;
	if(m>wg.ny) return 0.0;
	if(l>wg.nx+1) return 0.0;
	if(l == 0) x0 = XYplane.x0;
	else       x0 = wg.hx(l-1);
	if(l == wg.nx+1) x1 = XYplane.x1;
	else             x1 = wg.hx(l);

	if(pol == TM)
	{
		pt = 0.5*(verlineint(l, m+1, HY, HY, x0, x1, wg.hy(m))
		         +verlineint(l, m,   HY, HY, x0, x1, wg.hy(m)));
		pt *= (wg.eps(l,m)-wg.eps(l,m+1));
		nrm = recint(HY, HY, XYplane);
		pt *= val_k0(wg.lambda)*val_k0(wg.lambda)/2.0/beta/nrm;
		return pt;
	}

	pt  = 0.5*(verlineint(l, m+1, EX, EX, x0, x1, wg.hy(m)) 
	          +verlineint(l, m,   EX, EX, x0, x1, wg.hy(m))); 
	pt += 0.5*verlineint(l, m+1, EY, EY, x0, x1, wg.hy(m))
	         *wg.eps(l,m+1)/wg.eps(l,m); 
	pt += 0.5*verlineint(l, m,   EY, EY, x0, x1, wg.hy(m))
	         *wg.eps(l,m)/wg.eps(l,m+1); 
	pt += 0.5*(verlineint(l, m+1, EZ, EZ, x0, x1, wg.hy(m)) 
	          +verlineint(l, m,   EZ, EZ, x0, x1, wg.hy(m))); 
	pt *= (wg.eps(l,m)-wg.eps(l,m+1));

	nrm  = recint(HY, EX, XYplane);
	nrm -= recint(HX, EY, XYplane);
	nrm *= 2.0;

	pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm;

	return pt;
}

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

/* overlap 0.5*intint mode_EX*inif_HY - mode_EY*inif_HX dxdy with 
   an initial field, with real HX and HY components, integrals are 
   evaluated by gaussian quadrature formulas, separately on the
   fundamental rectangles, inside the domain rectangle r */

/* helper functions */
double EIMode::ovleval(double (*inif)(Fcomp cp, double x, double y), 
		     double x, double y) const
{
	double fex, fhy, fey, fhx;
	fhy = inif(HY, x, y);
	fex = 0; if(fhy != 0.0) fex = field(EX, x, y);
	fhx = inif(HX, x, y);
	fey = 0; if(fhx != 0.0) fey = field(EY, x, y);
 	return 0.5*(fex*fhy-fey*fhx); 
}

/* based on: Numerical Recipes in C --- The Art of Scientific Computing
   Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */ 
double EIMode::ovlpix(double (*inif)(Fcomp cp, double x, double y), 
                    double x0, double x1, double y, int numx) const
{
	int j;
	double xr,xm,dx,s;
	static double x[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	double ovl = 0.0; 
	double a, b, step;
	int i;
	step = (x1-x0)/((double) numx);
	for(i=0; i<=numx-1; ++i)
	{
		a = x0+((double) i)*step;
		b = a+step;
		xm=0.5*(b+a);
		xr=0.5*(b-a);
		s=0.0;
		for (j=1;j<=5;j++) {
			dx=xr*x[j];
			s += w[j]*(ovleval(inif, xm+dx, y)
				  +ovleval(inif, xm-dx, y));
		}
		ovl += s*xr;
	}
	return ovl;
}

double EIMode::ovlpint(double (*inif)(Fcomp cp, double x, double y), Rect r,
		       int numx, int numy) const
{
	double ovl = 0.0;
	double a, b, step;
	int i;
	int j;
	double yr,ym,dy,s;
	static double y[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	step = (r.y1-r.y0)/((double) numy);
	for(i=0; i<=numy-1; ++i)
	{
		a = r.y0+((double) i)*step;
		b = a+step;
		ym=0.5*(b+a);
		yr=0.5*(b-a);
		s=0.0;
		for (j=1;j<=5;j++) {
			dy=yr*y[j];
			s += w[j]*(ovlpix(inif, r.x0, r.x1, ym+dy, numx)
				  +ovlpix(inif, r.x0, r.x1, ym-dy, numx));
		}
		ovl += s*yr;
	}
	return ovl;
}

#define HDIST 1.0e-8
double EIMode::ovlp(double (*inif)(Fcomp cp, double x, double y), 
		  Rect r, int numx, int numy) const
{
	double ovl = 0.0; 
	double xb, yl;
	double xt, yr;
	int l, m, dum;
	Rect rp;
	if(r.x0 > r.x1) { xb=r.x1; r.x1=r.x0; r.x0=xb; }
	if(r.y0 > r.y1) { yl=r.y1; r.y1=r.y0; r.y0=yl; }
	yl = r.y0;
	while(fabs(r.y1-yl) > HDIST)
	{
		wg.rectidx(0.0, yl+HDIST/2.0, dum, m);
		yr = wg.rectbounds(0,m).y1;
		if(r.y1<yr) yr = r.y1;
		xb = r.x0;	
		while(fabs(r.x1-xb) > HDIST)
		{
			wg.rectidx(xb+HDIST/2.0, 0.0, l, dum);
			xt = wg.rectbounds(l,0).x1;
			if(r.x1<xt) xt = r.x1;
			rp.x0 = xb; rp.y0 = yl; rp.x1 = xt; rp.y1 = yr;
			ovl += ovlpint(inif, rp, numx, numy); 
			xb = xt;
		}
		yl = yr;
	}
	return ovl;
}

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

/* write single component to MATLAB .m file  
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR  
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters
   pltype: 'C': contour plot 
           'S': surface plot
           'I': intensity image 
           'N': field + mesh only, no plot commands (default) */
void EIMode::plot(Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		  char ext0, char ext1, char pltype) const
{
	FILE *dat;
	char name[13] = "t_______.m";
	double minf, maxf;
	name[1] = polchr(pol);
	name[2] = afochr(foa);
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	name[7] = pltype;
	fprintf(stderr, "%c%c >> %s\n", fldchr(cp), cpchr(cp), name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "VEIMS mode profile");
	maxf = 1.0;
	minf = -1.0;
	switch(fldchr(cp))
	{
	case 'E': 
		switch(foa)
		{
		case ORG: case REP: case IMP:
			maxf = fabs(maxE);
			if(fabs(minE)>maxf) maxf=fabs(minE);
			minf = -maxf;
			break;
		case MOD:  
			maxf = fabs(maxE);
			if(fabs(minE)>maxf) maxf=fabs(minE);
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxE*maxE;
			if(minE*minE>maxf) maxf=minE*minE;
			minf = 0.0;
			break;
		}
		break;
	case 'H': 
		switch(foa)
		{
		case ORG: case REP: case IMP:
			maxf = fabs(maxH);
			if(fabs(minH)>maxf) maxf=fabs(minH);
			minf = -maxf;
			break;
		case MOD:  
			maxf = fabs(maxH);
			if(fabs(minH)>maxf) maxf=fabs(minH);
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxH*maxH;
			if(minH*minH>maxf) maxf=minH*minH;
			minf = 0.0;
			break;
		}
		break;
	case 'S': 
		switch(foa)
		{
		case ORG: case REP: case MOD:
			maxf = maxS;
			minf = 0.0;
			break;
		case SQR:  // not reasonable
			maxf = maxS*maxS;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;
			break;
		}
		break;
	}
	if(pltype == 'I') mlout_gengeoxy(dat, st, disp); 
	else              mlout_geo(dat, wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;		
	fld = field(disp, npx, npy, cp, foa);
	mlout_fld(dat, npx, npy, cp, fld);
	name[8] = 0; 
	switch(pltype)
	{
		case 'C': mlout_contour(name, dat, cp, foa); break;
		case 'S': mlout_surface(name, dat, cp, foa); break;
		case 'I': 
			mlout_image(name, dat, cp, foa, minf, maxf); 
			if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
			else                         mlout_magcolormap(dat);
			break;
		default: break;
	}
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}


/* write all components to MATLAB .m file, surface plots  
   ext0, ext1: filename id characters
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh */
void EIMode::acplot(Rect disp, int npx, int npy, char ext0, char ext1) const
{
	FILE *dat;
	Dmatrix fld;		
	char name[13] = "t___A.m";
	name[1] = polchr(pol);
	name[2] = ext0;
	name[3] = ext1;
	fprintf(stderr, "Ex - Hz >> %s\n", name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "VEIMS mode profile");
	name[5] = 0;  
	mlout_geo(dat, wg, 0.0, 1.0);
	mlout_meshxy(dat, disp, npx, npy);
	fld = field(disp, npx, npy, EX, ORG);
	mlout_fld(dat, npx, npy, EX, fld);
	fld = field(disp, npx, npy, EY, ORG);
	mlout_fld(dat, npx, npy, EY, fld);
	fld = field(disp, npx, npy, EZ, ORG);
	mlout_fld(dat, npx, npy, EZ, fld);
	fld = field(disp, npx, npy, HX, ORG);
	mlout_fld(dat, npx, npy, HX, fld);
	fld = field(disp, npx, npy, HY, ORG);
	mlout_fld(dat, npx, npy, HY, fld);
	fld = field(disp, npx, npy, HZ, ORG);
	mlout_fld(dat, npx, npy, HZ, fld);

	mlout_acmfile(name, dat, minE, maxE, minH, maxH);
	mlout_print(dat, name, 'e');
	
	fclose(dat);
	return;
}

/* write profile section along line (x0,y0) -> (x1,y1) to MATLAB .m file
	cp: EX - HZ, SZ
	ext0, ext1: filename id characters
	np: number of output points
	pltype: 'L': geometry information + plot commands (default)
	        'V': field, mesh, and plot command, to be included into *L.m 
   --- at present implemented for horizontal or vertical lines only ---
*/
#define HDISTF 1.0e-8
void EIMode::secplot(Fcomp cp, char ext0, char ext1, int np, char pltype,
  	             double x0, double y0, double x1, double y1) const
{
	FILE *dat;
	int i, j;
	double x, y, dx, dy, t;
	double xbd, ybd;
	double epsold, epsnew;
	double minf, maxf, f;
	int nsec, nbd;
	char name[13] = "t______.m";
	char ori;
	name[1] = polchr(pol);
	name[2] = fldchr(cp);
	name[3] = cpchr(cp);
	name[4] = ext0;
	name[5] = ext1;
	if(pltype != 'V') pltype = 'L';
	name[6] = pltype;

	dx = fabs(x1-x0);
	dy = fabs(y1-y0);
	
	if(dx > 1.0e-10 && dy > 1.0e-10)
	{
		fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> file \n", 
	                        fldchr(cp), cpchr(cp), x0, y0, x1, y1);
		fprintf(stderr, 
                " not implemented for tilted lines, use EIMode.writesec().\n"); 
		return;
	}
	if(dx < 1.0e-10 && dy < 1.0e-10)
	{
		fprintf(stderr, "EIMode.secplot(): Nothing to do !\n");
		return;
	}
	if(dx < 1.0e-10) ori = 'h';
	else             ori = 'v';
	
	if(ori == 'v') { if(x0 > x1) { t = x0; x0 = x1; x1 = t; }}
	else           { if(y0 > y1) { t = y0; y0 = y1; y1 = t; }}

	fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> %s\n", 
	                fldchr(cp), cpchr(cp), x0, y0, x1, y1, name);

	Dmatrix fld(wg.nx+2+wg.ny+2, np+1);
	Dmatrix pos(wg.nx+2+wg.ny+2, np+1);
	Ivector nump(wg.nx+2+wg.ny+2);
	Dvector of(np+1);
	Dvector op(np+1);
	nsec = 0;
	Dvector bd(wg.nx+1+wg.ny+1);
	Dvector ri(wg.nx+2+wg.ny+2);
	nbd  = 0;

	dx = x1-x0; dy = y1-y0;
	nump(nsec) = 0;
	x = x0;
	y = y0;
	epsold = wg.eps(x, y);  	
	ri(nbd) = sqrt(epsold);
	if(ori=='h') pos(nsec, nump(nsec)) = y;  
	else         pos(nsec, nump(nsec)) = x;  
	fld(nsec, nump(nsec)) = field(cp, x, y);  
	++nump(nsec);
	for(i=1; i<=np; ++i)
	{
		x = x0+i*dx/np;
		y = y0+i*dy/np;
		epsnew = wg.eps(x, y);
		if(fabs(epsnew-epsold)<1.0e-10)
		{
			if(ori=='h') pos(nsec, nump(nsec)) = y;  
			else         pos(nsec, nump(nsec)) = x;  
			fld(nsec, nump(nsec)) = field(cp, x, y);  
			++nump(nsec);
		}
		else
		{
			if(ori=='h') { xbd = x0; 
                                       ybd = (wg.rectbounds(x, y)).y0; }
			else         { xbd = (wg.rectbounds(x, y)).x0; 
                                       ybd = y0; }
			if(ori=='h') pos(nsec, nump(nsec)) = ybd;  
			else         pos(nsec, nump(nsec)) = xbd;  
			fld(nsec, nump(nsec)) = 	
				field(cp, xbd-dx*HDISTF, ybd-dy*HDISTF);  
			++nump(nsec);
			if(ori=='h') bd(nbd) = ybd; 
			else         bd(nbd) = xbd; 
			++nbd;
			epsold = epsnew;
			ri(nbd) = sqrt(epsold);
			++nsec;
			nump(nsec) = 0;
			if(ori=='h') pos(nsec, nump(nsec)) = ybd;  
			else         pos(nsec, nump(nsec)) = xbd;  
			fld(nsec, nump(nsec)) = 	
				field(cp, xbd+dx*HDISTF, ybd+dy*HDISTF);  
			++nump(nsec);
			if(ori=='h') pos(nsec, nump(nsec)) = y;  
			else         pos(nsec, nump(nsec)) = x;  
			fld(nsec, nump(nsec)) = field(cp, x, y);  
			++nump(nsec);
		}
	}
	++nsec;

	dat = fopen(name, "w+");
	mlout_title(dat, name, "VEIMS mode field section");
	name[7] = 0;
	minf = -1.0;
	maxf = 1.0;
	if(pltype == 'L')
	{
		switch(fldchr(cp))
		{
			case 'E': minf = minE; maxf = maxE; break;
			case 'H': minf = minH; maxf = maxH; break;
			case 'S': minf = minS; maxf = maxS; break;
		}
		mlout_geo(dat, wg, minf, maxf); 
		mlout_Lsecgeo(dat, x0, y0, x1, y1, nbd, bd, ri);
	}
	minf = fld(0, 0);
	maxf = fld(0, 0);
	
	for(j=0; j<=nsec-1; ++j)
	{
		for(i=0; i<=nump(j)-1; ++i)
		{
			f = fld(j, i);
			if(f < minf) minf = f;
			if(f > maxf) maxf = f;
			of(i) = f;
			op(i) = pos(j, i);
		}
		if(pltype == 'L')
			mlout_sec1D(dat, cp, dig10(j), dig1(j), ' ', ' ', 
                                    nump(j), of, op); 
		else
			mlout_sec1D(dat, cp, dig10(j), dig1(j), ext0, ext1, 
                                    nump(j), of, op); 
	}
	if(pltype == 'L') 
	{
		mlout_Lsecplot(name, dat, ori, minf, maxf, cp, nbd, nsec); 
		mlout_print(dat, name, 'e');
	}
	else mlout_Vsecplot(dat, cp, nsec, ext0, ext1); 
	fclose(dat);
	return;
}


/* write single component to MATLAB .m file, fancy style :-)  
   cp: EX - HZ, SZ
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters */
void EIMode::fplot(Fcomp cp, Rect disp, int npx, int npy, 
                   char ext0, char ext1) const
{
	FILE *dat;
	int np, l, m;
	double x0, x1, xp, y0, y1, yp;
	int numc;
	char name[13] = "t_____F.m";
	name[1] = polchr(pol);
	name[2] = fldchr(cp);
	name[3] = cpchr(cp);
	name[4] = ext0;
	name[5] = ext1;
	fprintf(stderr, "%c%c >> %s\n", fldchr(cp), cpchr(cp), name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "VEIMS mode profile :-)");
	name[7] = 0; 

	switch(fldchr(cp))
	{
		case 'E': mlout_geo(dat, wg, minE, maxE); break;
		case 'H': mlout_geo(dat, wg, minH, maxH); break;
		case 'S': mlout_geo(dat, wg, minS, maxS); break;
	}
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = field(disp, npx, npy, cp, ORG);
	mlout_fld(dat, npx, npy, cp, fld);

	Dvector fv;
	numc = 0;
	for(l=0; l<=wg.nx+1; ++l)
	{
		if(l==0) x0 = disp.x0; 	
                else     x0 = wg.hx(l-1);
		if(l==wg.nx+1) x1 = disp.x1; 	
                else           x1 = wg.hx(l);
		for(m=0; m<=wg.ny; ++m)
		{
			if(fabs(wg.n(l,m)-wg.n(l,m+1)) > 1.0e-10)
			{
				yp = wg.hy(m);
				np = (int)(((double) npx)
                                     *(x1-x0)/(disp.x1-disp.x0)); 
				if(np >= 2)
				{
				fv = field(cp, np, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				fv = field(cp, np, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	for(m=0; m<=wg.ny+1; ++m)
	{
		if(m==0) y0 = disp.y0; 	
                else     y0 = wg.hy(m-1);
		if(m==wg.ny+1) y1 = disp.y1; 	
                else           y1 = wg.hy(m);
		for(l=0; l<=wg.nx; ++l)
		{
			if(fabs(wg.n(l,m)-wg.n(l+1,m)) > 1.0e-10)
			{
				xp = wg.hx(l);
				np = (int)(((double)npy)
                                     *(y1-y0)/(disp.y1-disp.y0)); 
				if(np >= 2)
				{
				fv = field(cp, np, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST); 
				mlout_sec(dat, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				fv = field(cp, np, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST);
				mlout_sec(dat, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	fv = field(cp, npx, disp.x0, disp.y0, disp.x1, disp.y0);
	mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, 
                  npx, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(cp, npy, disp.x1, disp.y0, disp.x1, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, 
                  npy, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(cp, npx, disp.x1, disp.y1, disp.x0, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, 
                  npx, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(cp, npy, disp.x0, disp.y1, disp.x0, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, 
                  npy, cp, dig10(numc), dig1(numc), fv);
	++numc;

	mlout_fancy(name, dat, cp, numc);  	
	mlout_print(dat, name, 'p');

	fclose(dat);
	return;
}

/* export full mode profile data ("all") into a viewer m-file 
   disp: the plot window   
   npx, npy: number of plot points in the x and y directions
   ext0, ext1: filename id characters */
void EIMode::viewer(Rect disp, int npx, int npy, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "prf__A.m";
	double wl = wg.lambda;
	name[3] = ext0;
	name[4] = ext1;
	double xbeg, xend, ybeg, yend;
	xbeg = disp.x0; xend = disp.x1;
	ybeg = disp.y0; yend = disp.y1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   xbeg, npx, xend, ybeg, npy, yend, name);
	dat = fopen(name, "w+");
	name[6] = 0; 
	mlout_viewertopxy(dat, name, pol, wl);
	mlout_gengeoxy(dat, st, disp);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix f;
	Fcomp cp;
	if(pol == TE)
	{
		cp = EX;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EZ;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtoim(dat, cp);
		cp = HX;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HY;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HZ;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtoim(dat, cp);
	}
	else
	{
		cp = EX;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EZ;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtoim(dat, cp);
		cp = HX;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HY;
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtore(dat, cp);
		mlout_0fldxy(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HZ;
		mlout_0fldxy(dat, cp);
		mlout_fldtore(dat, cp);
		f = field(disp, npx, npy, cp, ORG); 
		mlout_fld(dat, npx, npy, cp, f);
		mlout_fldtoim(dat, cp);
	}
	Dmatrix n(npx, npy);
	double dx = (xend-xbeg)/((double)(npx-1));
	double dy = (yend-ybeg)/((double)(npy-1));
	double x, y;
	y = ybeg;
	for(int j=0; j<=npy-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) { n(i, j) = st.n(x, y); x += dx; }
		y += dy;
	}
	mlout_n(dat, npx, npy, n);
	mlout_fldviewerxy(dat, name);
	fclose(dat);
	return;
}

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

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

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

/* copy constructor */
EIModeArray::EIModeArray(const EIModeArray& ma)
	    : num(ma.num)
{
	mvec = new EIMode[num];
	EIMode* ap = ma.mvec;
	EIMode* mp =    mvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}

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

}

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

}

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

/* add a mode */
void EIModeArray::add(EIMode m)
{
	EIMode *avec;
	avec = new EIMode[num+1];
	EIMode* ap = avec;
	EIMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	*ap = m;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	++num;
	return;
}

/* delete a mode entry */
void EIModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		eimodeerror("EIModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	EIMode *avec;
	avec = new EIMode[num-1];
	EIMode* ap = avec;
	EIMode* mp = mvec;
	for(int j=0; j<=i-1; ++j) *ap++ = *mp++;
	mp++;
	for(int j=i+1; j<=num-1; ++j) *ap++ = *mp++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	--num;
	return;
}

/* add an entire EIModeArray nma */ 
void EIModeArray::merge(EIModeArray& ma)
{
	EIMode *avec;
	avec = new EIMode[num+ma.num];
	EIMode* ap = avec;
	EIMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	EIMode* np = ma.mvec;
	for(int i=0; i<=ma.num-1; ++i) *ap++ = *np++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	num += ma.num;
	return;
}

/* sort the array by propagation constants, highest first */
void EIModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	EIMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		maxi = j; 
		maxb = (mvec[j]).beta;
		for(k=j+1; k<=num-1; ++k)
		{
			if(maxb<(mvec[k]).beta)
			{
				maxb = (mvec[k]).beta;
				maxi = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t;
	}
	return;
}

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

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

/* field superposition at point (x, y, z), 
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp:   EX, EY, EZ, HX, HY, HZ, SZ */
Complex EIModeArray::field(Cvector amp, Cvector pert, Fcomp cp, 
                           double x, double y, double z) const
{
	Complex s, iphase;
	double f, sum;
	Cvector a;
	int i, j, k;
	a = amp;	
	if(amp.nel < num)
	{
		a = Cvector(num);
		a.init(CC0);
		for(i=0; i<=amp.nel-1; ++i) a(i) = amp(i);
	}
	Cvector cbet(num);
	for(i=0; i<=num-1; ++i) cbet(i) = pert(i)+(mvec[i]).beta;
	if(z != 0.0)
	{
		for(i=0; i<=num-1; ++i) 
		{
			iphase = CCI*cbet(i)*(-z);
			a(i) = a(i)*exp(iphase);
		}
	}
	s = CC0;
	if(cp != SZ)
	{
		for(i=0; i<=num-1; ++i)
		{
			f = (mvec[i]).field(cp, x, y);  
			if(cp != EZ && cp != HZ) 
				s = s+a(i)*f;
			else
			{
				a(i) = CCI*a(i);
				s = s+a(i)*f;
			}
			
		}
		return s;
	}
	sum = 0.0;
	Dvector fEx(num); 
	Dvector fEy(num); 
	Dvector fHx(num); 
	Dvector fHy(num); 
	for(j=0; j<=num-1; ++j) fEx(j) = (mvec[j]).field(EX, x, y);
	for(j=0; j<=num-1; ++j) fEy(j) = (mvec[j]).field(EY, x, y);
	for(j=0; j<=num-1; ++j) fHx(j) = (mvec[j]).field(HX, x, y);
	for(j=0; j<=num-1; ++j) fHy(j) = (mvec[j]).field(HY, x, y);
	for(j=0; j<=num-1; ++j)
	{
		for(k=0; k<=num-1; ++k)
		{
			f = (a(j).re*a(k).re + a(j).im*a(k).im);
			sum += f*(fEx(j)*fHy(k)-fEy(j)*fHx(k));
		}
	}
	s = Complex(0.5*sum);
	return s;
}

/* evaluate component cp on a rectangular npx x npy mesh on area disp 
   at position z
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ 
   foa: MOD, SQR, REP, IMP */
Dmatrix EIModeArray::field(Cvector amp, Cvector pert, double z,
                           Rect disp, int npx, int npy, 
                           Fcomp cp, Afo foa) const
{
	double x, y;
	double dx, dy;
	int i, j;
	if(npx <= 1) eimodeerror("field: npx <= 1");
	if(npy <= 1) eimodeerror("field: npy <= 1");
	Dmatrix f(npx, npy);
	Complex ft;
	double ftd;
	dx = (disp.x1-disp.x0)/(npx-1);
	dy = (disp.y1-disp.y0)/(npy-1);
	for(i=0; i<=npx-1; ++i)
	{
		x = disp.x0+i*dx;	
		for(j=0; j<=npy-1; ++j)
		{
			y = disp.y0+j*dy;	
			ft = field(amp, pert, cp, x, y, z);
			ftd = 0.0; 
			switch(foa)
			{
			case MOD: case ORG: 
				ftd = ft.abs(); 
				break;
			case SQR: 
				ftd = ft.sqabs(); 
				break;
			case REP:  
				ftd = ft.re;
				break;
			case IMP:  
				ftd = ft.im;
				break;
			}
			f(i, j) = ftd;
		}
	}
	return f;
} 

/* store nump values of component cp between (x0, y0) and (x1, y1) in a vector 
   at position z
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ 
   foa: MOD, SQR, REP, IMP */
Dvector EIModeArray::field(Cvector amp, Cvector pert, double z, 
				Fcomp cp, Afo foa, int nump, 
                                double x0, double y0, 
				double x1, double y1) const
{
	double x, y;
	double dx, dy;
	int j;
	if(nump <= 1) eimodeerror("field: nump <= 1");
	Dvector f(nump);
	Complex ft; 
	double ftd;
	dx = (x1-x0)/(nump-1);
	dy = (y1-y0)/(nump-1);
	for(j=0; j<=nump-1; ++j)
	{
		x = x0+j*dx;	
		y = y0+j*dy;	
		ft = field(amp, pert, cp, x, y, z);
		ftd = 0.0;
		switch(foa)
		{
		case MOD: case ORG: 
			ftd = ft.abs(); 
			break;
		case SQR: 
			ftd = ft.sqabs(); 
			break;
		case REP:  
			ftd = ft.re;
			break;
		case IMP:  
			ftd = ft.im;
			break;
		}
		f(j) = ftd;
	}
	return f;
}

/* - Three segment coupler, power transfer evaluation -------------------- 
   imode: input mode, excites the modes this(j) with rel. power 1 at z=0 
   omode: output mode, relative power is returned 
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) 
(!) this(j) are assumed to be normalized. */

/* weights for the power evaluation:
   w = ( omode | this(m) ) ( this(m) | imode ) */
double EIModeArray::pweight(const EIMode& imode, int m,
                              const EIMode& omode) 
{
	return lscalprod(omode, (mvec[m]))*lscalprod((mvec[m]), imode);
}
   
/* single relative power level for device length l */    
double EIModeArray::iopower(const EIMode& imode, const EIMode& omode, 
                              Cvector pert, double l)
{
	Complex a, am, cbet;
	int m;
	a=CC0;
	for(m=0; m<=num-1; ++m)
	{
		cbet = pert(m)+(mvec[m]).beta;
		cbet = CCI*cbet*(-l);
		am = exp(cbet)*pweight(imode, m, omode);
		a = a+am;
	}
	return a.sqabs();
}

/* output power for numl devices of lengths between lmin and lmax */    
Dvector EIModeArray::iopower(const EIMode& imode, const EIMode& omode, 
                               Cvector pert, 
                               int numl, double lbeg, double lend)
{
	Complex a, am, cbet;
	a=CC0;
	Dvector w(num);
	double l, dl;
	int li, m;
	if(numl <= 0) eimodeerror("iopower: numl <= 0");
	Dvector p(numl);
	if(numl == 1) 
	{
		p(0) = iopower(imode, omode, pert, 0.5*(lbeg+lend));
		return p;
	}
	for(m=0; m<=num-1; ++m)
	{
		w(m) = pweight(imode, m, omode);	
	}
	dl = (lend-lbeg)/(numl-1.0);
	for(li=0; li<=numl-1; ++li)
	{
		a=CC0;
		l = lbeg+li*dl;
		for(m=0; m<=num-1; ++m)
		{
			cbet = pert(m)+(mvec[m]).beta;
			cbet = CCI*cbet*(-l);
			am = exp(cbet)*w(m);
			a = a+am;
		}
		p(li) = a.sqabs();
	}
	return p;
}

/* ... write to file */    
void EIModeArray::writeiopower(const EIMode& imode, const EIMode& omode, 
                                Cvector pert, 
                                int numl, double lbeg, double lend,
                                char ext0, char ext1)
{
	Dvector p;
	int i;
	double l, dl;
	char name[13] = "__pow__.xyf";
	FILE *dat;
	if(numl <= 0) eimodeerror("iopower: numl <= 0");
	switch((mvec[0]).pol)
	{
		case TE: name[0] = 't'; name[1] = 'e'; break;
		case TM: name[0] = 't'; name[1] = 'm'; break;
	}
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "P(L = %g -> %g) >> %s\n", lbeg, lend, name);
	p = iopower(imode, omode, pert, numl, lbeg, lend); 
	dat = fopen(name, "w+");
	if(p.nel == 1)
	{
		fprintf(dat, "%g %.10g\n", 0.5*(lbeg+lend), p(0));
		fclose(dat);
		return;
	}
	dl = (lend-lbeg)/(p.nel-1.0);
	for(i=0; i<=p.nel-1; ++i)
	{
		l = lbeg+i*dl;
		fprintf(dat, "%g %.10g\n", l, p(i));
	}
	fclose(dat);
	return;
}

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

/* write single component of the interference field at position z 
   to MATLAB .m file  
   amp: complex amplitudes at z=0
	amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR, REP, IMP  
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters
   pltype: 'C': contour plot
           'S': surface plot
           'I': intensity image 
           'N': field + mesh only, no plot commands (default) */
void EIModeArray::plot(Cvector amp, Cvector pert, double z, 
                       Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		       char ext0, char ext1, char pltype) const
{
	FILE *dat;
	char name[13] = "int______.m";
	double rz;
	double max, f, maxsum, minf, maxf; 
	int i;
	rz = floor(z*10.0+0.5)/10.0;
	name[3] = afochr(foa);
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	name[6] = ext0;
	name[7] = ext1;
	name[8] = pltype;
	fprintf(stderr, "%c%c(z=%g) >> %s\n", 
                        fldchr(cp), cpchr(cp), z, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxf = 1.0;
	minf = -1.0;
	switch(fldchr(cp))
	{
	case 'E':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((mvec[i]).maxE)*amp(i).abs();
			f = fabs((mvec[i]).minE)*amp(i).abs();
			if(max < f) max = f;
			maxsum += max;
		}
		switch(foa)
		{
		case MOD: 
			maxf = maxsum;
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case ORG: case REP: case IMP:
			maxf = maxsum;
			minf = -maxsum;
			break;
		}
		break;
	case 'H':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((mvec[i]).maxH)*amp(i).abs();
			f = fabs((mvec[i]).minH)*amp(i).abs();
			if(max < f) max = f;
			maxsum += max;
		}
		switch(foa)
		{
		case MOD: 
			maxf = maxsum;
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case ORG: case REP: case IMP:
			maxf = maxsum;
			minf = -maxsum;
			break;
		}
		break;
	case 'S':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
			f = fabs((mvec[i]).minS)*amp(i).abs()*amp(i).abs();
			if(max < f) max = f;
			maxsum += sqrt(max);
		}
		switch(foa)
		{
		case MOD: case ORG: case REP: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case SQR:  // not reasonable 
			maxf = maxsum*maxsum*maxsum*maxsum;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;	
			break;
		}
		break;
	}
	if(pltype == 'I') mlout_gengeoxy(dat, (mvec[0]).st, disp); 
	else              mlout_geo(dat, (mvec[0]).wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = field(amp, pert, z, disp, npx, npy, cp, foa);
	mlout_fld(dat, npx, npy, cp, fld);
	name[9] = 0; 
	switch(pltype)
	{
		case 'C': 
			mlout_contour(name, dat, cp, foa); 
			mlout_annotatezpos(dat, rz);
			break;
		case 'S': 
			mlout_surface(name, dat, cp, foa); 
			break;
		case 'I': 
			mlout_image(name, dat, cp, foa, minf, maxf); 
			if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
			else                         mlout_magcolormap(dat);
			mlout_annotatezpos(dat, rz);
			break;
		default: break;
	}
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}

/* write interference field at position z to MATLAB .m file, fancy style :-)  
   amp: complex amplitudes at z=0
	amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters */
#define HDIST 1.0e-8
void EIModeArray::fplot(Cvector amp, Cvector pert, double z,
                        Rect disp, int npx, int npy, 
                        char ext0, char ext1) const
{
	FILE *dat;
	int np, l, m;
	double x0, x1, xp, y0, y1, yp;
	int numc;
	double max, f, maxsum, minf, maxf; 
	int i;
	WgCrs wg;
	char name[13] = "intfSz__F.m";
	name[6] = ext0;
	name[7] = ext1;
	fprintf(stderr, "Sz(z=%g) >> %s\n", z, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field  :-)");
	name[9] = 0; 
	wg = (mvec[0]).wg;
	maxsum = 0.0;
	for(i=0; i<=num-1; ++i)
	{
		max = fabs((mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
		f = fabs((mvec[i]).minS)*amp(i).abs()*amp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum*maxsum;
	minf = 0.0;
	mlout_geo(dat, wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = field(amp, pert, z, disp, npx, npy, SZ, ORG);
	mlout_fld(dat, npx, npy, SZ, fld);

	Dvector fv;
	numc = 0;
	for(l=0; l<=wg.nx+1; ++l)
	{
		if(l==0) x0 = disp.x0; 	
                else     x0 = wg.hx(l-1);
		if(l==wg.nx+1) x1 = disp.x1; 	
                else           x1 = wg.hx(l);
		for(m=0; m<=wg.ny; ++m)
		{
			if(fabs(wg.n(l,m)-wg.n(l,m+1)) > 1.0e-10)
			{
				yp = wg.hy(m);
				np = (int)(((double) npx)
                                     *(x1-x0)/(disp.x1-disp.x0)); 
				if(np >= 2)
				{
				fv = field(amp, pert, z, SZ, ORG, np, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				fv = field(amp, pert, z, SZ, ORG, np, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	for(m=0; m<=wg.ny+1; ++m)
	{
		if(m==0) y0 = disp.y0; 	
                else     y0 = wg.hy(m-1);
		if(m==wg.ny+1) y1 = disp.y1; 	
                else           y1 = wg.hy(m);
		for(l=0; l<=wg.nx; ++l)
		{
			if(fabs(wg.n(l,m)-wg.n(l+1,m)) > 1.0e-10)
			{
				xp = wg.hx(l);
				np = (int)(((double)npy)
                                     *(y1-y0)/(disp.y1-disp.y0)); 
				if(np >= 2)
				{
				fv = field(amp, pert, z, SZ, ORG, np, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST); 
				mlout_sec(dat, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				fv = field(amp, pert, z, SZ, ORG, np, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST);
				mlout_sec(dat, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	fv = field(amp, pert, z, SZ, ORG, 
                       npx, disp.x0, disp.y0, disp.x1, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, 
                  npx, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(amp, pert, z, SZ, ORG, 
                       npy, disp.x1, disp.y0, disp.x1, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, 
                  npy, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(amp, pert, z, SZ, ORG, 
                       npx, disp.x1, disp.y1, disp.x0, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, 
                  npx, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = field(amp, pert, z, SZ, ORG, 
                       npy, disp.x0, disp.y1, disp.x0, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, 
                  npy, SZ, dig10(numc), dig1(numc), fv);
	++numc;

	mlout_fancy(name, dat, SZ, numc);  	
	mlout_print(dat, name, 'p');

	fclose(dat);
	return;
}


/* animate the interference field  
   amp: complex amplitudes at z=0
        amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR, REP, IMP 
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   pltype: 'C': contour plot 
           'S': surface plot
           'I': intensity image 
           'F': fancy plot, cp and foa are set to SZ, MOD (default)
   z0, z1, numz: .m files are generated for z-positions 
                  z0+i*(z1-z0)/numz, i=0, ..., numz-1 */
void EIModeArray::movie(Cvector amp, Cvector pert,  
                        Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		        char pltype,
                        int numz, double z0, double z1) const
{
	int i;
	double z;
	FILE *dat;
	char name[13] = "int___plM.m";
	if(numz >= 100) numz = 100;
	if(pltype == 'F') { cp = SZ; foa = ORG; }
	name[3] = afochr(foa);
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	for(i=0; i<=numz-1; ++i)
	{
		z = z0+i*(z1-z0)/numz; 
		switch(pltype)
		{
		case 'C': case 'S': case 'I':
			plot(amp, pert, z, cp, foa, 
                             disp, npx, npy, dig10(i), dig1(i), pltype);
			break;
		default:
			fplot(amp, pert, z, 
                              disp, npx, npy, dig10(i), dig1(i));
			break;
		}
	}
	fprintf(stderr, ">> %s\n", name);	
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field animation");
	name[8] = pltype;
	name[9] = 0;
	mlout_play(dat, name, numz);
	fclose(dat);
	return;
}


/* interference pattern on the horizontal plane 
   ybeg <= y <= yend, zbeg <= z <= zend at position x
   image plot corresponding to the squareroot of the local intensity (SZ)
   amp: complex amplitudes at z=0
        amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   npy, npz: number of points in output mesh
   ext0, ext1: filename id characters */
void EIModeArray::prop(Cvector amp, Cvector pert, double x,
                         double ybeg, double yend, int npy,
                         double zbeg, double zend, int npz,
                         char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "propSz__I.m";
	double max, f, maxsum, minf, maxf; 
	double y, dy, z, dz;
	Complex ph;
	int i, j, k, m;
	name[6] = ext0;
	name[7] = ext1;
	fprintf(stderr, "prop(z=%g -> %g) >> %s\n", zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxsum = 0.0;
	for(i=0; i<=num-1; ++i)
	{
		max = fabs((mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
		f = fabs((mvec[i]).minS)*amp(i).abs()*amp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum;
	minf = 0.0;
	mlout_geoyz(dat, (mvec[0]).wg, minf, maxf);
	mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz);

	Dmatrix fld(npy, npz);
	dy = (yend-ybeg)/(npy-1);
	dz = (zend-zbeg)/(npz-1);

	Dmatrix pmprofEx(num, npy);
	Dmatrix pmprofEy(num, npy);
	Dmatrix pmprofHx(num, npy);
	Dmatrix pmprofHy(num, npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		for(m=0; m<=num-1; ++m)
		{
			pmprofEx(m, i) = (mvec[m]).field(EX, x, y);
			pmprofEy(m, i) = (mvec[m]).field(EY, x, y);
			pmprofHx(m, i) = (mvec[m]).field(HX, x, y);
			pmprofHy(m, i) = (mvec[m]).field(HY, x, y);
		}
	}
	Cvector cbet(num);
	Cvector ma(num);
	for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(mvec[m]).beta;
	for(j=0; j<=npz-1; ++j)
	{
		z = zbeg+j*dz;
		for(m=0; m<=num-1; ++m) 
		{
			ph = CCI*cbet(m)*(-z);
			ma(m) = amp(m)*exp(ph);
		}
		for(i=0; i<=npy-1; ++i) 
		{
			y = ybeg+i*dy;
			f = 0.0;
			for(m=0; m<=num-1; ++m) 
			{
				for(k=0; k<=num-1; ++k)
				{
		f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im)
		    *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i));
				}
			}
			fld(i, j) = sqrt(0.5*fabs(f));
		}

	}
	mlout_fld(dat, npy, npz, SZ, fld);
	name[9] = 0; 
	mlout_propimage(name, dat); 
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* Three segment coupler, interference pattern on the horizontal plane 
   ybeg <= y <= yend, zbeg <= z <= zend at elevation x
   image plot corresponding to the squareroot of the local intensity (SZ)
       z < 0: intensity profile of the input modes
   0 < z < l: mode interference pattern
       l < z: intensity profiles of the output modes, 
              scaled by the output amplitudes
(!) this(j) are assumed to be normalized.
   iamp: complex input mode amplitudes
   imodes: input modes, excite the modes this(j) at z=0, 
           belonging to well separated port waveguides 
   iwg: the geometry of the combined input waveguides
   omodes: output modes, belonging to well separated port waveguides 
   owg: the geometry of the combined output waveguides
	iwg and owg are for displaying purposes only
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   npy, npz: number of points in output mesh
   ext0, ext1: filename id characters */
#define ZDIST 1.0e-8
void EIModeArray::ioprop(Cvector iamp, 
                           const EIModeArray& imodes, WgCrs iwg,
                           const EIModeArray& omodes, WgCrs owg,
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "propSz__I.m";
	double max, f, maxsum, minf, maxf; 
	int i, j, m, k;
	double y, z, dy, dz;
	Complex a, aj;
	Complex ph;
	double s;
	name[6] = ext0;
	name[7] = ext1;
	if(zbeg >= zend) eimodeerror("ioprop: invalid z-range");
	fprintf(stderr, "ioprop(z=%g -> %g) >> %s\n", zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxsum = 0.0;
	for(i=0; i<=imodes.num-1; ++i)
	{
		max = fabs(imodes(i).maxS)*iamp(i).abs()*iamp(i).abs();
		f = fabs(imodes(i).minS)*iamp(i).abs()*iamp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum;
	minf = 0.0;
	mlout_iopropgeo(dat, iwg, (mvec[0]).wg, owg, l, minf, maxf);
	mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz);
	s = 0.0;
	for(i=0; i<=imodes.num-1; ++i)
	{
		fprintf(stderr, "    A^I_%d = %g + i %g,  |A|^2: %g\n", 
                                i, iamp(i).re, iamp(i).im, iamp(i).sqabs());
		s += iamp(i).sqabs();

	}
	fprintf(stderr, "   total input power: %g\n", s); 
	Cvector pamp(num);
	for(i=0; i<=num-1; ++i)
	{
		a = CC0;	
		for(j=0; j<=imodes.num-1; ++j)
		{
			aj = iamp(j)*lscalprod((mvec[i]), imodes(j));
			a = a+aj;
		}
		pamp(i) = a;
	}
	s=0.0;
	for(i=0; i<=num-1; ++i)
	{
		fprintf(stderr, "    A^II_%d = %g + i %g,  |A|^2: %g\n", 
                                 i, pamp(i).re, pamp(i).im, pamp(i).sqabs());
		s += pamp(i).sqabs();
	}
	fprintf(stderr, "   propagating power: %g\n", s); 
	Cvector oamp(omodes.num);
	for(i=0; i<=omodes.num-1; ++i)
	{
		a = CC0;	
		for(j=0; j<=num-1; ++j)
		{
			aj = pamp(j)*lscalprod(omodes(i), (mvec[j]));
			ph = pert(j)+(mvec[j]).beta; 
			ph = CCI*ph*(-l);
			aj = aj*exp(ph);
			a = a+aj;
		}
		oamp(i) = a;
	} 
	s=0.0;
	for(i=0; i<=omodes.num-1; ++i)
	{
		fprintf(stderr, "    A^III_%d = %g + i %g,  |A|^2: %g\n", 
                                 i, oamp(i).re, oamp(i).im, oamp(i).sqabs());
		s += oamp(i).sqabs();
	}
	fprintf(stderr, "   total output power: %g\n", s); 

	Dmatrix fld(npy, npz);
	dy = (yend-ybeg)/(npy-1);
	dz = (zend-zbeg)/(npz-1);

	Dvector improf(npy);
	Dvector omprof(npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		improf(i) = 0.0;
		for(m=0; m<=imodes.num-1; ++m)
		{
			improf(i) += imodes(m).field(SZ, x, y)
		  	              *iamp(m).sqabs();
		}
		omprof(i) = 0.0;
		for(m=0; m<=omodes.num-1; ++m)
		{
			omprof(i) += omodes(m).field(SZ, x, y)
		  	              *oamp(m).sqabs();
		}
	}
	Dmatrix pmprofEx(num, npy);
	Dmatrix pmprofEy(num, npy);
	Dmatrix pmprofHx(num, npy);
	Dmatrix pmprofHy(num, npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		for(m=0; m<=num-1; ++m)
		{
			pmprofEx(m, i) = (mvec[m]).field(EX, x, y);
			pmprofEy(m, i) = (mvec[m]).field(EY, x, y);
			pmprofHx(m, i) = (mvec[m]).field(HX, x, y);
			pmprofHy(m, i) = (mvec[m]).field(HY, x, y);
		}
	}
	Cvector cbet(num);
	Cvector ma(num);
	for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(mvec[m]).beta;
	for(j=0; j<=npz-1; ++j)
	{
		z = zbeg+j*dz;
		if(z<= 0.0)
		{
			for(i=0; i<=npy-1; ++i) fld(i, j) = sqrt(improf(i));
		}
		else
		{
			if(z>=l)
			{
				for(i=0; i<=npy-1; ++i) 
					fld(i, j) = sqrt(omprof(i));
			}
			else
			{
				for(m=0; m<=num-1; ++m) 
				{
					ph = CCI*cbet(m)*(-z);
					ma(m) = pamp(m)*exp(ph);
				}
				for(i=0; i<=npy-1; ++i) 
				{
					y = ybeg+i*dy;
					f = 0.0;
					for(m=0; m<=num-1; ++m) 
					{
						for(k=0; k<=num-1; ++k)
						{
		f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im)
		    *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i));
						}
					}
					fld(i, j) = sqrt(0.5*fabs(f));
				}
			}

		}
	}
	mlout_fld(dat, npy, npz, SZ, fld);
	name[9] = 0; 
	mlout_iopropimage(name, dat); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}

/* ... as before, but with a single input mode imode only */
void EIModeArray::ioprop(const EIMode& imode, 
                           const EIModeArray& omodes, WgCrs owg,
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const 
{
	EIModeArray imodes;
	imodes.add(imode);
	Cvector iamp(1);
	iamp(0) = CC1;
	ioprop(iamp, imodes, imode.wg, 
                     omodes, owg, pert, x, l, 
                     ybeg, yend, npy, zbeg, zend, npz, ext0, ext1);
	return;
} 

/* ... as before, with single input and output modes imode and omode */
void EIModeArray::ioprop(const EIMode& imode, 
                           const EIMode& omode, 
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const
{
	EIModeArray omodes;
	omodes.add(omode);
	ioprop(imode, omodes, omode.wg, pert, x, l, 
               ybeg, yend, npy, zbeg, zend, npz, ext0, ext1);
	return;
} 

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

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

/* helper function */
double twomprodint(const EIMode& mode1, int l1, int m1, Fcomp cp1, 
	           const EIMode& mode2, int l2, int m2, Fcomp cp2,
                   Rect r)
{
	FldorDer vfod1, vfod2, hfod1, hfod2;
	double fac1, fac2;
	mode1.fldrep(cp1, fac1, vfod1, hfod1, l1, m1);
	mode2.fldrep(cp2, fac2, vfod2, hfod2, l2, m2);
	return fac1*fac2
	       *mode1.vmp.integrate(vfod1, mode2.vmp, vfod2, Interval(r.x0, r.x1))
	       *mode1.hmp.integrate(hfod1, mode2.hmp, hfod2, Interval(r.y0, r.y1));
}

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

/* scalar product 0.5\int\int(E_1x H_2y - E_1y H_2x)dxdy of two modes */
double scalprod(const EIMode& m1, const EIMode& m2) 
{
	Rect r;
	r = XYplane; 	  
/*
	int l, m;
	double s;
	if(m1.pol == TE && m2.pol == TE)
		return 0.5*m2.beta*val_invommu0(m1.wg.lambda)
		       *twomrecint(m1, EY, m2, HY, r);	

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

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

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

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

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

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

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

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

/* --- VEIMS mode solver, helper functions, mode identification ------- */

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

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

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

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

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

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


/* solve nonstandard  effective VEIMS problem,
   determine all modes with squared propagation constants between 
   bqmin and bqmax, for the problem specified by wg, pol, 
   boundaries of type bdtb at bpb, of type bdtt at bpt; 
   verb: 0: no control output, 1: progress, 2: all, 
   ec: VEIMS effective permittivity */
void findmodesX(Waveguide wg, Polarization pol, 
               Boundary_type bdtb, double bpb,
               Boundary_type bdtt, double bpt,
               double bqmin, double bqmax, ModeArray& ma, int verb, 
               Dvector ec)
{
	double k0 = val_k0(wg.lambda);
	double tbq, elim; 
	int cl, split, nsm;

// fprintf(stderr, "\nFM_I( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
// fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//                  bdtb, bpb, bdtt, bpt);
// fprintf(stderr, "      bqmin = %g, bqmax = %g )\n", bqmin, bqmax);
// verb = 2; 
// wg.write(stderr);
	
	if(bqmax < bqmin+bqSep(bqmax)) return;	
//	if(bqmax < bqmin+bqSep(bqmax)) eimodeerror("findmodesX: bqmin, bqmax");	
	cl = wg.checksymmetry();
	split = 1;
	if(cl <= 1 || cl >= wg.nx) split = 0;
	if(bdtb != bdtt) split = 0;
	if(bdtb != OPEN && fabs((wg.hx(0)-bpb)-(bpt-wg.hx(wg.nx)))>COMPTOL_HX) 
		split = 0;
	if(split == 1)
	{
		Waveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		ModeArray sm; 
		double cp = (wg.hx(cl-1)+wg.hx(cl))/2.0;
		if(verb == 1) fprintf(stderr, "s");
		if(verb == 2) fprintf(stderr, "sym ");
		sm.clear();
		findmodesX(wg0, pol, bdtb, bpb, LIMN, cp, 
		          bqmin, bqmax, sm, verb, ec); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				;
			else sm(j).expand(); 
		}
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "a");
		if(verb == 2) fprintf(stderr, "asy ");
		sm.clear();
		findmodesX(wg0, pol, bdtb, bpb, LIMD, cp, 
		          bqmin, bqmax, sm, verb, ec); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				sm(j).mirror(cp);
			else sm(j).expand(); 
		}
		ma.merge(sm);
		ma.sort();
		return;
	}	
	elim = MAXREFIND*MAXREFIND; 
	cl = -1;
	for(int l=1; l<=wg.nx; ++l)
	{ 
		double te = wg.decoupepseff(l);
		if(te < elim) 
		{
			elim = te;
			cl = l;
		}
	}
	tbq = k0*k0*elim;
	if(tbq < bqmin) tbq = bqmin;
	if(bqmax > tbq+bqSep(tbq))
	{
		Waveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		ModeArray sm;
		if(verb == 1) fprintf(stderr, "d-");
		if(verb == 2) fprintf(stderr, "d^(0, %d) ", cl);
		sm.clear();
		findmodesX(wg0, pol, bdtb, bpb, OPEN, AWAY, 
		          tbq, bqmax, sm, verb, ec); 
		for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "d+");
		if(verb == 2) fprintf(stderr, "d^(%d, %d) ", cl, wg.nx+1);
		sm.clear();
		findmodesX(wg1, pol, OPEN, -AWAY, bdtt, bpt, 
		          tbq, bqmax, sm, verb, ec); 
		for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
		ma.merge(sm);
		if(tbq > bqmin+bqSep(bqmin))
		{
			if(verb == 1) fprintf(stderr, "dr");
			if(verb == 2) fprintf(stderr, "d_(0, %d) ", wg.nx+1);
			sm.clear();
			findmodesX(wg, pol, bdtb, bpb, bdtt, bpt, 
				  bqmin, tbq, sm, verb, ec); 
			ma.merge(sm);
		}
		ma.sort();
		return;
	}
	if(bdtb != OPEN)
	{
		elim = wg.lbdecepseff(bpb); 
		tbq = k0*k0*elim;
		if(tbq < bqmin) tbq = bqmin;
		if(bqmax > tbq+bqSep(tbq))
		{
			ModeArray sm;
			if(verb == 1) fprintf(stderr, "l^");
			if(verb == 2) fprintf(stderr, "l^ ");
			sm.clear();
			findmodesX(wg, pol, OPEN, -AWAY, bdtt, bpt, 
				  tbq, bqmax, sm, verb, ec); 
			for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
			ma.merge(sm);
			if(bqmin < tbq-bqSep(tbq))
			{
				if(verb == 1) fprintf(stderr, "l_");
				if(verb == 2) fprintf(stderr, "l_ ");
				sm.clear();
				findmodesX(wg, pol, bdtb, bpb, bdtt, bpt, 
					  bqmin, tbq, sm, verb, ec); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt != OPEN)
	{
		elim = wg.ubdecepseff(bpt); 
		tbq = k0*k0*elim;
		if(tbq < bqmin) tbq = bqmin;
		if(bqmax > tbq+bqSep(tbq))
		{
			ModeArray sm;
			if(verb == 1) fprintf(stderr, "u^");
			if(verb == 2) fprintf(stderr, "u^ ");
			sm.clear();
			findmodesX(wg, pol, bdtb, bpb, OPEN, AWAY, 
				  tbq, bqmax, sm, verb, ec); 
			for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
			ma.merge(sm);
			if(bqmin < tbq-bqSep(tbq))
			{
				if(verb == 1) fprintf(stderr, "u_");
				if(verb == 2) fprintf(stderr, "u_ ");
				sm.clear();
				findmodesX(wg, pol, bdtb, bpb, bdtt, bpt, 
					  bqmin, tbq, sm, verb, ec); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt==OPEN && bdtb==OPEN)
	{
		double emin = wg.defaultepseffmin();	
		double emax = wg.defaultepseffmax();	
		tbq = k0*k0*emin;
		if(bqmin < tbq) bqmin = tbq;
		tbq = k0*k0*emax;
		if(bqmax > tbq) bqmax = tbq;
		if(bqmax < bqmin+bqSep(bqmin)) return;
	}

	double mineps = wg.n.min();
	double t = ec.min();	
	if(t < mineps) mineps = t;
	if(mineps > 0.0)
	{
		double bq, dbq;
		int zmax, zmin;
		Mode m(wg, ec, bdtb, bpb, bdtt, bpt);
		bq = bqmax-bqTol(bqmax);
		zmin = m.nummodesabove(bq);
		bq = bqmin+bqTol(bqmin);
		zmax = m.nummodesabove(bq);
		if(zmin >= zmax) return;
//		fprintf(stderr, "\nFM_F( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
//		fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//		                 bdtb, bpb, bdtt, bpt);
//		fprintf(stderr, "      bqmin = %.15g, zmax = %d, bqmax = %g, zmin= %d)\n", 
//				 bqmin, zmax, bqmax, zmin);
//		m.travresinspect(bqmin+bqTol(bqmin), bqmax-bqTol(bqmax));

		Ivector found(zmax);
		found.init(0);	
		Dvector bqf(zmax);
		bqf.init(0.0);

		trap(ma, m, bqmin, bqmax, 2.0, 0, 1, bqf, found, verb);
		
		int fa, mn, i, newv;
		double bqup, bqdn, sf;
		fa = 0;
		while(fa == 0)
		{
			fa = 1; 
			mn = zmax-1;
			while(fa == 1 && mn>= zmin)
			{
				if(found(mn) != 0) --mn;
				else
				{
					fa = 0;
					i=mn+1; 
					while(i<zmax && found(i) == 0) ++i;
					if(i >= zmax) bqdn = bqmin;
					else          bqdn = bqf(i);
					i=mn-1;
					while(i>zmin-1 && found(i) == 0) --i;
					if(i <= zmin-1) bqup = bqmax;
					else 	  bqup = bqf(i);
					dbq = bqup-bqdn;
					sf = 1.0;
					while((newv = trap(ma, m, 
						      bqdn+bqTol(bqdn), 
						      bqup-bqTol(bqup), sf, 1, 1,
						      bqf, found, verb)) == 0 
					      && sf < 4.0) sf += 0.5;
					if(newv >= 1) goto fnd;
					if(mn < zmax-1 && found(mn+1) != 0)
					{
					     bqdn = bqf(mn+1);
					     sf = 1.0;
					     tbq = bqdn+dbq/100.0;
					     while((newv = trap(ma, m, bqdn, tbq, sf, 
								 1, 1, bqf, found, verb)
						    ) == 0 && sf < 4.0) sf += 0.5;
					}
					if(newv >= 1) goto fnd;
					if(mn > zmin && found(mn-1) != 0)
					{
					     bqup = bqf(mn-1);
					     sf = 1.0;
					     tbq = bqup-dbq/100.0;
					     while((newv = trap(ma, m, tbq, bqup, sf,
								 1, 1, bqf, found, verb)
						    ) == 0 && sf < 4.0) sf += 0.5;
					}
					if(newv >= 1) goto fnd;

					if(mn == zmax-1)
					{
					     bqdn = bqmin;
					     sf = 1.0;
					     tbq = bqdn+dbq/100.0;
					     while((newv = trap(ma, m, bqdn, tbq, sf, 
								 1, 1, bqf, found, verb)
						    ) == 0 && sf < 4.0) sf += 0.5;
					}
					if(newv >= 1) goto fnd;
					if(mn == zmin)
					{
					     bqup = bqmax;
					     sf = 1.0;
					     tbq = bqup-dbq/100.0;
					     while((newv = trap(ma, m, tbq, bqup, sf,
								 1, 1, bqf, found, verb)
						    ) == 0 && sf < 4.0) sf += 0.5;
					}
					if(newv >= 1) goto fnd;

					fprintf(stderr, " T%c%d ?\n", 
							polCHR(m.pol), mn);
					m.travresinspect(bqmin, bqmax);
					m.wg.write(stderr);
					m.wg.plot('e','r');
					for(int j=0; j<=ma.num-1; ++j) 
					{
						fprintf(stderr, " < %s > ", ma(j).ids);
						ma(j).writeprofile(dig10(j), dig1(j));
					}
					eimodeerror("findmodesX, pos: missing mode");
fnd: 					;;
				}
			}
		}
	}
	else    // min epsilon negative 
	{
		Mode m(wg, ec, bdtb, bpb, bdtt, bpt);
		if(bqmin >= bqmax) return;
//		fprintf(stderr, "\nFM_F( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
//		fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//		                 bdtb, bpb, bdtt, bpt);
//		fprintf(stderr, "      bqmin = %.15g, bqmax = %g)\n", 
//				 bqmin, bqmax);
//		m.travresinspect(bqmin+bqTol(bqmin), bqmax-bqTol(bqmax));

		trapX(ma, m, bqmin, bqmax, 3.0, 0, 1, verb);

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

/* --- VEIMS mode solver ---------------------------------------------- */

/* VEIMS algorithm: solve the effective 1D problem */
int veims_hsolve(Waveguide hwg, Dvector hwg_ec, double eemin, double eemax, ModeArray& ma)         
{
	double emin = hwg.defaultepseffmin();	
	double emax = hwg.defaultepseffmax();	
	if(emin < 0.0) emin = 0.0;
	if(emin < eemin) emin = eemin;
	if(emax > eemax) emax = eemax;
	double k0 = val_k0(hwg.lambda);
	double bqmin = k0*k0*emin;
	double bqmax = k0*k0*emax;

	bqmin += SLAMS_BQSEP;
	ma.clear();
	if(bqmin >= bqmax) return 0;
	findmodesX(hwg, TM, OPEN, -AWAY, OPEN, AWAY, bqmin, bqmax, ma, 0, hwg_ec);
//	findmodesX(hwg, TM, OPEN, -AWAY, OPEN, AWAY, bqmin, bqmax, ma, 2, hwg_ec);
	ma.sort();
	for(int j=0; j<=ma.num-1; ++j)
	{
		if(ma(j).special==0 && ma(j).ord!=j) 
			eimodeerror("veims_hsolve: mode set corrupted");
		ma(j).setids(j);
	}
	return ma.num;
}


/* VEIMS guided mode analysis, returns the number of found modes
   wg, st: the structure under consideration
   pol: polarization type, TE or TM
   rsn: nmber of the reference slice; -1: determined automatically
   ma: the modes found (output)
   quiet == 1: suppress log output */
int veims(SegWgCrs st, Polarization pol, int rsn, EIModeArray& ma, int quiet)
{
	st.consistency();
	double lambda = st(0).lambda;
	double k0 = val_k0(lambda);
	WgCrs wg = st.wgcrs();
	Waveguide rs;
	int ri = 0;
	double mn = 0.0;
	int nxp;

	if(rsn < 0)
	{
		for(int j=1; j<=wg.ny; ++j) 
		{
			ModeArray a;
			nxp = modeanalysis(st(j), pol, a, 1);
			if(nxp >= 1)
			{
				if(a(0).neff > mn)	
				{
					ri = j;		
					mn = a(0).neff;
				}
			}
		}
		if(ri<=0 || ri>=wg.ny+1)
		{ 
//			eimodeerror("solver: reference slice identification");
			ma.clear();
			return 0;
		}
	}
	else
	{ 
		ri = rsn; 
	}
	rs = st(ri);

	ModeArray xpa;
	nxp = modeanalysis(rs, pol, xpa, 1);

if(quiet != 1)
{
	fprintf(stderr, 
		"\n------------- Metric - VEIMS --------------- 2010 ---\n");
	switch(pol)
	{
		case TE: fprintf(stderr, "TE, "); break;
		case TM: fprintf(stderr, "TM, "); break;
	}
	fprintf(stderr, "Lambda: %.10g, K0: %g, ", lambda, k0);
	fprintf(stderr, "RSlice: %d,  #VModes: %d\n", ri, nxp);
	fprintf(stderr, 
		"-----------------------------------------------------\n");
	fprintf(stderr, "    Nx: %d  Ny: %d\n", wg.nx, wg.ny);
	fprintf(stderr, "    Hx: ");
	for(int j=0; j<=wg.nx; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
	fprintf(stderr, "\n");
	fprintf(stderr, "    Hy: ");
	for(int j=0; j<=wg.ny; ++j) fprintf(stderr, "%6.10g  ", wg.hy(j));
	fprintf(stderr, "\n");
	fprintf(stderr, "     N: ");
	for(int j=wg.nx+1; j>=0; --j) 
	{
		for(int k=0; k<=wg.ny+1; ++k) fprintf(stderr, "%6.10g  ", wg.n(j,k));
		if(j>0) fprintf(stderr, "\n        ");
	}
	fprintf(stderr, "\n");
	fprintf(stderr, 
		"-----------------------------------------------------\n");
}
	ModeArray tma, ttma;
	ttma.clear();
	modeanalysis(st(0),       TE, tma, 1); ttma.merge(tma);
	modeanalysis(st(0),       TM, tma, 1); ttma.merge(tma);
	modeanalysis(st(st.ny+1), TE, tma, 1); ttma.merge(tma);
	modeanalysis(st(st.ny+1), TM, tma, 1); ttma.merge(tma);
	ttma.sort();
	double leakylim = -1.0; 
	if(ttma.num >= 1) leakylim = ttma(0).neff;
	tma.clear(); ttma.clear();

	ma.clear();
	for(int xpi=0; xpi<=nxp-1; ++xpi)	
	{
		if(quiet != 1) fprintf(stderr, "v(%d):\n", xpi);
		EIMode m;
		m.wg = wg;
		m.st = st;
		m.pol = pol;
		m.k0 = k0;
		m.rsegi = ri;
		Mode xp = xpa(xpi);
		m.vmp = xp; 

		Waveguide hwg(wg.ny);
		hwg.hx = wg.hy;
		hwg.lambda = wg.lambda;
		hwg.special = 1;
		Dvector hwg_q(wg.ny+2);
		Dvector hwg_ec(wg.ny+2);
		for(int s=0; s<=wg.ny+1; ++s)
		{
			if(rs.nx > wg.nx) eimodeerror("veims: nx");
			if(pol == TE)
			{
				double na = 0.0; 
				double no = 0.0;
				for(int l=0; l<=wg.nx+1; ++l)
				{ 
					double x0, x1, xc;
					if(l==0) x0 = -AWAY; else x0 = wg.hx(l-1);
					if(l==wg.nx+1) x1 = AWAY; else x1 = wg.hx(l);
					xc = 0.5*(x0+x1);
					double i = xp.integrate(FLD, FLD, Interval(x0, x1));
					no += i;
					na += i*(wg.eps(l, s) - rs.eps(xc));
				}
				hwg.n(s) = xp.neff*xp.neff+na/no;
				hwg_q(s) = 1.0;
				hwg_ec(s) = hwg.n(s);
			}
			else // TM
			{
				double x_1_x      = 0.0;
				double x_1de_x    = 0.0;
				double x_edeq_x   = 0.0;
				double xp_1de_xp  = 0.0;
				double xp_edeq_xp = 0.0;
				
				for(int l=0; l<=wg.nx+1; ++l)
				{ 
					double x0, x1, xc;
					if(l==0) x0 = -AWAY; else x0 = wg.hx(l-1);
					if(l==wg.nx+1) x1 = AWAY; else x1 = wg.hx(l);
					xc = 0.5*(x0+x1);
					double i;
					i = xp.integrate(FLD, FLD, Interval(x0, x1));
					x_1_x      += i;
					x_1de_x    += i/rs.eps(xc);
					x_edeq_x   += i*wg.eps(l, s)/rs.eps(xc)/rs.eps(xc);
					i = xp.integrate(DER, DER, Interval(x0, x1));
					xp_1de_xp  += i/rs.eps(xc);
					xp_edeq_xp += i*wg.eps(l, s)/rs.eps(xc)/rs.eps(xc);
				}
				hwg_q(s) = xp_1de_xp/xp_edeq_xp;
				hwg_ec(s) = xp.neff*xp.neff*hwg_q(s)+x_1_x/x_1de_x*(1.0-hwg_q(s)) ;
				hwg.n(s) = hwg_ec(s)*x_edeq_x/x_1de_x;
			}
		}
		m.hwg = hwg;
// 		if(quiet != 1) hwg.write(stderr);
		m.hwg_q = hwg_q;
		m.hwg_ec = hwg_ec;
		if(quiet != 1) 
		{
			if(pol == TE)
			{
				fprintf(stderr, "eps_eff:  ");
				for(int s=0; s<=wg.ny+1; ++s) fprintf(stderr, "%6.4g  ", hwg.eps(s));
				fprintf(stderr, "\n");
			}
			else
			{
				fprintf(stderr, "eps_eff:  ");
				for(int s=0; s<=wg.ny+1; ++s) fprintf(stderr, "%6.4g  ", hwg.eps(s));
				fprintf(stderr, "\n");
				fprintf(stderr, "eps_c:    ");
				for(int s=0; s<=wg.ny+1; ++s) fprintf(stderr, "%6.4g  ", hwg_ec(s));
				fprintf(stderr, "\n");
				fprintf(stderr, "q:        ");
				for(int s=0; s<=wg.ny+1; ++s) fprintf(stderr, "%6.4g  ", hwg_q(s));
				fprintf(stderr, "\n");
			}
		}
		
		double eemin = xp.wg.defaultepseffmin();
		double eemax = xp.wg.defaultepseffmax();
		ModeArray hma;
		veims_hsolve(hwg, hwg_ec, eemin, eemax, hma); 
		for(int j=0; j<=hma.num-1; ++j)
		{
			if(hma(j).neff*hma(j).neff > eemin)
			{
				m.beta = hma(j).beta;	
				m.neff = hma(j).neff;
				m.hmp = hma(j);
				m.normalize(1.0);
				m.setfieldmax();
				m.leaky = 0;
				if(m.neff <= leakylim) m.leaky = 1;
				ma.add(m);
				if(quiet != 1)
				{
					fprintf(stderr, "T%c_(%d, %d), neff = %g", 
						polCHR(pol), xpi, j, m.neff);
					if(m.leaky == 1) fprintf(stderr, " (leaky !)");
					fprintf(stderr, "\n");
				}
				
			}
		}
		if(quiet != 1) fprintf(stderr, "   - - -\n");
	}
	if(quiet != 1) fprintf(stderr, "\n");
	return ma.num;	
}
// with automatically determined number of the reference slice
int veims(SegWgCrs st, Polarization pol, EIModeArray& ma, int quiet)
{
	return veims(st, pol, -1, ma, quiet);
}
int veims(SegWgCrs st, Polarization pol, EIModeArray& ma)
{
        return veims(st, pol, -1, ma, 0);
}
int veims(WgCrs wg, Polarization pol, EIModeArray& ma, int quiet)
{
	return veims(SegWgCrs(wg), pol, -1, ma, quiet);
}
int veims(WgCrs wg, Polarization pol, EIModeArray& ma)
{
        return veims(wg, pol, ma, 0);
}