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

/*
 * fbmode.cpp
 * Floquet-Bloch-mode analysis of rectangular periodic dielectric waveguides
 * 2-D solver based on BEP / QUEP expansions with Dirichlet boundary conditions
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"slamode.h"
#include"slamarr.h"
#include"slams.h"
#include"bepfld.h"    
#include"quepfld.h"    
#include"matlvis.h"    
#include"fbmode.h"

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

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

/* map propagation coordinate to original period */ 
double FBMode::wrap(double z) const  
{
	if(z > st.hz(st.nz)) return z-floor((z-st.hz(0))/L)*L;
	if(z < st.hz(0)) return z+floor((st.hz(st.nz)-z)/L)*L;
	return z; 
}

/* FB mode profile at position x, z, 
   cp: EX, EY, EZ, HX, HY, HZ */
Complex FBMode::fbmprof(Fcomp cp, double x, double z) const
{
	double zw = wrap(z);
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ:
		if(type == 'B') return fldB.field(cp, x, zw)*expi(beta*(zw-st.hz(0)));
		else            return fldQ.field(cp, x, zw)*expi(beta*(zw-st.hz(0)));
		break;
	default:
		fbmodeerror("fbmprof: cp");
		break;
	}
	return CC0;
}

/* e.m. field at point x, z 
   cp: EX, EY, EZ, HX, HY, HZ */
Complex FBMode::field(Fcomp cp, double x, double z) const
{
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ:
		if(st.hz(0) <= z && z <= st.hz(st.nz))
		{
			if(type == 'B') return fldB.field(cp, x, z);
			else            return fldQ.field(cp, x, z);
		}
		return fbmprof(cp, x, z)*expi(-beta*(z-st.hz(0)));
		break;
	default:
		fbmodeerror("field: cp not implemented");
		break;
	}
	return CC0;
}


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

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

/* ... np equidistant field values on the straight line 
   between (x0, z0) and (x1, z1) */
Cvector FBMode::field(Fcomp cp, double x0, double z0, 
                                double x1, double z1, int np) const
{
	if(np <= 1) fbmodeerror("discretized field: np"); 
	Cvector fld(np);
	double dx = (x1-x0)/(np-1);
	double dz = (z1-z0)/(np-1);
	double x = x0; 
	double z = z0;
	for(int j=0; j<=np-1; ++j)
	{
		fld(j) = field(cp, x, z);
		x += dx;
		z += dz;
	}
	return fld;
}
Dvector FBMode::field(Fcomp cp, Afo foa, double x0, double z0, 
                                         double x1, double z1, int np) const
{
	Cvector f = field(cp, x0, z0, x1, z1, np);
	return realize(f, foa);
}

/* optical power (per lateral unit length) carried by the mode,
   x-integral of Sz evaluated at the first interface,
   +1 / -1 for normalized modes */
double FBMode::power() const
{
	if(type == 'B') return fldB.Pout(SBRIGHT)-fldB.Pout(SBLEFT);
	else            return fldQ.Pout(RIGHT)-fldQ.Pout(LEFT);
}

/* polarization character of the mode, relevant for vM = 1 only, 
   relative incoming / outgoing power carried by TE fields;
   return value > 0.5 / < 0.5: a more TE-like / more TM-like mode */
double FBMode::polcharacter() const
{
	if(type == 'B') return fabs(fldB.Pout(SBRIGHT, TE)-fldB.Pout(SBLEFT, TE))
	                       / fabs(fldB.Pout(SBRIGHT)-fldB.Pout(SBLEFT));
	return fabs(fldQ.Pout(RIGHT, TE)-fldQ.Pout(LEFT, TE))
	        / fabs(fldQ.Pout(RIGHT)-fldQ.Pout(LEFT));
}


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

/* preparation of physical field plots:
   selection of the global phase / snapshot time 0 within one period */

/* multiply the entire field by a phase factor exp(i phi) */
void FBMode::adjustphase(double phi)
{
	if(type == 'B') fldB.adjustphase(phi); 
	else            fldQ.adjustphase(phi); 
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   basic field component (cp = EY for TE, cp = HY for TM) exhibits the 
   maximum field at position (x, z) */
void FBMode::adjustphase(double x, double z)
{
	if(type == 'B') fldB.adjustphase(x, z); 
	else            fldQ.adjustphase(x, z); 
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   basic field component (cp = EY for TE, cp = HY for TM) exhibits an 
   overall field maximum (-> time snapshot of an extremal state  
   in case of a configuration with partially standing waves)   
   ... determine the maximum on a regular mesh of numx * mumz points
       on a reasonable computational window */
void FBMode::adjustphase(int numx, int numz)
{
	if(type == 'B') fldB.adjustphase(numx, numz);
	else            fldQ.adjustphase(numx, numz);
	return;
}

/* ... as before, with (coarse) defaults for numx, numz */
void FBMode::adjustphase()
{
	if(type == 'B') fldB.adjustphase();
	else            fldQ.adjustphase();
	return;
}

/* ------------------------------------------------------------------- */
/* visualization of the light propagation, output: MATLAB scripts */

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

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

/* mode interference plots, fancy style surface ;-)
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void FBMode::fplot(Fcomp cp, 
                   double xbeg, double xend, double zbeg, double zend,
                   int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "pl_____F.m";
	name[2] = 'f';
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "fplot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	SegWgStruct tst = st.unwrap(zbeg, zend);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
        mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Interval dwx(xbeg, xend);
	Interval dwz(zbeg, zend);
	Dmatrix fld = field(cp, ORG, dwx, npx, dwz, npz);
	mlout_fld(dat, npx, npz, cp, fld);
	Circuit stc = tst.circuit();
	int nsm = 2*(stc.nz+2)*(stc.nx+1)+2*(stc.nx+2)*(stc.nz+1)+4;
	Dvector bx0, bx1, bz0, bz1;
	bx0 = bx1 = bz0 = bz1 = Dvector(nsm);
	Ivector bsi, bnp;
	bsi = bnp = Ivector(nsm);
	Dvector bf(2*(2*(stc.nz+2)*npx+2*(stc.nx+2)*npz)+2*npx+2*npz);
	double x0, x1, z0, z1, xp, zp;
	Dvector fv;
	int nc = 0;
	int si = 0;
	for(int l=0; l<=stc.nx+1; ++l)
	{
		if(l==0)
		{
			if(xbeg < stc.hx(0)) x0 = xbeg; 	
			else                x0 = stc.hx(0)-COMPTOL_HX;
		}
                else x0 = stc.hx(l-1);
		if(l==stc.nx+1)
		{
			if(xend > stc.hx(stc.nx)) x1 = xend; 	
			else                    x1 = stc.hx(stc.nx)+COMPTOL_HX; 	
		}
                else x1 = stc.hx(l);
		Interval r(x0, x1);
		int o = dwx.clip(r);
		if(o)
		{
		    for(int s=0; s<=stc.nz; ++s)
		    {
			if(fabs(stc.n(l,s)-stc.n(l,s+1)) > COMPTOL_N)
			{
			    zp = stc.hz(s);
			    if(dwz.in(zp))
			    {
				int np=(int)(((double) npx)*r.len()/dwx.len());
				if(np >= 2)
				{
				fv = field(cp, ORG,
			  	           r.x0+HDIST, zp-HDIST, 
				           r.x1-HDIST, zp-HDIST, np); 
				bx0(nc) = r.x0+HDIST; bz0(nc) = zp-HDIST;   
				bx1(nc) = r.x1-HDIST; bz1(nc) = zp-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				fv = field(cp, ORG,
			  	           r.x0+HDIST, zp+HDIST, 
				           r.x1-HDIST, zp+HDIST, np); 
				bx0(nc) = r.x0+HDIST; bz0(nc) = zp+HDIST;   
				bx1(nc) = r.x1-HDIST; bz1(nc) = zp+HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				}
			    }
			}
		    }
		}
	}
	for(int s=0; s<=stc.nz+1; ++s)
	{
		if(s==0)
		{
			if(zbeg < stc.hz(0)) z0 = zbeg; 	
			else                z0 = stc.hz(0)-COMPTOL_HX;
		}
                else z0 = stc.hz(s-1);
		if(s==stc.nz+1)
		{
			if(zend > stc.hz(stc.nz)) z1 = zend; 	
			else                    z1 = stc.hz(stc.nz)+COMPTOL_HX; 	
		}
                else z1 = stc.hz(s);
		Interval r(z0, z1);
		int o = dwz.clip(r);
		if(o)
		{
		    for(int l=0; l<=stc.nx; ++l)
		    {
			if(fabs(stc.n(l,s)-stc.n(l+1,s)) > COMPTOL_N)
			{
			    xp = stc.hx(l);
			    if(dwx.in(xp))
			    {
				int np=(int)(((double)npz)*r.len()/dwz.len()); 
				if(np >= 2)
				{
				fv = field(cp, ORG, 
			  	           xp-HDIST, r.x0+HDIST, 
				           xp-HDIST, r.x1-HDIST, np); 
				bx0(nc) = xp-HDIST; bz0(nc) = r.x0+HDIST;   
				bx1(nc) = xp-HDIST; bz1(nc) = r.x1-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				fv = field(cp, ORG,
			  	           xp+HDIST, r.x0+HDIST, 
				           xp+HDIST, r.x1-HDIST, np); 
				bx0(nc) = xp+HDIST; bz0(nc) = r.x0+HDIST;   
				bx1(nc) = xp+HDIST; bz1(nc) = r.x1-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				}
			    }
			}
		    }
		}
	}
	fv = field(cp, ORG, xbeg, zbeg, xend, zbeg, npx);
	bx0(nc) = xbeg; bz0(nc) = zbeg; bx1(nc) = xend; bz1(nc) = zbeg;
	bnp(nc) = npx; bsi(nc) = si; bf.setsubvector(fv, si); si += npx;
	++nc;
	fv = field(cp, ORG, xbeg, zend, xend, zend, npx);
	bx0(nc) = xbeg; bz0(nc) = zend; bx1(nc) = xend; bz1(nc) = zend;
	bnp(nc) = npx; bsi(nc) = si; bf.setsubvector(fv, si); si += npx;
	++nc;
	fv = field(cp, ORG, xbeg, zbeg, xbeg, zend, npz);
	bx0(nc) = xbeg; bz0(nc) = zbeg; bx1(nc) = xbeg; bz1(nc) = zend;
	bnp(nc) = npz; bsi(nc) = si; bf.setsubvector(fv, si); si += npz;
	++nc;
	fv = field(cp, ORG, xend, zbeg, xend, zend, npz);
	bx0(nc) = xend; bz0(nc) = zbeg; bx1(nc) = xend; bz1(nc) = zend;
	bnp(nc) = npz; bsi(nc) = si; bf.setsubvector(fv, si); si += npz;
	++nc;
	mlout_bddata(dat, cp, nc, bsi, bnp, bf, bx0, bz0, bx1, bz1);
	name[8] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_fancy(name, desc, dat, cp, nc);  	
	mlout_print(dat, name, 'p');
	fclose(dat);
	return;
}

/* animation of the light propagation
   the frames show the EY component (TE) resp. HY component (TM), at
   ntfr equidistant times over one time period
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void FBMode::movie(double xbeg, double xend, double zbeg, double zend,
                   int npx, int npz, int ntfr, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "an____M.m";
	double tmax, dt;
	double lambda = st(0).lambda;
	Fcomp cp = principalcomp(pol);
        name[2] = fldchr(cp);
        name[3] = cpchr(cp);      
	name[4] = ext0;
	name[5] = ext1;
	tmax = val_period_T(lambda);
	dt = tmax/ntfr;
	fprintf(stderr, 
	   "movie(%d x [%g, %g] x [%g, %g]), T=%g fs, dt=%g fs >> %s\n", 
	   ntfr, xbeg, xend, zbeg, zend, tmax, dt, name);
	SegWgStruct tst = st.unwrap(zbeg, zend);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, tst, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix fld = field(cp, Interval(xbeg, xend), npx, 
	                        Interval(zbeg, zend), npz); 
	double famp = fld.max();
	mlout_fld(dat, npx, npz, cp, fld.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npz, cp, fld.im());
	mlout_fldtoim(dat, cp);
	name[7] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_genmovie(cp, name, desc, dat, ntfr, dt, val_omega(lambda), famp); 
	fclose(dat);
	return;
}
		
/* animation of the light propagation, fancy style ;-)
   the frames show the EY component (TE) resp. HY component (TM), at
   ntfr equidistant times over one time period
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void FBMode::fmovie(double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, int ntfr, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "fa____M.m";
	double tmax, dt;
	double lambda = st(0).lambda;
	Fcomp cp = principalcomp(pol);
        name[2] = fldchr(cp);
        name[3] = cpchr(cp);      
	name[4] = ext0;
	name[5] = ext1;
	tmax = val_period_T(lambda);
	dt = tmax/ntfr;
	fprintf(stderr, 
	   "fmovie(%d x [%g, %g] x [%g, %g]), T=%g fs, dt=%g fs >> %s\n", 
	   ntfr, xbeg, xend, zbeg, zend, tmax, dt, name);
	SegWgStruct tst = st.unwrap(zbeg, zend);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, tst, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Interval dwx(xbeg, xend);
	Interval dwz(zbeg, zend);
	Cmatrix fld = field(cp, dwx, npx, dwz, npz); 
	double famp = fld.max();
	mlout_fld(dat, npx, npz, cp, fld.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npz, cp, fld.im());
	mlout_fldtoim(dat, cp);
	Circuit stc = tst.circuit();
	int nsm = 2*(stc.nz+2)*(stc.nx+1)+2*(stc.nx+2)*(stc.nz+1)+4;
	Dvector bx0, bx1, bz0, bz1;
	bx0 = bx1 = bz0 = bz1 = Dvector(nsm);
	Ivector bsi, bnp;
	bsi = bnp = Ivector(nsm);
	Cvector bf(2*(2*(stc.nz+2)*npx+2*(stc.nx+2)*npz)+2*npx+2*npz);
	double x0, x1, z0, z1, xp, zp;
	Cvector fv;
	int nc = 0;
	int si = 0;
	for(int l=0; l<=stc.nx+1; ++l)
	{
		if(l==0)
		{
			if(xbeg < stc.hx(0)) x0 = xbeg; 	
			else                x0 = stc.hx(0)-COMPTOL_HX;
		}
                else x0 = stc.hx(l-1);
		if(l==stc.nx+1)
		{
			if(xend > stc.hx(stc.nx)) x1 = xend; 	
			else                    x1 = stc.hx(stc.nx)+COMPTOL_HX; 	
		}
                else x1 = stc.hx(l);
		Interval r(x0, x1);
		int o = dwx.clip(r);
		if(o)
		{
		    for(int s=0; s<=stc.nz; ++s)
		    {
			if(fabs(stc.n(l,s)-stc.n(l,s+1)) > COMPTOL_N)
			{
			    zp = stc.hz(s);
			    if(dwz.in(zp))
			    {
				int np=(int)(((double) npx)*r.len()/dwx.len());
				if(np >= 2)
				{
				fv = field(cp, 
			  	           r.x0+HDIST, zp-HDIST, 
				           r.x1-HDIST, zp-HDIST, np); 
				bx0(nc) = r.x0+HDIST; bz0(nc) = zp-HDIST;   
				bx1(nc) = r.x1-HDIST; bz1(nc) = zp-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				fv = field(cp,
			  	           r.x0+HDIST, zp+HDIST, 
				           r.x1-HDIST, zp+HDIST, np); 
				bx0(nc) = r.x0+HDIST; bz0(nc) = zp+HDIST;   
				bx1(nc) = r.x1-HDIST; bz1(nc) = zp+HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				}
			    }
			}
		    }
		}
	}
	for(int s=0; s<=stc.nz+1; ++s)
	{
		if(s==0)
		{
			if(zbeg < stc.hz(0)) z0 = zbeg; 	
			else                z0 = stc.hz(0)-COMPTOL_HX;
		}
                else z0 = stc.hz(s-1);
		if(s==stc.nz+1)
		{
			if(zend > stc.hz(stc.nz)) z1 = zend; 	
			else                    z1 = stc.hz(stc.nz)+COMPTOL_HX; 	
		}
                else z1 = stc.hz(s);
		Interval r(z0, z1);
		int o = dwz.clip(r);
		if(o)
		{
		    for(int l=0; l<=stc.nx; ++l)
		    {
			if(fabs(stc.n(l,s)-stc.n(l+1,s)) > COMPTOL_N)
			{
			    xp = stc.hx(l);
			    if(dwx.in(xp))
			    {
				int np=(int)(((double)npz)*r.len()/dwz.len()); 
				if(np >= 2)
				{
				fv = field(cp, 
			  	           xp-HDIST, r.x0+HDIST, 
				           xp-HDIST, r.x1-HDIST, np); 
				bx0(nc) = xp-HDIST; bz0(nc) = r.x0+HDIST;   
				bx1(nc) = xp-HDIST; bz1(nc) = r.x1-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				fv = field(cp,
			  	           xp+HDIST, r.x0+HDIST, 
				           xp+HDIST, r.x1-HDIST, np); 
				bx0(nc) = xp+HDIST; bz0(nc) = r.x0+HDIST;   
				bx1(nc) = xp+HDIST; bz1(nc) = r.x1-HDIST;
				bnp(nc) = np; bsi(nc) = si; 
				bf.setsubvector(fv, si); si += np;
				++nc;
				}
			    }
			}
		    }
		}
	}
	fv = field(cp, xbeg, zbeg, xend, zbeg, npx);
	bx0(nc) = xbeg; bz0(nc) = zbeg; bx1(nc) = xend; bz1(nc) = zbeg;
	bnp(nc) = npx; bsi(nc) = si; bf.setsubvector(fv, si); si += npx;
	++nc;
	fv = field(cp, xbeg, zend, xend, zend, npx);
	bx0(nc) = xbeg; bz0(nc) = zend; bx1(nc) = xend; bz1(nc) = zend;
	bnp(nc) = npx; bsi(nc) = si; bf.setsubvector(fv, si); si += npx;
	++nc;
	fv = field(cp, xbeg, zbeg, xbeg, zend, npz);
	bx0(nc) = xbeg; bz0(nc) = zbeg; bx1(nc) = xbeg; bz1(nc) = zend;
	bnp(nc) = npz; bsi(nc) = si; bf.setsubvector(fv, si); si += npz;
	++nc;
	fv = field(cp, xend, zbeg, xend, zend, npz);
	bx0(nc) = xend; bz0(nc) = zbeg; bx1(nc) = xend; bz1(nc) = zend;
	bnp(nc) = npz; bsi(nc) = si; bf.setsubvector(fv, si); si += npz;
	++nc;
	mlout_bddata(dat, cp, nc, bsi, bnp, bf, bx0, bz0, bx1, bz1);
	name[7] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_genfmovie(cp, name, desc, dat, nc, ntfr, dt, val_omega(lambda), famp);
	fclose(dat);
	return;
}

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

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

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

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

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

}

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

}

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

/* add a mode */
void FBModeArray::add(FBMode m)
{
	FBMode *avec;
	avec = new FBMode[num+1];
	FBMode* ap = avec;
	FBMode* 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 FBModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		fbmodeerror("FBModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	FBMode *avec;
	avec = new FBMode[num-1];
	FBMode* ap = avec;
	FBMode* 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 FBModeArray nma */ 
void FBModeArray::merge(FBModeArray& ma)
{
	FBMode *avec;
	avec = new FBMode[num+ma.num];
	FBMode* ap = avec;
	FBMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	FBMode* 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 FBModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	FBMode 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;
}

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

/* spectral reduction to avoid eigensolver failures:
   lowest number of spectral terms allowed */
int FBNUMMXMIN = 20;

/* filter erroneous modal solutions with large power attenuation
   due to a nonnegligible imaginary part in the exp(i*beta*L) eigenvalue,
   att = |exp(i*beta*L)|, 
   discard modes where  where |att-1| > FBMAXATTLVL */
double FBMAXATTLVL = 0.02; 

/* selection of valid FB modes:
   after projection onto "vertical" directional slab modes associated 
   with the lower- and uppermost layers in the structures, with LIMD and LIMN
   boundary conditions on one period, FB modes are selected as valid 
   if the outward modes in those sets are assigned a relative directional 
   power below FBVOPLVL;
   FBOPLVL = 0.0: only modes "below the light line" are considered */
double FBVOPLVL = 1.0e-4; 

/* FB-mode solver, guided mode analysis,
   st: the structure under consideration, including wavelength
   pol: polarization type, TE or TM
   cw: computational window for the BEP expansion
   Mx: number of spectral discetization terms
   red = 0, 1: prevent, permit spectral reduction to avoid convergence
               failures of the eigenvalue solver; limit: Mx >= NUMMXMIN 
   ma: the FB modes found (output)
   quiet == 1: suppress log output
   return value: >=0: the number of found modes,
                  <0: failed */
int fbmsolve(SegWgStruct st, Polarization pol, Interval cw, int Mx, int red,
             FBModeArray& ma, int quiet)
{
	if(cw.x0 > cw.x1) fbmodeerror("fbmsolve: computational window");
	double lambda = st(0).lambda;
	double k0 = val_k0(lambda);
	if(st.nz < 3) fbmodeerror("fbmsolve: illegal structure (1)");  
	if(st(0).equal(st(1)) != 1) fbmodeerror("fbmsolve: illegal structure (2)");  
	if(st(0).equal(st(st.nz)) != 1) fbmodeerror("fbmsolve: illegal structure (3)");  
	if(st(0).equal(st(st.nz+1)) != 1) fbmodeerror("fbmsolve: illegal structure (4)");  

	int done = 0;
	int nummx = Mx;
	if(nummx < FBNUMMXMIN) fbmodeerror("fbmsolve: Mx < FBNUMMXMIN");  
	double e_min = 1.0;
	for(int j=0; j<=st.nz+1; ++j)
	{
		if(st(j).defaultepseffmin() > e_min) 
			e_min = st(j).defaultepseffmin();
	}	
	double n_min = sqrt(e_min);
	double L = st.hz(st.nz)-st.hz(0);
	
	if(quiet != 1)
	{
		Circuit rc = st.circuit();
fprintf(stderr, "\n------------- FBMode, BEP ------------------------------------------ '11 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d\n", cw.x0, cw.x1, Mx);
fprintf(stderr,   "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", rc.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=rc.nx; ++j) fprintf(stderr, "%6.10g  ", rc.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", rc.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=rc.nz; ++j) fprintf(stderr, "%6.10g  ", rc.hz(j));
fprintf(stderr, "\n");
if(rc.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=rc.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=rc.nz+1; ++j) fprintf(stderr, "%6.10g  ", rc.n(i, j));
fprintf(stderr,"\n");
}
fprintf(stderr, "L: %g  ", L);
fprintf(stderr, "lambda: %.10g  k0: %g", rc.lambda, k0);
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
if(FBVOPLVL <= 0.0)
{
	fprintf(stderr, "guided FBModes: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g\n", L*n_min/lambda); 
}
else
{
	fprintf(stderr, "guided FBModes: rel. vert. Output < %g\n", FBVOPLVL); 
	fprintf(stderr, "(light cone: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g)\n", L*n_min/lambda); 
}
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	while(done == 0 && nummx >= FBNUMMXMIN)
	{
		
		BepField bfld(st, 0, pol, cw, nummx, 1);
		if(bfld.s.nz < 3) fbmodeerror("fbmsolve: BEP initialization"); 
		Cmatrix M0mbf, M0mbb, M0pff, M0pfb;
		Cmatrix MNmbf, MNmbb, MNpff, MNpfb;
		Cmatrix Kff, Kfb, Kbf, Kbb;
		Cmatrix *Iff, *Ifb, *Ibf, *Ibb;
		Cvector *TX;
		Iff = new Cmatrix[bfld.s.nz];
		Ifb = new Cmatrix[bfld.s.nz];
		Ibf = new Cmatrix[bfld.s.nz];
		Ibb = new Cmatrix[bfld.s.nz];        
		TX  = new Cvector[bfld.s.nz];        
		Cvector T1, TN;
		Cmatrix D0, D1, DN, DNp;
		bfld.quepsetup(M0mbf, M0mbb, M0pff, M0pfb, T1, D0, D1,
			    MNmbf, MNmbb, MNpff, MNpfb, TN, DN, DNp,
			    Kff, Kfb, Kbf, Kbb, Iff, Ifb, Ibf, Ibb, TX);     

		Cmatrix kff, kfb, kbf, kbb;
		kff = diagmult(TN, diagmult(Kff, T1));
		kbf = diagmult(T1, diagmult(Kbf, T1));
		kfb = diagmult(TN, diagmult(Kfb, TN));
		kbb = diagmult(T1, diagmult(Kbb, TN));

		int n = kff.r;
		Cmatrix I; I.unit(n);	
		Cmatrix Br(2*n, 2*n), Bl(2*n, 2*n);

		Br.init(CC0);
		Br.setsubmatrix(  I, 0, 0);
		Br.setsubmatrix(kbf, n, 0);
		Br.setsubmatrix(kbb, n, n);
		Bl.init(CC0);
		Bl.setsubmatrix(kff, 0, 0);
		Bl.setsubmatrix(kfb, 0, n);
		Bl.setsubmatrix(  I, n, n);

		Cmatrix B;
		Br.inverse();
		B = mult(Br, Bl);

		Cvector ev;
		int info;
		fprintf(stderr, "EVS: ");
		ev = B.geneigenv(info);

		if(info == 0)
		{
			done = 1;
			fprintf(stderr, "Ok.\n");
			
			ma.clear();
			for(int k=0; k<=ev.nel-1; ++k)
			{
				Complex l;
				l = ev(k);
				if(fabs(l.abs()-1.0) < FBMAXATTLVL)
				{ 
				if(FBVOPLVL > 0.0 || fabs(l.args()/2.0/PI) >  L/lambda*n_min) 
				{

				FBMode m;
				m.st = st;	
				m.pol = pol;
				m.L = L;
				m.k0 = k0;
				m.beta = -l.args()/L;
				m.neff = m.beta/k0;
				m.vM = 0;
				m.ky = CC0;
				m.type = 'B';
				m.fldB = bfld;
				
				if(fabs(m.neff) > 1.0e-9)
				{
				
				Cvector f0, bN; 
				f0 = B.col(k).subvector(n, 0);
				bN = B.col(k).subvector(n, n);

				Cvector tf, tb, sol;
				tf = diagmult(T1, f0);
				tb = diagmult(TN, bN);
				m.fldB.setf(1, SBLEFT, f0);
				m.fldB.setf(1, SBRIGHT, tf);
				sol = add(mult(Kbf, tf), mult(Kbb, tb));
				m.fldB.setb(1, SBRIGHT, sol);
				m.fldB.setb(1, SBLEFT,  diagmult(T1, sol));
				for(int j=2; j<=m.fldB.s.nz-1; ++j)
				{
					sol = add(mult(Iff[j], tf), mult(Ifb[j], tb));
					m.fldB.setf(j, SBLEFT,  sol);
					m.fldB.setf(j, SBRIGHT, diagmult(TX[j], sol));
					sol = add(mult(Ibf[j], tf), mult(Ibb[j], tb));
					m.fldB.setb(j, SBRIGHT, sol);
					m.fldB.setb(j, SBLEFT,  diagmult(TX[j], sol));
				}
				sol = add(mult(Kff, tf), mult(Kfb, tb));
				m.fldB.setf(m.fldB.s.nz, SBLEFT,  sol);
				m.fldB.setf(m.fldB.s.nz, SBRIGHT, diagmult(TN, sol));
				m.fldB.setb(m.fldB.s.nz, SBRIGHT, bN);
				m.fldB.setb(m.fldB.s.nz, SBLEFT, tb);

				m.fldB.setf(0, SBRIGHT, m.fldB.f(1, SBLEFT));	
				m.fldB.setb(0, SBRIGHT, m.fldB.b(1, SBLEFT));	
				m.fldB.setf(m.fldB.s.nz+1, SBLEFT, m.fldB.f(m.fldB.s.nz, SBRIGHT));	
				m.fldB.setb(m.fldB.s.nz+1, SBLEFT, m.fldB.b(m.fldB.s.nz, SBRIGHT));	

				double pnf = 1.0/sqrt(fabs(m.power()));
				for(int j=0; j<=m.fldB.s.nz+1; ++j)
				{
					m.fldB.setf(j, SBLEFT,  mult(m.fldB.f(j, SBLEFT),  pnf));
					m.fldB.setf(j, SBRIGHT, mult(m.fldB.f(j, SBRIGHT), pnf));
					m.fldB.setb(j, SBLEFT,  mult(m.fldB.b(j, SBLEFT),  pnf));
					m.fldB.setb(j, SBRIGHT, mult(m.fldB.b(j, SBRIGHT), pnf));
				}
				
				if(FBVOPLVL > 0.0)
				{
					Interval hcw(st.hz(0), st.hz(st.nz));
					Circuit tst;
					tst = st.circuit();
					Waveguide twg, ll, ul; 
					int tnz = tst.nz;
					twg = Waveguide(tnz, tst.hz, lambda, tst.n.row(0));
					ll = Waveguide(tnz-2, twg.hx.subvector(tnz-1, 1), 
					               lambda, twg.n.subvector(tnz, 1));
					twg = Waveguide(tnz, tst.hz, lambda, tst.n.row(tst.nx+1));
					ul = Waveguide(tnz-2, twg.hx.subvector(tnz-1, 1), 
					               lambda, twg.n.subvector(tnz, 1));
					ModeArray vma, lsd, usd;
					lsd.clear();
					modespectrum(ll, pol, hcw, LIMD, nummx, vma, 1);
					lsd.merge(vma);
					modespectrum(ll, pol, hcw, LIMN, nummx, vma, 1);
					lsd.merge(vma);
					lsd.restrict2prop();

					usd.clear();
					modespectrum(ul, pol, hcw, LIMD, nummx, vma, 1);
					usd.merge(vma);
					modespectrum(ul, pol, hcw, LIMN, nummx, vma, 1);
					usd.merge(vma);
					usd.restrict2prop();
				
					double lp = tst.hx(0)-0.8*(tst.hx(0)-cw.x0);
					double up = tst.hx(tst.nx)+0.8*(cw.x1-tst.hx(tst.nx));

					Cvector lamp(lsd.num);
					lamp.init(CC0);
					Cvector uamp(usd.num);
					uamp.init(CC0);
					
					for(int s=1; s<=st.nz; ++s)
					{
						Interval sec = st.segment(s);
						Cmatrix cff, cfb, cbf, cbb;
						
						lsd.crossoverlap(m.fldB(s), lp, sec, cff, cfb, cbf, cbb);  	
						lamp.addeq(mult(cbf, m.fldB.f(s, SBLEFT)));
						lamp.addeq(mult(cbb, m.fldB.b(s, SBRIGHT)));

						usd.crossoverlap(m.fldB(s), up, sec, cff, cfb, cbf, cbb);  	
						uamp.addeq(mult(cff, m.fldB.f(s, SBLEFT)));
						uamp.addeq(mult(cfb, m.fldB.b(s, SBRIGHT)));
					}
					double nop = m.fldB.Pout(SBRIGHT)+m.fldB.Pout(SBLEFT);
					double pdn = lamp.sqnorm()/nop;
					double pup = uamp.sqnorm()/nop;
					if(pdn+pup < FBVOPLVL) 
					{
						fprintf(stderr, "[%d]  neff = %g  b*L/2/pi = %g  att = %g  vopl = %g  p = %g.\n",
							k, m.neff, l.args()/2.0/PI, l.abs(), pdn+pup, m.power());
						ma.add(m);	
					}
				}
				else
				{
					fprintf(stderr, "[%d]  neff = %g  b*L/2/pi = %g  att = %g  p = %g",
						k, m.neff, l.args()/2.0/PI, l.abs(), m.power());
					if(fabs(B.col(k).norm()-1.0) <1.0e-10)
						fprintf(stderr, " .\n");
					else
						fprintf(stderr, " ?\n");
					ma.add(m);	
				}
				}
				}
				}
			}
			ma.sort();
		}
		else 
		{
			if(red == 1)
			{
				fprintf(stderr, "Failed, Mx: %d ", nummx);
				nummx = (int) ceil(((double) nummx)*0.9);
				fprintf(stderr, "-> %d, retry.\n", nummx);
			}
			else 
			{
				fprintf(stderr, "Failed.\n");
				done = 1;
			}
		}
		delete[] Iff;
		delete[] Ifb;
		delete[] Ibf;
		delete[] Ibb;
		delete[] TX;
	}
	return ma.num;
}

int fbmsolve(SegWgStruct st, Polarization pol, Interval cw, int Mx,
             FBModeArray& ma, int quiet)
{
	return fbmsolve(st, pol, cw, Mx, 1, ma, quiet);
}

/* -- vBEP - implementation --------------------------------------------- */

/* include common lateral wavenumber ky */
int fbmsolve(SegWgStruct st, Polarization pol, 
             Complex ky, 
             Interval cw, int Mx, int red,
             FBModeArray& mam, int quiet)
{
	if(cw.x0 > cw.x1) fbmodeerror("fbmsolve: computational window");
	double lambda = st(0).lambda;
	double k0 = val_k0(lambda);
	if(st.nz < 3) fbmodeerror("fbmsolve: illegal structure (1)");  
	if(st(0).equal(st(1)) != 1) fbmodeerror("fbmsolve: illegal structure (2)");  
	if(st(0).equal(st(st.nz)) != 1) fbmodeerror("fbmsolve: illegal structure (3)");  
	if(st(0).equal(st(st.nz+1)) != 1) fbmodeerror("fbmsolve: illegal structure (4)");  
	if(fabs(ky.im) > 0.0) fbmodeerror("fbmsolve: non-real ky");  

	int done = 0;
	int nummx = Mx;
	if(nummx < FBNUMMXMIN) fbmodeerror("fbmsolve: Mx < FBNUMMXMIN");  
	double e_min = 1.0;
	for(int j=0; j<=st.nz+1; ++j)
	{
		if(st(j).defaultepseffmin() > e_min) 
			e_min = st(j).defaultepseffmin();
	}	
	e_min = e_min-ky.re*ky.re/k0/k0;
	double n_min = 0.0;
	if(e_min > 0.0) n_min = sqrt(e_min);
	double L = st.hz(st.nz)-st.hz(0);
	
	if(quiet != 1)
	{
		Circuit rc = st.circuit();
fprintf(stderr, "\n------------- FBMode, vBEP ----------------------------------------- '11 ---\n");
fprintf(stderr, "VEC  ");
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d, ky: ", cw.x0, cw.x1, Mx); ky.nice(stderr);
fprintf(stderr, "\n");
fprintf(stderr,   "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", rc.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=rc.nx; ++j) fprintf(stderr, "%6.10g  ", rc.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", rc.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=rc.nz; ++j) fprintf(stderr, "%6.10g  ", rc.hz(j));
fprintf(stderr, "\n");
if(rc.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=rc.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=rc.nz+1; ++j) fprintf(stderr, "%6.10g  ", rc.n(i, j));
fprintf(stderr,"\n");
}
fprintf(stderr, "L: %g  ", L);
fprintf(stderr, "lambda: %.10g  k0: %g", rc.lambda, k0);
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
if(FBVOPLVL <= 0.0)
{
	fprintf(stderr, "guided FBModes: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g\n", L*n_min/lambda); 
}
else
{
	fprintf(stderr, "guided FBModes: rel. vert. Output < %g\n", FBVOPLVL); 
	fprintf(stderr, "(light cone: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g)\n", L*n_min/lambda); 
}
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	while(done == 0 && nummx >= FBNUMMXMIN)
	{
		
		BepField bfld(st);
		bfld.discretize(ky, 0, cw, nummx, 1);
		if(bfld.s.nz < 3) fbmodeerror("fbmsolve: BEP initialization"); 
		Cmatrix M0mbf, M0mbb, M0pff, M0pfb;
		Cmatrix MNmbf, MNmbb, MNpff, MNpfb;
		Cmatrix Kff, Kfb, Kbf, Kbb;
		Cmatrix *Iff, *Ifb, *Ibf, *Ibb;
		Cvector *TX;
		Iff = new Cmatrix[bfld.s.nz];
		Ifb = new Cmatrix[bfld.s.nz];
		Ibf = new Cmatrix[bfld.s.nz];
		Ibb = new Cmatrix[bfld.s.nz];        
		TX  = new Cvector[bfld.s.nz];        
		Cvector T1, TN;
		Cmatrix D0, D1, DN, DNp;
		bfld.quepsetup(M0mbf, M0mbb, M0pff, M0pfb, T1, D0, D1,
			    MNmbf, MNmbb, MNpff, MNpfb, TN, DN, DNp,
			    Kff, Kfb, Kbf, Kbb, Iff, Ifb, Ibf, Ibb, TX);     

		Cmatrix kff, kfb, kbf, kbb;
		kff = diagmult(TN, diagmult(Kff, T1));
		kbf = diagmult(T1, diagmult(Kbf, T1));
		kfb = diagmult(TN, diagmult(Kfb, TN));
		kbb = diagmult(T1, diagmult(Kbb, TN));

		int n = kff.r;
		Cmatrix I; I.unit(n);	
		Cmatrix Br(2*n, 2*n), Bl(2*n, 2*n);

		Br.init(CC0);
		Br.setsubmatrix(  I, 0, 0);
		Br.setsubmatrix(kbf, n, 0);
		Br.setsubmatrix(kbb, n, n);
		Bl.init(CC0);
		Bl.setsubmatrix(kff, 0, 0);
		Bl.setsubmatrix(kfb, 0, n);
		Bl.setsubmatrix(  I, n, n);

		Cmatrix B;
		Br.inverse();
		B = mult(Br, Bl);

		Cvector ev;
		int info;
		fprintf(stderr, "EVS: ");
		ev = B.geneigenv(info);

		if(info == 0)
		{
			done = 1;
			fprintf(stderr, "Ok.\n");
			
			mam.clear();
			for(int k=0; k<=ev.nel-1; ++k)
			{
				Complex l;
				l = ev(k);
/*
fprintf(stderr, "(%d) l = ", k); 
l.nice(stderr);
fprintf(stderr, " |l| = %g", l.abs()); 
fprintf(stderr, "\n"); 
*/
				if(fabs(l.abs()-1.0) < FBMAXATTLVL)
				{ 
				if(FBVOPLVL > 0.0 || fabs(l.args()/2.0/PI) >  L/lambda*n_min) 
				{

				FBMode m;
				m.st = st;	
				m.pol = pol;
				m.L = L;
				m.k0 = k0;
				m.beta = -l.args()/L;
				m.neff = m.beta/k0;
				m.vM = 1;
				m.ky = ky;
				m.type = 'B';
				m.fldB = bfld;
				
				if(fabs(m.neff) > 1.0e-9)
				{
				
				Cvector f0, bN; 
				f0 = B.col(k).subvector(n, 0);
				bN = B.col(k).subvector(n, n);

				Cvector tf, tb, sol;
				tf = diagmult(T1, f0);
				tb = diagmult(TN, bN);
				m.fldB.setf(1, SBLEFT, f0);
				m.fldB.setf(1, SBRIGHT, tf);
				sol = add(mult(Kbf, tf), mult(Kbb, tb));
				m.fldB.setb(1, SBRIGHT, sol);
				m.fldB.setb(1, SBLEFT,  diagmult(T1, sol));
				for(int j=2; j<=m.fldB.s.nz-1; ++j)
				{
					sol = add(mult(Iff[j], tf), mult(Ifb[j], tb));
					m.fldB.setf(j, SBLEFT,  sol);
					m.fldB.setf(j, SBRIGHT, diagmult(TX[j], sol));
					sol = add(mult(Ibf[j], tf), mult(Ibb[j], tb));
					m.fldB.setb(j, SBRIGHT, sol);
					m.fldB.setb(j, SBLEFT,  diagmult(TX[j], sol));
				}
				sol = add(mult(Kff, tf), mult(Kfb, tb));
				m.fldB.setf(m.fldB.s.nz, SBLEFT,  sol);
				m.fldB.setf(m.fldB.s.nz, SBRIGHT, diagmult(TN, sol));
				m.fldB.setb(m.fldB.s.nz, SBRIGHT, bN);
				m.fldB.setb(m.fldB.s.nz, SBLEFT, tb);

				m.fldB.setf(0, SBRIGHT, m.fldB.f(1, SBLEFT));	
				m.fldB.setb(0, SBRIGHT, m.fldB.b(1, SBLEFT));	
				m.fldB.setf(m.fldB.s.nz+1, SBLEFT, m.fldB.f(m.fldB.s.nz, SBRIGHT));	
				m.fldB.setb(m.fldB.s.nz+1, SBLEFT, m.fldB.b(m.fldB.s.nz, SBRIGHT));	

				double pnf = 1.0/sqrt(fabs(m.power()));
				for(int j=0; j<=m.fldB.s.nz+1; ++j)
				{
					m.fldB.setf(j, SBLEFT,  mult(m.fldB.f(j, SBLEFT),  pnf));
					m.fldB.setf(j, SBRIGHT, mult(m.fldB.f(j, SBRIGHT), pnf));
					m.fldB.setb(j, SBLEFT,  mult(m.fldB.b(j, SBLEFT),  pnf));
					m.fldB.setb(j, SBRIGHT, mult(m.fldB.b(j, SBRIGHT), pnf));
				}
				
				if(FBVOPLVL > 0.0)
				{
					Interval hcw(st.hz(0), st.hz(st.nz));
					Circuit tst;
					tst = st.circuit();
					Waveguide twg, ll, ul; 
					int tnz = tst.nz;
					twg = Waveguide(tnz, tst.hz, lambda, tst.n.row(0));
					ll = Waveguide(tnz-2, twg.hx.subvector(tnz-1, 1), 
					               lambda, twg.n.subvector(tnz, 1));
					twg = Waveguide(tnz, tst.hz, lambda, tst.n.row(tst.nx+1));
					ul = Waveguide(tnz-2, twg.hx.subvector(tnz-1, 1), 
					               lambda, twg.n.subvector(tnz, 1));
					ModeArray vma, lsd, usd;
					lsd.clear();
					modespectrum(ll, TE, hcw, LIMD, nummx, vma, 1);
					lsd.merge(vma);
					modespectrum(ll, TE, hcw, LIMN, nummx, vma, 1);
					lsd.merge(vma);
					modespectrum(ll, TM, hcw, LIMD, nummx, vma, 1);
					lsd.merge(vma);
					modespectrum(ll, TM, hcw, LIMN, nummx, vma, 1);
					lsd.merge(vma);
					lsd.restrict2prop();
					lsd.vectorify(ky);

					usd.clear();
					modespectrum(ul, TE, hcw, LIMD, nummx, vma, 1);
					usd.merge(vma);
					modespectrum(ul, TE, hcw, LIMN, nummx, vma, 1);
					usd.merge(vma);
					modespectrum(ul, TM, hcw, LIMD, nummx, vma, 1);
					usd.merge(vma);
					modespectrum(ul, TM, hcw, LIMN, nummx, vma, 1);
					usd.merge(vma);
					usd.restrict2prop();
					usd.vectorify(ky);
				
					double lp = tst.hx(0)-0.8*(tst.hx(0)-cw.x0);
					double up = tst.hx(tst.nx)+0.8*(cw.x1-tst.hx(tst.nx));

					Cvector lamp(lsd.num);
					lamp.init(CC0);
					Cvector uamp(usd.num);
					uamp.init(CC0);
					
					for(int s=1; s<=st.nz; ++s)
					{
						Interval sec = st.segment(s);
						Cmatrix cff, cfb, cbf, cbb;
						
						lsd.crossoverlap(m.fldB(s), lp, sec, cff, cfb, cbf, cbb);  	
						lamp.addeq(mult(cbf, m.fldB.f(s, SBLEFT)));
						lamp.addeq(mult(cbb, m.fldB.b(s, SBRIGHT)));

						usd.crossoverlap(m.fldB(s), up, sec, cff, cfb, cbf, cbb);  	
						uamp.addeq(mult(cff, m.fldB.f(s, SBLEFT)));
						uamp.addeq(mult(cfb, m.fldB.b(s, SBRIGHT)));
					}
					double nop = m.fldB.Pout(SBRIGHT)+m.fldB.Pout(SBLEFT);
					double pdn = lamp.sqnorm()/nop;
					double pup = uamp.sqnorm()/nop;
					if(pdn+pup < FBVOPLVL) 
					{
						fprintf(stderr, "(%d)  neff = %g  b*L/2/pi = %g  att = %g  vopl = %g  p = %g.\n",
							k, m.neff, l.args()/2.0/PI, l.abs(), pdn+pup, m.power());
						mam.add(m);	
					}
				}
				else
				{
					fprintf(stderr, "(%d)  neff = %g  b*L/2/pi = %g  att = %g  p = %g",
						k, m.neff, l.args()/2.0/PI, l.abs(), m.power());
					if(fabs(B.col(k).norm()-1.0) <1.0e-10)
						fprintf(stderr, " .\n");
					else
						fprintf(stderr, " ?\n");
					mam.add(m);	
				}
				}
				}
				}
			}
			mam.sort();
		}
		else 
		{
			if(red == 1)
			{
				fprintf(stderr, "Failed, Mx: %d ", nummx);
				nummx = (int) ceil(((double) nummx)*0.9);
				fprintf(stderr, "-> %d, retry.\n", nummx);
			}
			else 
			{
				fprintf(stderr, "Failed.\n");
				done = 1;
			}
		}
		delete[] Iff;
		delete[] Ifb;
		delete[] Ibf;
		delete[] Ibb;
		delete[] TX;
	}
	return mam.num;
}

/* ... with ky given by the incidence of mode of polarization pol, of order oder ord,
   coming in at border bdr, at angle theta (^o) */
int fbmsolve(SegWgStruct st, Polarization pol, 
             SBorder bdr, int ord, double theta,
             Interval cw, int Mx, int red,
             FBModeArray& mam, int quiet)
{
	Waveguide ewg;
	switch(bdr) 
	{
	case SBLEFT:
		ewg = st(0);
		break;
	case SBRIGHT:
		ewg = st(st.nz+1);
		break;
	}
	ModeArray ima;
	modeanalysis(ewg, pol, ima, 1);
	if(ord < 0 || ord >= ima.num) fbmodeerror("fbmsolve, vBEP, setup, mode order");
	double N = ima(ord).neff;
	double b = ima(ord).beta;
	if(theta <= -90.0 || theta >= 90.0) fbmodeerror("fbmsolve, vBEP, setup, theta");
	Complex ky = Complex(b*sin(theta/180.0*PI), 0.0);
	if(quiet != 1)
	{
fprintf(stderr, "\n----------------------------------------------------------------------------\n");
fprintf(stderr, "in: %s, neff = %g, beta = %g; angle = %g^o", 
                ima(ord).ids, N, b, theta);
	}
	return fbmsolve(st, pol, ky, cw, Mx, red, mam, quiet);
}


/* -- vQUEP - implementation --------------------------------------------- */

/* include common lateral wavenumber ky */
int fbmsolveQ(SegWgStruct st, Polarization pol, 
             Complex ky, 
             Interval cwx, int Mx, int Mz,
             FBModeArray& mam, int quiet)
{
	if(cwx.x0 > cwx.x1) fbmodeerror("fbmsolve: computational window");
	double lambda = st(0).lambda;
	double k0 = val_k0(lambda);
	if(st.nz < 3) fbmodeerror("fbmsolve: illegal structure (1)");  
	if(st(0).equal(st(1)) != 1) fbmodeerror("fbmsolve: illegal structure (2)");  
	if(st(0).equal(st(st.nz)) != 1) fbmodeerror("fbmsolve: illegal structure (3)");  
	if(st(0).equal(st(st.nz+1)) != 1) fbmodeerror("fbmsolve: illegal structure (4)");  
	if(fabs(ky.im) > 0.0) fbmodeerror("fbmsolve: non-real ky");  

	double e_min = 1.0;
	for(int j=0; j<=st.nz+1; ++j)
	{
		if(st(j).defaultepseffmin() > e_min) 
			e_min = st(j).defaultepseffmin();
	}	
	e_min = e_min-ky.re*ky.re/k0/k0;
	double n_min = 0.0;
	if(e_min > 0.0) n_min = sqrt(e_min);
	Interval cwz(st.hz(0), st.hz(st.nz));
	double L = cwz.len();
	
	if(quiet != 1)
	{
		Circuit rc = st.circuit();
fprintf(stderr, "\n------------- FBMode, vQUEP ---------------------------------------- '22 ---\n");
fprintf(stderr, "VEC  ");
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d, ", cwx.x0, cwx.x1, Mx); 
fprintf(stderr, "z in (%.10g, %.10g), Mz: %d, ", cwz.x0, cwz.x1, Mz);
fprintf(stderr, "ky/k0: "); (ky/k0).nice(stderr);
fprintf(stderr, "\n");
fprintf(stderr,   "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", rc.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=rc.nx; ++j) fprintf(stderr, "%6.10g  ", rc.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", rc.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=rc.nz; ++j) fprintf(stderr, "%6.10g  ", rc.hz(j));
fprintf(stderr, "\n");
if(rc.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=rc.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=rc.nz+1; ++j) fprintf(stderr, "%6.10g  ", rc.n(i, j));
fprintf(stderr,"\n");
}
fprintf(stderr, "L: %g  ", L);
fprintf(stderr, "lambda: %.10g  k0: %g", rc.lambda, k0);
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
if(FBVOPLVL <= 0.0)
{
	fprintf(stderr, "guided FBModes: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g\n", L*n_min/lambda); 
}
else
{
	fprintf(stderr, "guided FBModes: rel. vert. Output < %g\n", FBVOPLVL); 
	fprintf(stderr, "(light cone: |neff| > %g  ", n_min); 
	fprintf(stderr, "|beta L / (2 pi)| > %g)\n", L*n_min/lambda); 
}
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}

	SegWgStruct qst = st.reduce();
	QuepField qfld(qst, ky, cwx, Mx, cwz, Mz, 1);

	if(qfld.h.s.nz < 3) fbmodeerror("fbmsolve: #hor. segments: h.s.nz < 3");  
	Cmatrix hM0mbf, hM0mbb, hM0pff, hM0pfb;
	Cmatrix hMNmbf, hMNmbb, hMNpff, hMNpfb;
	Cmatrix hKff, hKfb, hKbf, hKbb;
        Cmatrix *hIff, *hIfb, *hIbf, *hIbb;
        Cvector *hTX;
	hIff = new Cmatrix[qfld.h.s.nz];
        hIfb = new Cmatrix[qfld.h.s.nz];
        hIbf = new Cmatrix[qfld.h.s.nz];
        hIbb = new Cmatrix[qfld.h.s.nz];        
        hTX  = new Cvector[qfld.h.s.nz];        
	Cvector hT1, hTN;
	Cmatrix hD0, hD1, hDN, hDNp;
	if(qfld.v.s.nz < 3) fbmodeerror("fbmsolve: #vert. segments: v.s.nz < 3");  
	Cmatrix vM0mbf, vM0mbb, vM0pff, vM0pfb;
	Cmatrix vMNmbf, vMNmbb, vMNpff, vMNpfb;
	Cmatrix vKff, vKfb, vKbf, vKbb;
        Cmatrix *vIff, *vIfb, *vIbf, *vIbb;
        Cvector *vTX;
	vIff = new Cmatrix[qfld.v.s.nz];
        vIfb = new Cmatrix[qfld.v.s.nz];
        vIbf = new Cmatrix[qfld.v.s.nz];
        vIbb = new Cmatrix[qfld.v.s.nz];        
        vTX  = new Cvector[qfld.v.s.nz];        
	Cvector vT1, vTN;
	Cmatrix vD0, vD1, vDN, vDNp;
	qfld.h.quepsetup(hM0mbf, hM0mbb, hM0pff, hM0pfb, hT1, hD0, hD1,
	            hMNmbf, hMNmbb, hMNpff, hMNpfb, hTN, hDN, hDNp,
	            hKff, hKfb, hKbf, hKbb, hIff, hIfb, hIbf, hIbb, hTX);     
	qfld.v.quepsetup(vM0mbf, vM0mbb, vM0pff, vM0pfb, vT1, vD0, vD1,  
	            vMNmbf, vMNmbb, vMNpff, vMNpfb, vTN, vDN, vDNp,
	            vKff, vKfb, vKbf, vKbb, vIff, vIfb, vIbf, vIbb, vTX);     
	fprintf(stderr, "QUEP (");
	fprintf(stderr, "%d x", qfld.h.s.nz+2);
	int rfmin, rfmax, rfn;
	rfmin = rfmax = qfld.h(0).num;
	for(int i=1; i<=qfld.h.s.nz+1; ++i)
	{
		rfn = qfld.h(i).num;
		if(rfmin > rfn) rfmin = rfn;
		if(rfmax < rfn) rfmax = rfn;
	}
	if(rfmin == rfmax) fprintf(stderr, " %d)*(", rfmin);
	else fprintf(stderr, " (%d-%d))*(", rfmin, rfmax);
	fprintf(stderr, "%d x", qfld.v.s.nz+2);
	rfmin = rfmax = qfld.v(0).num;
	for(int i=1; i<=qfld.v.s.nz+1; ++i)
	{
		rfn = qfld.v(i).num;
		if(rfmin > rfn) rfmin = rfn;
		if(rfmax < rfn) rfmax = rfn;
	}
	if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin);
	else fprintf(stderr, " (%d-%d)) ", rfmin, rfmax);

	Cmatrix vmat;
	Cmatrix rmat;

	int dh, dvf, dvb;
	dvf = qfld.v(1).num;
	dvb = qfld.v(qfld.v.s.nz).num;

	Cmatrix hR0mff, hR0mfb, hR0mbf, hR0mbb; 
	dh  = qfld.h(0).num;
	vmat = qfld.v.crossoverlap(qfld.h(0), qfld.h.s.hz(0), 
               vT1, vTN, vKff, vKfb, vKbf, vKbb, vIff, vIfb, vIbf, vIbb); 
	rmat = mult(hD0,  vmat);  
	hR0mff = rmat.submatrix(dh, dvf, 0,  0);
	hR0mfb = rmat.submatrix(dh, dvb, 0,  dvf);
	hR0mbf = rmat.submatrix(dh, dvf, dh, 0);
	hR0mbb = rmat.submatrix(dh, dvb, dh, dvf);

	fprintf(stderr, ".");
	
	Cmatrix hR0pff, hR0pfb, hR0pbf, hR0pbb; 
	if(qfld.h(1).equal(qfld.h(0)) != 1)
	{
		dh  = qfld.h(1).num;
		vmat = qfld.v.crossoverlap(qfld.h(1), qfld.h.s.hz(0), 
		     vT1, vTN, vKff, vKfb, vKbf, vKbb, vIff, vIfb, vIbf, vIbb);
	}
	rmat = mult(hD1,  vmat); 
	hR0pff = rmat.submatrix(dh, dvf, 0,  0);
	hR0pfb = rmat.submatrix(dh, dvb, 0,  dvf);
	hR0pbf = rmat.submatrix(dh, dvf, dh, 0);
	hR0pbb = rmat.submatrix(dh, dvb, dh, dvf);

	fprintf(stderr, ".");

	Cmatrix hRNmff, hRNmfb, hRNmbf, hRNmbb; 
	dh  = qfld.h(qfld.h.s.nz).num;
	vmat = qfld.v.crossoverlap(qfld.h(qfld.h.s.nz), qfld.h.s.hz(qfld.h.s.nz), 
               vT1, vTN, vKff, vKfb, vKbf, vKbb, vIff, vIfb, vIbf, vIbb); 
	rmat = mult(hDN,  vmat); 
	hRNmff = rmat.submatrix(dh, dvf, 0,  0);
	hRNmfb = rmat.submatrix(dh, dvb, 0,  dvf);
	hRNmbf = rmat.submatrix(dh, dvf, dh, 0);
	hRNmbb = rmat.submatrix(dh, dvb, dh, dvf);

	fprintf(stderr, ".");

	Cmatrix hRNpff, hRNpfb, hRNpbf, hRNpbb; 
	if(qfld.h(qfld.h.s.nz+1).equal(qfld.h(qfld.h.s.nz)) != 1)
	{
		dh  = qfld.h(qfld.h.s.nz+1).num;
		vmat = qfld.v.crossoverlap(qfld.h(qfld.h.s.nz+1), qfld.h.s.hz(qfld.h.s.nz), 
		     vT1, vTN, vKff, vKfb, vKbf, vKbb, vIff, vIfb, vIbf, vIbb);
	}
	rmat = mult(hDNp, vmat); 
	hRNpff = rmat.submatrix(dh, dvf, 0,  0);
	hRNpfb = rmat.submatrix(dh, dvb, 0,  dvf);
	hRNpbf = rmat.submatrix(dh, dvf, dh, 0);
	hRNpbb = rmat.submatrix(dh, dvb, dh, dvf);

	fprintf(stderr, ".");

	int dv, dhf, dhb;
	dhf = qfld.h(1).num;
	dhb = qfld.h(qfld.h.s.nz).num;

	Cmatrix vR0mff, vR0mfb, vR0mbf, vR0mbb; 
	dv  = qfld.v(0).num;
	vmat = qfld.h.crossoverlap(qfld.v(0), qfld.v.s.hz(0), 
               hT1, hTN, hKff, hKfb, hKbf, hKbb, hIff, hIfb, hIbf, hIbb); 
	rmat = mult(vD0,  vmat); 
	vR0mff = rmat.submatrix(dv, dhf, 0,  0);
	vR0mfb = rmat.submatrix(dv, dhb, 0,  dhf);
	vR0mbf = rmat.submatrix(dv, dhf, dv, 0);
	vR0mbb = rmat.submatrix(dv, dhb, dv, dhf);

	fprintf(stderr, ".");

	Cmatrix vR0pff, vR0pfb, vR0pbf, vR0pbb; 
	if(qfld.v(1).equal(qfld.v(0)) != 1)
	{
		dv  = qfld.v(1).num;
		vmat = qfld.h.crossoverlap(qfld.v(1), qfld.v.s.hz(0), 
		     hT1, hTN, hKff, hKfb, hKbf, hKbb, hIff, hIfb, hIbf, hIbb);
	}
	rmat = mult(vD1,  vmat); 
	vR0pff = rmat.submatrix(dv, dhf, 0,  0);
	vR0pfb = rmat.submatrix(dv, dhb, 0,  dhf);
	vR0pbf = rmat.submatrix(dv, dhf, dv, 0);
	vR0pbb = rmat.submatrix(dv, dhb, dv, dhf);

	fprintf(stderr, ".");

	Cmatrix vRNmff, vRNmfb, vRNmbf, vRNmbb; 
	dv  = qfld.v(qfld.v.s.nz).num;
	vmat = qfld.h.crossoverlap(qfld.v(qfld.v.s.nz), qfld.v.s.hz(qfld.v.s.nz), 
               hT1, hTN, hKff, hKfb, hKbf, hKbb, hIff, hIfb, hIbf, hIbb); 
	rmat = mult(vDN,  vmat); 
	vRNmff = rmat.submatrix(dv, dhf, 0,  0);
	vRNmfb = rmat.submatrix(dv, dhb, 0,  dhf);
	vRNmbf = rmat.submatrix(dv, dhf, dv, 0);
	vRNmbb = rmat.submatrix(dv, dhb, dv, dhf);

	fprintf(stderr, ".");

	Cmatrix vRNpff, vRNpfb, vRNpbf, vRNpbb; 
	if(qfld.v(qfld.v.s.nz+1).equal(qfld.v(qfld.v.s.nz)) != 1)
	{
		dv  = qfld.v(qfld.v.s.nz+1).num;
		vmat = qfld.h.crossoverlap(qfld.v(qfld.v.s.nz+1), qfld.v.s.hz(qfld.v.s.nz), 
		     hT1, hTN, hKff, hKfb, hKbf, hKbb, hIff, hIfb, hIbf, hIbb);
	}
	rmat = mult(vDNp, vmat); 
	vRNpff = rmat.submatrix(dv, dhf, 0,  0);
	vRNpfb = rmat.submatrix(dv, dhb, 0,  dhf);
	vRNpbf = rmat.submatrix(dv, dhf, dv, 0);
	vRNpbb = rmat.submatrix(dv, dhb, dv, dhf);

	vmat.strip();
	rmat.strip();
	fprintf(stderr, ".:");

	Cmatrix u, qs, a; 
	int df1, dbn, du1, ddn; 
	df1 = qfld.h(1).num;
	dbn = qfld.h(qfld.h.s.nz).num;
	du1 = qfld.v(1).num;
	ddn = qfld.v(qfld.v.s.nz).num;
	Cmatrix q(df1+dbn+du1+ddn, df1+dbn+du1+ddn);
	
	u.unit(df1);  
	a = add(hM0mbf, mult(diagmult(hM0mbb, hT1), diagmult(hKbf, hT1)));
	qs = sub(u, mult(hM0pfb, a)); 
	q.setsubmatrix(qs, 0, 0);
	a = mult(hM0pfb, mult(diagmult(hM0mbb, hT1), diagmult(hKbb, hTN)));
	qs = mult(a, -1.0);
	q.setsubmatrix(qs, 0, df1);
	qs = sub(hR0pff, mult(hM0pfb, hR0mbf));
	q.setsubmatrix(qs, 0, df1+dbn);
	qs = sub(hR0pfb, mult(hM0pfb, hR0mbb));
	q.setsubmatrix(qs, 0, df1+dbn+du1);

	fprintf(stderr, ".");

	u.unit(dbn);  
	a = mult(hMNmbf, mult(diagmult(hMNpff, hTN), diagmult(hKff, hT1)));
	qs = mult(a, -1.0);
	q.setsubmatrix(qs, df1, 0);
	a = add(hMNpfb, mult(diagmult(hMNpff, hTN), diagmult(hKfb, hTN)));
	qs = sub(u, mult(hMNmbf, a)); 
	q.setsubmatrix(qs, df1, df1);
	qs = sub(hRNmbf, mult(hMNmbf, hRNpff));
	q.setsubmatrix(qs, df1, df1+dbn);
	qs = sub(hRNmbb, mult(hMNmbf, hRNpfb));
	q.setsubmatrix(qs, df1, df1+dbn+du1);

	fprintf(stderr, ".");

	u.unit(du1);  
	qs = sub(vR0pff, mult(vM0pfb, vR0mbf));
	q.setsubmatrix(qs, df1+dbn, 0);
	qs = sub(vR0pfb, mult(vM0pfb, vR0mbb));
	q.setsubmatrix(qs, df1+dbn, df1);
	a = add(vM0mbf, mult(diagmult(vM0mbb, vT1), diagmult(vKbf, vT1)));
	qs = sub(u, mult(vM0pfb, a)); 
	q.setsubmatrix(qs, df1+dbn, df1+dbn);
	a = mult(vM0pfb, mult(diagmult(vM0mbb, vT1), diagmult(vKbb, vTN)));
	qs = mult(a, -1.0);
	q.setsubmatrix(qs, df1+dbn, df1+dbn+du1);

	fprintf(stderr, ".");

	u.unit(ddn);  
	qs = sub(vRNmbf, mult(vMNmbf, vRNpff));
	q.setsubmatrix(qs, df1+dbn+du1, 0);
	qs = sub(vRNmbb, mult(vMNmbf, vRNpfb));
	q.setsubmatrix(qs, df1+dbn+du1, df1);
	a = mult(vMNmbf, mult(diagmult(vMNpff, vTN), diagmult(vKff, vT1)));
	qs = mult(a, -1.0);
	q.setsubmatrix(qs, df1+dbn+du1, df1+dbn);
	a = add(vMNpfb, mult(diagmult(vMNpff, vTN), diagmult(vKfb, vTN)));
	qs = sub(u, mult(vMNmbf, a)); 
	q.setsubmatrix(qs, df1+dbn+du1, df1+dbn+du1);

	u.strip();
	a.strip();
	qs.strip();
	fprintf(stderr, ".:");

	q.inverse();
	fprintf(stderr, ".:");

	Cmatrix M(df1+dbn+du1+ddn, df1+dbn+du1+ddn);
	M.init(CC0);
	M.setsubmatrix(hM0pff, 0, 0);
	M.setsubmatrix(hMNmbb, df1, df1);
	M.setsubmatrix(vM0pff, df1+dbn, df1+dbn);
	M.setsubmatrix(vMNmbb, df1+dbn+du1, df1+dbn+du1);

	Cmatrix W(df1+dbn+du1+ddn, df1+dbn+du1+ddn);
	W.init(CC0);

	W.setsubmatrix(mult(diagmult(hMNpff, hTN), diagmult(hKff, hT1)),              0, 0);
	W.setsubmatrix(add(hMNpfb, mult(diagmult(hMNpff, hTN), diagmult(hKfb, hTN))), 0, df1);
	W.setsubmatrix(hRNpff,                                                        0, df1+dbn);
	W.setsubmatrix(hRNpfb,                                                        0, df1+dbn+du1);

	W.setsubmatrix(add(hM0mbf, mult(diagmult(hM0mbb, hT1), diagmult(hKbf, hT1))), df1, 0);
	W.setsubmatrix(mult(diagmult(hM0mbb, hT1), diagmult(hKbb, hTN)),              df1, df1);
	W.setsubmatrix(hR0mbf,                                                        df1, df1+dbn);
	W.setsubmatrix(hR0mbb,                                                        df1, df1+dbn+du1);

	W.setsubmatrix(vRNpff,                                                        df1+dbn, 0);
	W.setsubmatrix(vRNpfb,                                                        df1+dbn, df1);
	W.setsubmatrix(mult(diagmult(vMNpff, vTN), diagmult(vKff, vT1)),              df1+dbn, df1+dbn);
	W.setsubmatrix(add(vMNpfb, mult(diagmult(vMNpff, vTN), diagmult(vKfb, vTN))), df1+dbn, df1+dbn+du1);

	W.setsubmatrix(vR0mbf,                                                        df1+dbn+du1, 0);
	W.setsubmatrix(vR0mbb,                                                        df1+dbn+du1, df1);
	W.setsubmatrix(add(vM0mbf, mult(diagmult(vM0mbb, vT1), diagmult(vKbf, vT1))), df1+dbn+du1, df1+dbn);
	W.setsubmatrix(mult(diagmult(vM0mbb, vT1), diagmult(vKbb, vTN)),              df1+dbn+du1, df1+dbn+du1);

	Cmatrix A;
	A = mult(W, mult(q, M));

	int n = df1;
	if(dbn != n) fbmodeerror("fbmsolve: mismatching discretizations");  

	Cmatrix Aff = A.submatrix(n, n, 0, 0);
	Cmatrix Afb = A.submatrix(n, n, 0, n);
	Cmatrix Abf = A.submatrix(n, n, n, 0);
	Cmatrix Abb = A.submatrix(n, n, n, n);

	Cmatrix I; I.unit(n);	
	Cmatrix Br(2*n, 2*n), Bl(2*n, 2*n);

	Br.init(CC0);
	Br.setsubmatrix(  I, 0, 0);
	Br.setsubmatrix(Abf, n, 0);
	Br.setsubmatrix(Abb, n, n);
	Bl.init(CC0);
	Bl.setsubmatrix(Aff, 0, 0);
	Bl.setsubmatrix(Afb, 0, n);
	Bl.setsubmatrix(  I, n, n);

	Cmatrix B;
	Br.inverse();
	B = mult(Br, Bl);

	fprintf(stderr, ".\n");

	Cvector ev;
	int info;
	fprintf(stderr, "EVS: ");
	ev = B.geneigenv(info);

/* - - - - - */

	if(info == 0)
	{
		fprintf(stderr, "Ok.\n");
		
		mam.clear();
		for(int k=0; k<=ev.nel-1; ++k)
		{
			Complex l;
			l = ev(k);
/*
fprintf(stderr, "(%d) l = ", k); 
l.nice(stderr);
fprintf(stderr, " |l| = %g", l.abs()); 
fprintf(stderr, "\n"); 
*/
			if(fabs(l.abs()-1.0) < FBMAXATTLVL)
			{ 
			if(FBVOPLVL > 0.0 || fabs(l.args()/2.0/PI) >  L/lambda*n_min) 
			{

			FBMode m;
			m.st = st;	
			m.pol = pol;
			m.L = L;
			m.k0 = k0;
			m.beta = -l.args()/L;
			m.neff = m.beta/k0;
			m.vM = 1;
			m.ky = ky;
			m.type = 'Q';
			m.fldQ = qfld;
			
			if(fabs(m.neff) > 1.0e-9)
			{

			Cvector f0, bN, u0, dN; 
			f0 = B.col(k).subvector(n, 0);
			bN = B.col(k).subvector(n, n);
			u0 = Cvector(du1); u0.init(CC0);
			dN = Cvector(ddn); dN.init(CC0);

			m.fldQ.h.setf(0, SBRIGHT, f0);
			m.fldQ.h.setf(0, SBLEFT, f0);
			m.fldQ.h.setb(m.fldQ.h.s.nz+1, SBLEFT, bN);
			m.fldQ.h.setb(m.fldQ.h.s.nz+1, SBRIGHT, bN);
			m.fldQ.v.setf(0, SBRIGHT, u0);
			m.fldQ.v.setf(0, SBLEFT, u0);
			m.fldQ.v.setb(m.fldQ.v.s.nz+1, SBLEFT, dN);
			m.fldQ.v.setb(m.fldQ.v.s.nz+1, SBRIGHT, dN);
/* - - - - - */
			Cvector rs;
			Cvector rhs(df1+dbn+du1+ddn);

			rs = mult(hM0pff, m.fldQ.h.f(0, SBRIGHT));
			rhs.setsubvector(rs, 0);
			rs = mult(hMNmbb, m.fldQ.h.b(m.fldQ.h.s.nz+1, SBLEFT));
			rhs.setsubvector(rs, df1);
			rs = mult(vM0pff, m.fldQ.v.f(0, SBRIGHT));
			rhs.setsubvector(rs, df1+dbn);
			rs = mult(vMNmbb, m.fldQ.v.b(m.fldQ.v.s.nz+1, SBLEFT));
			rhs.setsubvector(rs, df1+dbn+du1);

			Cvector sol;
			sol = mult(q, rhs);

			rhs.strip();

			Cvector f1, bn, u1, dn;
			f1 = sol.subvector(df1, 0);
			bn = sol.subvector(dbn, df1);
			u1 = sol.subvector(du1, df1+dbn);
			dn = sol.subvector(ddn, df1+dbn+du1);

			sol = mult(add(hMNpfb, mult(diagmult(hMNpff, hTN), 
						    diagmult(hKfb, hTN))), bn);
			sol.addeq(mult(mult(diagmult(hMNpff, hTN), diagmult(hKff, hT1)), f1));
			sol.addeq(mult(hRNpff, u1));
			sol.addeq(mult(hRNpfb, dn));
			m.fldQ.h.setf(m.fldQ.h.s.nz+1, SBLEFT, sol);

			sol = mult(add(hM0mbf, mult(diagmult(hM0mbb, hT1), 
						    diagmult(hKbf, hT1))), f1);
			sol.addeq(mult(mult(diagmult(hM0mbb, hT1), diagmult(hKbb, hTN)), bn));
			sol.addeq(mult(hR0mbf, u1));
			sol.addeq(mult(hR0mbb, dn));
			m.fldQ.h.setb(0, SBRIGHT, sol);

			sol = mult(add(vMNpfb, mult(diagmult(vMNpff, vTN), 
						    diagmult(vKfb, vTN))), dn);
			sol.addeq(mult(mult(diagmult(vMNpff, vTN), diagmult(vKff, vT1)), u1));
			sol.addeq(mult(vRNpff, f1));
			sol.addeq(mult(vRNpfb, bn));
			m.fldQ.v.setf(m.fldQ.v.s.nz+1, SBLEFT, sol);

			sol = mult(add(vM0mbf, mult(diagmult(vM0mbb, vT1), 
						    diagmult(vKbf, vT1))), u1);
			sol.addeq(mult(mult(diagmult(vM0mbb, vT1), diagmult(vKbb, vTN)), dn));
			sol.addeq(mult(vR0mbf, f1));
			sol.addeq(mult(vR0mbb, bn));
			m.fldQ.v.setb(0, SBRIGHT, sol);

			Cvector tf, tb;
			tf = diagmult(hT1, f1);
			tb = diagmult(hTN, bn);
			m.fldQ.h.setf(1, SBLEFT, f1);
			m.fldQ.h.setf(1, SBRIGHT, tf);
			sol = add(mult(hKbf, tf), mult(hKbb, tb));
			m.fldQ.h.setb(1, SBRIGHT, sol);
			m.fldQ.h.setb(1, SBLEFT,  diagmult(hT1, sol));
			for(int j=2; j<=m.fldQ.h.s.nz-1; ++j)
			{
				sol = add(mult(hIff[j], tf), mult(hIfb[j], tb));
				m.fldQ.h.setf(j, SBLEFT,  sol);
				m.fldQ.h.setf(j, SBRIGHT, diagmult(hTX[j], sol));
				sol = add(mult(hIbf[j], tf), mult(hIbb[j], tb));
				m.fldQ.h.setb(j, SBRIGHT, sol);
				m.fldQ.h.setb(j, SBLEFT,  diagmult(hTX[j], sol));
			}
			sol = add(mult(hKff, tf), mult(hKfb, tb));
			m.fldQ.h.setf(m.fldQ.h.s.nz, SBLEFT,  sol);
			m.fldQ.h.setf(m.fldQ.h.s.nz, SBRIGHT, diagmult(hTN, sol));
			m.fldQ.h.setb(m.fldQ.h.s.nz, SBRIGHT, bn);
			m.fldQ.h.setb(m.fldQ.h.s.nz, SBLEFT, tb);

			tf = diagmult(vT1, u1);
			tb = diagmult(vTN, dn);
			m.fldQ.v.setf(1, SBLEFT, u1);
			m.fldQ.v.setf(1, SBRIGHT, tf);
			sol = add(mult(vKbf, tf), mult(vKbb, tb));
			m.fldQ.v.setb(1, SBRIGHT, sol);
			m.fldQ.v.setb(1, SBLEFT,  diagmult(vT1, sol));
			for(int j=2; j<=m.fldQ.v.s.nz-1; ++j)
			{
				sol = add(mult(vIff[j], tf), mult(vIfb[j], tb));
				m.fldQ.v.setf(j, SBLEFT,  sol);
				m.fldQ.v.setf(j, SBRIGHT, diagmult(vTX[j], sol));
				sol = add(mult(vIbf[j], tf), mult(vIbb[j], tb));
				m.fldQ.v.setb(j, SBRIGHT, sol);
				m.fldQ.v.setb(j, SBLEFT,  diagmult(vTX[j], sol));
			}
			sol = add(mult(vKff, tf), mult(vKfb, tb));
			m.fldQ.v.setf(m.fldQ.v.s.nz, SBLEFT,  sol);
			m.fldQ.v.setf(m.fldQ.v.s.nz, SBRIGHT, diagmult(vTN, sol));
			m.fldQ.v.setb(m.fldQ.v.s.nz, SBRIGHT, dn);
			m.fldQ.v.setb(m.fldQ.v.s.nz, SBLEFT, tb);

/* - - - - - */
			double pnf = 1.0/sqrt(fabs(m.power()));
			for(int j=0; j<=m.fldQ.h.s.nz+1; ++j)
			{
				m.fldQ.h.setf(j, SBLEFT,  mult(m.fldQ.h.f(j, SBLEFT),  pnf));
				m.fldQ.h.setf(j, SBRIGHT, mult(m.fldQ.h.f(j, SBRIGHT), pnf));
				m.fldQ.h.setb(j, SBLEFT,  mult(m.fldQ.h.b(j, SBLEFT),  pnf));
				m.fldQ.h.setb(j, SBRIGHT, mult(m.fldQ.h.b(j, SBRIGHT), pnf));
			}
			for(int j=0; j<=m.fldQ.v.s.nz+1; ++j)
			{
				m.fldQ.v.setf(j, SBLEFT,  mult(m.fldQ.v.f(j, SBLEFT),  pnf));
				m.fldQ.v.setf(j, SBRIGHT, mult(m.fldQ.v.f(j, SBRIGHT), pnf));
				m.fldQ.v.setb(j, SBLEFT,  mult(m.fldQ.v.b(j, SBLEFT),  pnf));
				m.fldQ.v.setb(j, SBRIGHT, mult(m.fldQ.v.b(j, SBRIGHT), pnf));
			}

			if(FBVOPLVL > 0.0)
			{
				double nop = m.fldQ.Pout(RIGHT)+m.fldQ.Pout(LEFT);
				double pdn = m.fldQ.Pout(BOTTOM)/nop;
				double pup = m.fldQ.Pout(TOP)/nop;
				if(pdn+pup < FBVOPLVL) 
				{
					fprintf(stderr, "(%d)  neff = %g  b*L/2/pi = %g  att = %g  vopl = %g  p = %g.\n",
						k, m.neff, l.args()/2.0/PI, l.abs(), pdn+pup, m.power());
					mam.add(m);	
				}
			}
			else
			{
				fprintf(stderr, "(%d)  neff = %g  b*L/2/pi = %g  att = %g  p = %g",
					k, m.neff, l.args()/2.0/PI, l.abs(), m.power());
				if(fabs(B.col(k).norm()-1.0) <1.0e-10)
					fprintf(stderr, " .\n");
				else
					fprintf(stderr, " ?\n");
				mam.add(m);	
			}
			}
			}
			}
		}
		mam.sort();
	}

	delete[] hIff;
        delete[] hIfb;
        delete[] hIbf;
        delete[] hIbb;
        delete[] hTX;

	delete[] vIff;
        delete[] vIfb;
        delete[] vIbf;
        delete[] vIbb;
        delete[] vTX;

	return mam.num;
}



/* ... with ky given by the incidence of mode of polarization pol, of order oder ord,
   coming in at border bdr, at angle theta (^o) */
int fbmsolveQ(SegWgStruct st, Polarization pol, 
             SBorder bdr, int ord, double theta,
             Interval cwx, int Mx, int Mz,
             FBModeArray& mam, int quiet)
{
	Waveguide ewg;
	switch(bdr) 
	{
	case SBLEFT:
		ewg = st(0);
		break;
	case SBRIGHT:
		ewg = st(st.nz+1);
		break;
	}
	ModeArray ima;
	modeanalysis(ewg, pol, ima, 1);
	if(ord < 0 || ord >= ima.num) fbmodeerror("fbmsolve, vQUEP, setup, mode order");
	double N = ima(ord).neff;
	double b = ima(ord).beta;
	if(theta <= -90.0 || theta >= 90.0) fbmodeerror("fbmsolve, vQUEP, setup, theta");
	Complex ky = Complex(b*sin(theta/180.0*PI), 0.0);
	if(quiet != 1)
	{
fprintf(stderr, "\n----------------------------------------------------------------------------\n");
fprintf(stderr, "in: %s, neff = %g, beta = %g; angle = %g^o", 
                ima(ord).ids, N, b, theta);
	}
	return fbmsolveQ(st, pol, ky, cwx, Mx, Mz, mam, quiet);
}