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

/*
 * bepfld.cpp
 * Bidirectional eigenmode expansion, sequences of arrays of modes
 */

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

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

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

/* initialize */
BepField::BepField()
{
	s = SegWgStruct(-1);
	marr = new ModeArray[s.nz+2];
	afl  = new Cvector[s.nz+2];
	afr  = new Cvector[s.nz+2];
	abl  = new Cvector[s.nz+2];
	abr  = new Cvector[s.nz+2];
	nums = s.nz+2;
	for(int i=0; i<=s.nz+1; ++i) marr[i].clear();
	BF_vM = 0;
}
BepField::BepField(SegWgStruct st)
{
	s = st;
	marr = new ModeArray[s.nz+2];
	afl  = new Cvector[s.nz+2];
	afr  = new Cvector[s.nz+2];
	abl  = new Cvector[s.nz+2];
	abr  = new Cvector[s.nz+2];
	nums = s.nz+2;
	for(int i=0; i<=s.nz+1; ++i) marr[i].clear();
	BF_vM = 0;
}
BepField::BepField(Circuit rc)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	SegWgStruct st(rc);
	(*this) = BepField(st);
}
BepField::BepField(Circuit rc, int opt)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	SegWgStruct st(rc, opt);
	(*this) = BepField(st);
}

/* full BEP initialization: setup the container for the quadridirectional 
   mode expansion, compute the spectral discretization 
   rc: the structure under investigation, 
   p:  fields of polarization p,
   wx: vertical computational window,
   Mx: number of expansion terms per vertical slice,
   suppress log output, if quiet == 1
   optimize structure, if opt == 1 */
BepField::BepField(Circuit rc, int opt, Polarization p, Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	if(quiet != 1)
	{
fprintf(stderr, "\n------------- BEP -------------------------------------------------- '04 ---\n");
fprintf(stderr, "T%c  ", polCHR(p));
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d\n", wx.x0, wx.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, "lambda: %.10g  k0: %g", rc.lambda, val_k0(rc.lambda));
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	(*this) = BepField(rc, opt);
	discretize(p, wx, Mx, quiet);
}
BepField::BepField(Circuit rc, Polarization p, Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc, 1, p, wx, Mx, 0);
}
BepField::BepField(Circuit rc, Polarization p, Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc, p, wx, Mx, 0);
}
BepField::BepField(SegWgStruct rc, Polarization p, 
                   Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), p, wx, Mx, quiet);
}
BepField::BepField(SegWgStruct rc, Polarization p, 
                   Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), p, wx, Mx, 0);
}
BepField::BepField(SegWgStruct rc, int opt, Polarization p, 
                   Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), opt, p, wx, Mx, quiet);
}


/* destroy */
BepField::~BepField()
{
	if(marr != NULL) delete[] marr;
	if(afl  != NULL) delete[] afl;
	if(afr  != NULL) delete[] afr;
	if(abl  != NULL) delete[] abl;
	if(abr  != NULL) delete[] abr;
	marr = NULL;
	afl  = NULL;
	afr  = NULL;
	abl  = NULL;
	abr  = NULL;
	nums = 0;
	s.clear();
}

/* copy constructor */
BepField::BepField(const BepField& mef)
{
	s = mef.s; 
	marr = new ModeArray[s.nz+2];
	nums = s.nz+2;
	for(int i=0; i<=s.nz+1; ++i) marr[i] = mef.marr[i]; 
	afl = new Cvector[s.nz+2];
	afr = new Cvector[s.nz+2];
	abl = new Cvector[s.nz+2];
	abr = new Cvector[s.nz+2];
	for(int i=0; i<=s.nz+1; ++i) afl[i] = mef.afl[i]; 
	for(int i=0; i<=s.nz+1; ++i) afr[i] = mef.afr[i]; 
	for(int i=0; i<=s.nz+1; ++i) abl[i] = mef.abl[i]; 
	for(int i=0; i<=s.nz+1; ++i) abr[i] = mef.abr[i]; 
	BF_pol    = mef.BF_pol;
	BF_lambda = mef.BF_lambda;
	BF_vM     = mef.BF_vM;
	BF_ky     = mef.BF_ky;
}
	
/* assignment */
BepField& BepField::operator=(const BepField& mef)
{
	if(this != &mef)
	{
		if(marr != NULL) delete[] marr;
		if(afl  != NULL) delete[] afl;
		if(afr  != NULL) delete[] afr;
		if(abl  != NULL) delete[] abl;
		if(abr  != NULL) delete[] abr;
		marr = NULL;
		afl  = NULL;
		afr  = NULL;
		abl  = NULL;
		abr  = NULL;
		nums = 0;
		s = mef.s; 
		marr = new ModeArray[s.nz+2];
		nums = s.nz+2;
		for(int i=0; i<=s.nz+1; ++i) marr[i] = mef.marr[i]; 
		afl = new Cvector[s.nz+2];
		afr = new Cvector[s.nz+2];
		abl = new Cvector[s.nz+2];
		abr = new Cvector[s.nz+2];
		for(int i=0; i<=s.nz+1; ++i) afl[i] = mef.afl[i]; 
		for(int i=0; i<=s.nz+1; ++i) afr[i] = mef.afr[i]; 
		for(int i=0; i<=s.nz+1; ++i) abl[i] = mef.abl[i]; 
		for(int i=0; i<=s.nz+1; ++i) abr[i] = mef.abr[i]; 
		BF_pol    = mef.BF_pol;
                BF_lambda = mef.BF_lambda;
                BF_vM     = mef.BF_vM;
                BF_ky     = mef.BF_ky;
	}
	return *this;
}

/* delete all segment entries */
void BepField::clear()
{
	s.clear();
	if(marr != NULL) delete[] marr; 
	if(afl  != NULL) delete[] afl;
	if(afr  != NULL) delete[] afr;
	if(abl  != NULL) delete[] abl;
	if(abr  != NULL) delete[] abr;
	marr = NULL;
	afl  = NULL;
	afr  = NULL;
	abl  = NULL;
	abr  = NULL;
	nums = 0;
}

/* set mode fields on all segments to ma */
void BepField::init(ModeArray ma)
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("init: SegWgStruct cleared" ); 
	if(nums < s.nz+2)  bepfielderror("init: structure mismatch"); 
	for(int i=0; i<=s.nz+1; ++i) marr[i] = ma; 
}

/* subscripting 
ModeArray& BepField::operator() (int seg)
{ 	
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("access: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("access ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("access ( > s.nz+1)");
	if(seg >= nums) bepfielderror("access: struct mismatch"); 
	return marr[seg]; 
}
ModeArray BepField::operator() (int seg) const 
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("access: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("access ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("access ( > s.nz+1)");
	if(seg >= nums) bepfielderror("access: struct mismatch"); 
	return marr[seg]; 
}
*/

/* input from FILE dat */
void BepField::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		bepfielderror("BepField: read");
	}
	else
	{
		if(marr != NULL) delete[] marr;
		if(afl  != NULL) delete[] afl;
		if(afr  != NULL) delete[] afr;
		if(abl  != NULL) delete[] abl;
		if(abr  != NULL) delete[] abr;
		marr = NULL;
		afl  = NULL;
		afr  = NULL;
		abl  = NULL;
		abr  = NULL;
		nums = 0;
		s.read(dat);
		if(s.nz<-1) bepfielderror("read: s.nz < -1");
		marr = new ModeArray[s.nz+2];
		nums = s.nz+2;
		for(int i=0; i<=s.nz+1; ++i) (marr[i]).read(dat);
		afl = new Cvector[s.nz+2];
		afr = new Cvector[s.nz+2];
		abl = new Cvector[s.nz+2];
		abr = new Cvector[s.nz+2];
		for(int i=0; i<=s.nz+1; ++i) (afl[i]).read(dat);
		for(int i=0; i<=s.nz+1; ++i) (afr[i]).read(dat);
		for(int i=0; i<=s.nz+1; ++i) (abl[i]).read(dat);
		for(int i=0; i<=s.nz+1; ++i) (abr[i]).read(dat);
	}
	check();
	return;
}

/* output to FILE dat */
void BepField::write(FILE *dat) const
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("write: SegWgStruct cleared" ); 
	comment("BepField", dat);
	s.write(dat);
	for(int i=0; i<=s.nz+1; ++i) (marr[i]).write(dat);
	comment("fl's", dat);
	for(int i=0; i<=s.nz+1; ++i) (afl[i]).write(dat);
	comment("fr's", dat);
	for(int i=0; i<=s.nz+1; ++i) (afr[i]).write(dat);
	comment("bl's", dat);
	for(int i=0; i<=s.nz+1; ++i) (abl[i]).write(dat);
	comment("br's", dat);
	for(int i=0; i<=s.nz+1; ++i) (abr[i]).write(dat);
	return;
}

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

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

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

	fprintf(stderr, ">> %s\n", name);
	dat = fopen(name, "w+");
	write(dat);
	fclose(dat);
	return;
}

/* expansion coefficients: specify amplitudes of forward traveling fields
   on segment seg */
void BepField::setf(int seg, SBorder pos, const Cvector& amp)
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("setf: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("setf ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("setf ( > s.nz+1)");
	if(marr[seg].num != amp.nel) bepfielderror("setf: dimensions"); 
	if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("setf: pos"); 
	if(seg >= nums) bepfielderror("setf: struct mismatch"); 
	if(pos == SBLEFT) afl[seg] = amp;
	else              afr[seg] = amp;
	return; 
}

/* expansion coefficients: specify amplitudes of backward traveling fields
   on segment seg */
void BepField::setb(int seg, SBorder pos, const Cvector& amp)
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("setb: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("setb ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("setb ( > s.nz+1)");
	if(marr[seg].num != amp.nel) bepfielderror("setb: dimensions "); 
	if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("setb: pos"); 
	if(seg >= nums) bepfielderror("setb: struct mismatch"); 
	if(pos == SBLEFT) abl[seg] = amp;
	else              abr[seg] = amp;
	return; 
}

/* expansion coefficients: retrieve amplitudes of forward traveling fields
   on segment seg */
Cvector BepField::f(int seg, SBorder pos) const
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("f: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("f ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("f ( > s.nz+1)");
	if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("f: pos"); 
	if(seg >= nums) bepfielderror("f: struct mismatch"); 
	if(pos == SBLEFT) return afl[seg];
	else              return afr[seg];
}

/* expansion coefficients: retrieve amplitudes of backward traveling fields
   on segment seg */
Cvector BepField::b(int seg, SBorder pos) const
{
	if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) 
		bepfielderror("b: SegWgStruct cleared" ); 
	if(seg < 0)      bepfielderror("b ( < 0)"); 
	if(seg > s.nz+1) bepfielderror("b ( > s.nz+1)");
	if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("b: pos"); 
	if(seg >= nums) bepfielderror("b: struct mismatch"); 
	if(pos == SBLEFT) return abl[seg];
	else              return abr[seg];
}

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

/* check the SegWgStruct for consistency with respect to
   polarization, wavelength, array dimensions; exit, if inconsistent */
void BepField::check()
{
	if(nums != s.nz+2) 
		bepfielderror("check: general array dimensions"); 
/*
	if(marr == NULL) bepfielderror("check: SegWgStruct cleared (marr)"); 
	if(afl  == NULL) bepfielderror("check: SegWgStruct cleared (afl)" ); 
	if(afr  == NULL) bepfielderror("check: SegWgStruct cleared (afr)" ); 
	if(abl  == NULL) bepfielderror("check: SegWgStruct cleared (abl)" ); 
	if(abr  == NULL) bepfielderror("check: SegWgStruct cleared (abr)" ); 
	for(int i=0; i<=s.nz+1; ++i)
	{
		if(afl[i].nel != marr[i].num)
			bepfielderror("check: array dimensions (afl)"); 
		if(afr[i].nel != marr[i].num)
			bepfielderror("check: array dimensions (afr)"); 
		if(abl[i].nel != marr[i].num)
			bepfielderror("check: array dimensions (abl)"); 
		if(abr[i].nel != marr[i].num)
			bepfielderror("check: array dimensions (abr)"); 
		if(afl[i].nel <= 0) bepfielderror("check: afl.nel == 0"); 
		if(afr[i].nel <= 0) bepfielderror("check: afr.nel == 0"); 
		if(abl[i].nel <= 0) bepfielderror("check: abl.nel == 0"); 
		if(abr[i].nel <= 0) bepfielderror("check: abr.nel == 0"); 
	}
*/
	BF_pol = (marr[0])(0).pol;
	BF_lambda = (marr[0])(0).wg.lambda;
	BF_vM = (marr[0])(0).vM;
        BF_ky = (marr[0])(0).ky;

	for(int i=0; i<=s.nz+1; ++i)
	{
		int mni = marr[i].num;
		Mode mij;
		for(int j=0; j<=mni-1; ++j)
		{
			mij = (marr[i])(j);
			if(mij.wg.lambda != BF_lambda) 
				bepfielderror("check: wavelengths mixed");
			if(BF_vM)
			{
			if(mij.vM != BF_vM) 
				bepfielderror("check: nonuniform vM");
			if((mij.ky-BF_ky).sqabs() > SLAMS_BQTOL)
				bepfielderror("check: nonuniform ky");
			}
			else
			{
			if(mij.pol != BF_pol) 
				bepfielderror("check: polarizations mixed");
			}
		}
	}
	return;
}

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

/* mode interference evaluation 
 * all modes on a segment are assumed to belong to the same waveguide ! */

/* field superposition at point x, z 
   cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex BepField::field(Fcomp cp, double x, double z) const
{
	if(s.nz == -1) return 
		marr[0].field(afl[0], abl[0], cp, x, z);
	int zi = s.segidx(z);
	if(zi <= 0) return 
		marr[0].field(afr[0], abr[0], cp, x, z-s.hz(0));
	if(zi >= s.nz+1) return 
		marr[s.nz+1].field(afl[s.nz+1], abl[s.nz+1], 
		                   cp, x, z-s.hz(s.nz));
	return marr[zi].field(afl[zi], abr[zi], cp, 
	                      x, z, s.hz(zi-1), s.hz(zi));
}

/* local field "strength",  id: one of mE, mH, qE, qH, mW */
double BepField::lfs(FSid id, double x, double z) const
{
	int zi = s.segidx(z);
	if(zi <= 0) return marr[0].lfs(afr[0], abr[0], id, x, z-s.hz(0));
	if(zi >= s.nz+1) return 
		marr[s.nz+1].lfs(afl[s.nz+1], abl[s.nz+1], 
		                 id, x, z-s.hz(s.nz));
	return marr[zi].lfs(afl[zi], abr[zi],
	                    id, x, z, s.hz(zi-1), s.hz(zi));
}

/* local electromagnetic energy density at x, z */
double BepField::endens(double x, double z) const
{
	return lfs(mW, x, z);
}

/* the six electric and magnetic components at x, z */
void BepField::emfield(double x, double z,
                       Complex& ex, Complex& ey, Complex& ez,
                       Complex& hx, Complex& hy, Complex& hz) const
{
	if(s.nz == -1)
	{ 
		marr[0].emfield(afl[0], abl[0], x, z, ex, ey, ez, hx, hy, hz);
		return;
	}
	int zi = s.segidx(z);
	if(zi <= 0)
	{ 
		marr[0].emfield(afr[0], abr[0], x, z-s.hz(0), ex, ey, ez, hx, hy, hz);
		return;
	}
	if(zi >= s.nz+1)
	{ 
		marr[s.nz+1].emfield(afl[s.nz+1], abl[s.nz+1], x, z-s.hz(s.nz),
	                             ex, ey, ez, hx, hy, hz);
		return; 
	}
	marr[zi].emfield(afl[zi], abr[zi], x, z, s.hz(zi-1), s.hz(zi),
	                 ex, ey, ez, hx, hy, hz);
	return;
}


/* the total directional power of the mode superposition on segment seg;
   assumptions:
   - mode orthorgonality, different modes of the same waveguide
   - modes are normalized
   - only forward or only backward propagating fields are present */
double BepField::dir_power(int seg, Propdir dir) const
{
	if(BF_vM)
	{ 
		double a, p = 0.0;
		if(seg < 0) seg = 0;
		if(seg > s.nz+1) seg = s.nz+1;
		ModeArray ma = marr[seg];
		int n = ma.num;
		switch(dir)
		{
		case FORW:
			for(int i=0; i<=n-1; ++i)
			{
				if(ma(i).cVppt == PROPAG) 
				{
					if(seg == 0) a = afr[0](i).sqabs();
					else         a = afl[seg](i).sqabs();
					p += a;
				}
			}
			return p;
			break;
		case BACK:
			for(int i=0; i<=n-1; ++i)
			{
				if(ma(i).cVppt == PROPAG) 
				{
					if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
					else              a = abr[seg](i).sqabs();
					p += a;
				}
			}
			return p;
			break;
		default: bepfielderror("dir_power, vM: dir"); break;
		}
		return 0.0;
	}
	double a, p = 0.0;
	if(seg < 0) seg = 0;
	if(seg > s.nz+1) seg = s.nz+1;
	ModeArray ma = marr[seg];
	int n = ma.num;
	switch(dir)
	{
	case FORW:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).ppt == PROPAG) 
			{
				if(seg == 0) a = afr[0](i).sqabs();
				else         a = afl[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	case BACK:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).ppt == PROPAG) 
			{
				if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
				else              a = abr[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	default: bepfielderror("dir_power: dir"); break;
	}
	return 0.0;
}
/* ... modes of polarization pol (vM) */
double BepField::dir_power(int seg, Polarization pol, Propdir dir) const
{
	if(BF_vM == 0) return dir_power(seg, dir);
	double a, p = 0.0;
	if(seg < 0) seg = 0;
	if(seg > s.nz+1) seg = s.nz+1;
	ModeArray ma = marr[seg];
	int n = ma.num;
	switch(dir)
	{
	case FORW:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).pol == pol && ma(i).cVppt == PROPAG) 
			{
				if(seg == 0) a = afr[0](i).sqabs();
				else         a = afl[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	case BACK:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).pol == pol && ma(i).cVppt == PROPAG) 
			{
				if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
				else              a = abr[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	default: bepfielderror("dir_power, pol: dir"); break;
	}
	return 0.0;
}

/* the guided directional power of the mode superposition on segment seg;
   assumptions:
   - mode orthorgonality, different modes of the same waveguide
   - modes are normalized
   - only forward or only backward propagating fields are present */
double BepField::dir_guided_power(int seg, Propdir dir) const
{
	if(BF_vM)
	{ 
		double a, p = 0.0;
		if(seg < 0) seg = 0;
		if(seg > s.nz+1) seg = s.nz+1;
		ModeArray ma = marr[seg];
		int n = ma.num;
		switch(dir)
		{
		case FORW:
			for(int i=0; i<=n-1; ++i)
			{
				if(ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG)
				{
					if(seg == 0) a = afr[0](i).sqabs();
					else         a = afl[seg](i).sqabs();
					p += a;
				}
			}
			return p;
			break;
		case BACK:
			for(int i=0; i<=n-1; ++i)
			{
				if(ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG)
				{
					if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
					else              a = abr[seg](i).sqabs();
					p += a;
				}
			}
			return p;
			break;
		default: bepfielderror("cVdir_guided_power: dir"); break;
		}
		return 0.0;
	}
	double a, p = 0.0;
	if(seg < 0) seg = 0;
	if(seg > s.nz+1) seg = s.nz+1;
	ModeArray ma = marr[seg];
	int n = ma.num;
	switch(dir)
	{
	case FORW:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).npcB > 0.0 && ma(i).ppt == PROPAG)
			{
				if(seg == 0) a = afr[0](i).sqabs();
				else         a = afl[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	case BACK:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).npcB > 0.0 && ma(i).ppt == PROPAG)
			{
				if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
				else              a = abr[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	default: bepfielderror("dir_guided_power: dir"); break;
	}
	return 0.0;
}
/* ... of polarization pol (vM) */
double BepField::dir_guided_power(int seg, Polarization pol, Propdir dir) const
{
	if(BF_vM == 0) return dir_guided_power(seg, dir);
	double a, p = 0.0;
	if(seg < 0) seg = 0;
	if(seg > s.nz+1) seg = s.nz+1;
	ModeArray ma = marr[seg];
	int n = ma.num;
	switch(dir)
	{
	case FORW:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).pol == pol && ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG)
			{
				if(seg == 0) a = afr[0](i).sqabs();
				else         a = afl[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	case BACK:
		for(int i=0; i<=n-1; ++i)
		{
			if(ma(i).pol == pol && ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG)
			{
				if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs();
				else              a = abr[seg](i).sqabs();
				p += a;
			}
		}
		return p;
		break;
	default: bepfielderror("dir_guided_power, pol: dir"); break;
	}
	return 0.0;
}


/* --- specification of input fields, additive ----------------------- */

/* mode m with amplitude A at border bdr */
void BepField::input(SBorder bdr, Complex A, const Mode& m)
{
	if(m.ppt != PROPAG) bepfielderror("input, EVANESC mode");
	Cvector amp;
	ModeArray ma;
	int n;
	switch(bdr)
	{
	case SBRIGHT:  
		amp = b(s.nz+1, SBLEFT); 
		ma = marr[s.nz+1];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				amp(j) -= overlap(ma(j), BACK, m, BACK);
		}
		setb(s.nz+1, SBLEFT, amp);
		return;
	case SBLEFT:  
		amp = f(0, SBRIGHT); 
		ma = marr[0];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				amp(j) += overlap(ma(j), FORW, m, FORW);
		}
		setf(0, SBRIGHT, amp);
		return;
	}
	return;
}
/* mode of order o with amplitude A at border bdr */
void BepField::input(SBorder bdr, Complex A, int o)
{
	Cvector amp;
	switch(bdr)
	{
	case SBRIGHT:  
		amp = b(s.nz+1, SBLEFT); 
		if(o<0 || o>amp.nel-1) bepfielderror("input: order, SBRIGHT");
		if((marr[s.nz+1])(o).ppt == EVANESC) 
			bepfielderror("input, EVANESC mode, SBRIGHT");
		amp(o) += A;
		setb(s.nz+1, SBLEFT, amp);
		return;
	case SBLEFT:  
		amp = f(0, SBRIGHT); 
		if(o<0 || o>amp.nel-1) bepfielderror("input: order, SBLEFT");
		if((marr[0])(o).ppt == EVANESC) 
			bepfielderror("input, EVANESC mode, SBLEFT");
		amp(o) += A;
		setf(0, SBRIGHT, amp);
		return;
	}
	return;
}
/* Gaussian beam gb, multiplied by amplitude A, at border bdr */
void BepField::input(SBorder bdr, const Gaussianbeam& gb, Complex A)
{
	Cvector amp;
	ModeArray ma;
	int n;
	switch(bdr)
	{
	case SBRIGHT:  
		amp = b(s.nz+1, SBLEFT); 
		ma = marr[s.nz+1];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				amp(j) -= A*(ma(j)).overlap(gb, BACK);
		}
		setb(s.nz+1, SBLEFT, amp);
		return;
	case SBLEFT:  
		amp = f(0, SBRIGHT); 
		ma = marr[0];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				amp(j) += A*(ma(j)).overlap(gb, FORW);
		}
		setf(0, SBRIGHT, amp);
		return;
	}
	return;
}
/* ... a normalized beam, unspecified phase */
void BepField::input(SBorder bdr, const Gaussianbeam& gb)
{
	input(bdr, gb, CC1);
	return;
}

/* vM: amplitude of incoming mode of order o passing border bdr, 
   with polarization pol, evaluated at propagation coordinate p */
Complex BepField::Ain(SBorder bdr, Polarization pol, int o, double p) const
{
	if(BF_vM == 0) bepfielderror("Ain, vM = 0, pol, pos: not implemented"); 
	int num, mx, mi; 
	double zw;
	Complex a, ef;
	switch(bdr)
	{
	case SBRIGHT:  
		num = marr[s.nz+1].num; 
		if(num/2 != ((int) ((double) num)/2.0))
			bepfielderror("Ain, pol, o, p; confused");
		mx = num/2;	
		if(o<0 || o>mx-1) 
			bepfielderror("Ain, pol, o, p; RIGHT: val o");
		if(pol == TE) mi = o;
		else          mi = mx+o;
		if(marr[s.nz+1](mi).pol != pol)
			bepfielderror("Ain, pol, o, p; RIGHT: pol");
		if(marr[s.nz+1](mi).ord != o)
			bepfielderror("Ain, pol, o, p; RIGHT: ord");
		a = (b(s.nz+1, SBLEFT))(mi); 
		zw = s.hz(s.nz);
		ef = (marr[s.nz+1](mi).efac(BACK, p-zw));
		return a*ef;
	case SBLEFT:  
		num = marr[0].num; 
		if(num/2 != ((int) ((double) num)/2.0))
			bepfielderror("Ain, pol, o, p; confused");
		mx = num/2;	
		if(o<0 || o>mx-1) 
			bepfielderror("Ain, pol, o, p; LEFT: val o");
		if(pol == TE) mi = o;
		else          mi = mx+o;
		if(marr[0](mi).pol != pol)
			bepfielderror("Ain, pol, o, p; LEFT: pol");
		if(marr[0](mi).ord != o)
			bepfielderror("Ain, pol, o, p; LEFT: ord");
		a = (f(0, SBRIGHT))(mi); 
		zw = s.hz(0);
		ef = (marr[0](mi).efac(FORW, p-zw));
		return a*ef;
	}
	return CC0;
}

/* --- output amplitudes & powers, assuming normalized modes --- */

/* amplitude of outgoing mode m at border bdr */
Complex BepField::Aout(SBorder bdr, const Mode& m) const 
{
	if(BF_vM)
	{
		if(m.vM == 0) bepfielderror("Aout, vM, m: vM = 0");
		if(m.cVppt != PROPAG) bepfielderror("Aout, vM, m: EVANESC mode");
		Complex a = CC0;
		Cvector amp;
		ModeArray ma;
		int n;
		switch(bdr)
		{
		case SBRIGHT:  
			amp = f(s.nz+1, SBLEFT); 
			ma = marr[s.nz+1];
			n = ma.num;
			for(int j=0; j<=n-1; ++j)
			{
				if(ma(j).cVppt == PROPAG)
				a += amp(j)*overlap(m, FORW, ma(j), FORW);
			}
			return a/overlap(m, FORW, m, FORW);
		case SBLEFT:  
			amp = b(0, SBRIGHT); 
			ma = marr[0];
			n = ma.num;
			for(int j=0; j<=n-1; ++j)
			{
				if(ma(j).cVppt == PROPAG)
				a -= amp(j)*overlap(m, BACK, ma(j), BACK);
			}
			return a/overlap(m, BACK, m, BACK);
		}
		return CC0;
	} 
	if(m.ppt != PROPAG) bepfielderror("Aout, EVANESC mode");
	Complex a = CC0;
	Cvector amp;
	ModeArray ma;
	int n;
	switch(bdr)
	{
	case SBRIGHT:  
		amp = f(s.nz+1, SBLEFT); 
		ma = marr[s.nz+1];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				a += amp(j)*overlap(m, FORW, ma(j), FORW);
		}
		return a/m.power();
	case SBLEFT:  
		amp = b(0, SBRIGHT); 
		ma = marr[0];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			if(ma(j).ppt == PROPAG)
				a -= amp(j)*overlap(m, BACK, ma(j), BACK);
		}
		return a/m.power();
	}
	return CC0;
}
/* amplitude of outgoing mode of order o at border bdr */
Complex BepField::Aout(SBorder bdr, int o) const 
{
	if(BF_vM) bepfielderror("Aout, o: BF_vM");
	switch(bdr)
	{
	case SBRIGHT:  
		if(o<0 || o>marr[s.nz+1].num-1) 
			bepfielderror("Aout: order, RIGHT");
		return (f(s.nz+1, SBLEFT))(o); 
	case SBLEFT:  
		if(o<0 || o>marr[0].num-1) 
			bepfielderror("Aout: order, LEFT");
		return (b(0, SBRIGHT))(o); 
	}
	return CC0;
}

/* amplitude of outgoing mode 
   of polarzation pol, order o at border bdr (vM) */
Complex BepField::Aout(SBorder bdr, Polarization pol, int o) const 
{
	if(BF_vM == 0) bepfielderror("Aout, vM = 0, pol: not implemented"); 
	switch(bdr)
	{
	case SBRIGHT: return Aout(bdr, pol, o, s.hz(s.nz));
	case SBLEFT:  return Aout(bdr, pol, o, s.hz(0));
	}
	return CC0;
}

/* vM: amplitude of outgoing mode of order o passing border bdr, 
   with polarization pol, evaluated at propagation coordinate p */
Complex BepField::Aout(SBorder bdr, Polarization pol, int o, double p) const
{
	if(BF_vM == 0) bepfielderror("Aout, vM = 0, pol, pos: not implemented"); 
	int num, mx, mi; 
	double zw;
	Complex a, ef;
	switch(bdr)
	{
	case SBRIGHT:  
		num = marr[s.nz+1].num; 
		if(num/2 != ((int) ((double) num)/2.0))
			bepfielderror("Aout, pol, o, p; confused");
		mx = num/2;	
		if(o<0 || o>mx-1) 
			bepfielderror("Aout, pol, o, p; RIGHT: val o");
		if(pol == TE) mi = o;
		else          mi = mx+o;
		if(marr[s.nz+1](mi).pol != pol)
			bepfielderror("Aout, pol, o, p; RIGHT: pol");
		if(marr[s.nz+1](mi).ord != o)
			bepfielderror("Aout, pol, o, p; RIGHT: ord");
		if(marr[s.nz+1](mi).cVppt != PROPAG) return CC0;
		a = (f(s.nz+1, SBLEFT))(mi); 
		zw = s.hz(s.nz);
		ef = (marr[s.nz+1](mi).efac(FORW, p-zw));
		return a*ef;
	case SBLEFT:  
		num = marr[0].num; 
		if(num/2 != ((int) ((double) num)/2.0))
			bepfielderror("Aout, pol, o, p; confused");
		mx = num/2;	
		if(o<0 || o>mx-1) 
			bepfielderror("Aout, pol, o, p; LEFT: val o");
		if(pol == TE) mi = o;
		else          mi = mx+o;
		if(marr[0](mi).pol != pol)
			bepfielderror("Aout, pol, o, p; LEFT: pol");
		if(marr[0](mi).ord != o)
			bepfielderror("Aout, pol, o, p; LEFT: ord");
		if(marr[0](mi).cVppt != PROPAG) return CC0;
		a = (b(0, SBRIGHT))(mi); 
		zw = s.hz(0);
		ef = (marr[0](mi).efac(BACK, p-zw));
		return a*ef;
	}
	return CC0;
}



/* output power associated with mode m at border bdr */
double BepField::Pout(SBorder bdr, const Mode& m) const 
{
	return Aout(bdr, m).sqabs(); 
}

/* output power associated with mode of order o at border bdr */
double BepField::Pout(SBorder bdr, int o) const 
{
	return Aout(bdr, o).sqabs();
}

/* output power associated with mode 
   of Polarization pol, order o at border bdr (vM) */
double BepField::Pout(SBorder bdr, Polarization pol, int o) const 
{
	return Aout(bdr, pol, o).sqabs();
}

/* guided output power crossing border bdr */
double BepField::Pgout(SBorder bdr) const 
{
	switch(bdr)
	{
	case SBRIGHT:  return dir_guided_power(s.nz+1, FORW);
	case SBLEFT:   return dir_guided_power(0, BACK);
	}
	return 0.0;
}

/* guided output power associated with polarization pol (vM)
   crossing border bdr */
double BepField::Pgout(SBorder bdr, Polarization pol) const 
{
	switch(bdr)
	{
	case SBRIGHT:  return dir_guided_power(s.nz+1, pol, FORW);
	case SBLEFT:   return dir_guided_power(0, pol, BACK);
	}
	return 0.0;
}


/* total output power crossing border bdr (guided and nonguided fields) */
double BepField::Pout(SBorder bdr) const 
{
	switch(bdr)
	{
	case SBRIGHT:  return dir_power(s.nz+1, FORW);
	case SBLEFT:   return dir_power(0, BACK);
	}
	return 0.0;
}

/* total output power associated with polarization pol
   crossing border bdr (guided and nonguided fields) */
double BepField::Pout(SBorder bdr, Polarization pol) const 
{
	switch(bdr)
	{
	case SBRIGHT:  return dir_power(s.nz+1, pol, FORW);
	case SBLEFT:   return dir_power(0, pol, BACK);
	}
	return 0.0;
}


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

/* total input power crossing border bdr (guided and nonguided fields) */
double BepField::Pin(SBorder bdr) const 
{
	switch(bdr)
	{
	case SBLEFT:   return dir_power(0, FORW);
	case SBRIGHT:  return dir_power(s.nz+1, BACK);
	}
	return 0.0;
}

/* --- "modal probing": projection of a mode onto the local expansion ---- */ 

/* directional amplitude of probe mode m at position z */
Complex BepField::locModAmp(double z, const Mode& m, Propdir dir) const 
{
	if(BF_vM) bepfielderror("locModAmp: vM not implemented");
	Complex a = CC0;
	Cvector famp, bamp;
	ModeArray ma;
	int n;
	int si = s.segidx(z);
	Interval seg = s.segment(z);
	if(si <= 0)
	{
		famp = f(0, SBRIGHT); 
		bamp = b(0, SBRIGHT); 
		ma = marr[0];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			a += famp(j)*overlap(m, dir, ma(j), FORW)
				    *ma(j).efac(FORW, z-seg.x1);
			a += bamp(j)*overlap(m, dir, ma(j), BACK)
				    *ma(j).efac(BACK, z-seg.x1);
		}
		return a/m.power();
	}
	if(si >= s.nz+1)
	{
		famp = f(s.nz+1, SBLEFT); 
		bamp = b(s.nz+1, SBLEFT); 
		ma = marr[s.nz+1];
		n = ma.num;
		for(int j=0; j<=n-1; ++j)
		{
			a += famp(j)*overlap(m, dir, ma(j), FORW)
				    *ma(j).efac(FORW, z-seg.x0);
			a += bamp(j)*overlap(m, dir, ma(j), BACK)
				    *ma(j).efac(BACK, z-seg.x0);
		}
		return a/m.power();
	}
	famp = f(si, SBLEFT); 
	bamp = b(si, SBRIGHT); 
	ma = marr[si];
	n = ma.num;
	for(int j=0; j<=n-1; ++j)
	{
		a += famp(j)*overlap(m, dir, ma(j), FORW)
		            *ma(j).efac(FORW, z-seg.x0);
		a += bamp(j)*overlap(m, dir, ma(j), BACK)
		            *ma(j).efac(BACK, z-seg.x1);
	}
	return a/m.power();
}

/* power associated with directional mode m at position z */
double BepField::locModPow(double z, const Mode& m, Propdir dir) const 
{
	if(BF_vM) bepfielderror("locModPow: vM not implemented");
	return locModAmp(z, m, dir).sqabs()*m.power(); 
}

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

/* compute the discretized bidirectional mode spectra:
   modes of polarization p,
   Mx expansion terms per vertical slice on a vertical 
   computational window cw */
void BepField::discretize(Polarization p, Interval cw, int Mx, int quiet)
{
	Waveguide wg;
	ModeArray ma;
	if(quiet != 1) fprintf(stderr, "SpecDisc: ");
	for(int i=0; i<=s.nz+1; ++i)
	{	
		int done = 0; 
		int j=0;
		while(j<i && done==0)
		{
			if(s(i).equal(s(j))==1)
			{
				marr[i] = marr[j];
				done = 1;
			}
			++j;
		}
		if(done==0)
		{
			wg = s(i);
			modespectrum(wg, p, cw, Mx, ma, 1);
			marr[i] = ma;
			if(quiet != 1) fprintf(stderr, ".");
		}
	}
	if(quiet != 1) fprintf(stderr, " Ok.\n");
	Cvector zero;
	zero = Cvector(marr[0].num); zero.init(CC0);
	setf(0, SBRIGHT, zero);
	zero = Cvector(marr[s.nz+1].num); zero.init(CC0);
	setb(s.nz+1, SBLEFT, zero);
	check();
	return;
}


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

/* mode expansion solver:
   given the structure with the accompanying mode expansion,
   where input mode amplitudes are specified by setf(0, f),
   find the remaining mode coefficients.
   Assumption: no backwards traveling fields on the last segment
   tm_m, tm_p: already computed arrays of transfermatrices 
               for the juntions in this object, 
               as output from transfermats() */
void BepField::bepsim(Cmatrix *tm_m, Cmatrix *tm_p)
{
	fprintf(stderr, "BEP (");
	if(s.nz < -1 ) bepfielderror("bepsim: nz < -1"); 
	if(s.nz+2 != nums) bepfielderror("bepsim: structure mismatch"); 
	fprintf(stderr, "%d x", s.nz+2);
	int rfmin, rfmax;
	if(afr[0].nel != marr[0].num) bepfielderror("bepsim: dimensions 0"); 
	abr[0] = Cvector(marr[0].num);
	afl[0] = Cvector(marr[0].num);
	abl[0] = Cvector(marr[0].num);
	rfmin = rfmax = marr[0].num;
	for(int i=1; i<=s.nz+1; ++i)
	{
		afl[i] = Cvector(marr[i].num);
		afr[i] = Cvector(marr[i].num);
		abl[i] = Cvector(marr[i].num);
		abr[i] = Cvector(marr[i].num);
		if(rfmin > marr[i].num) rfmin = marr[i].num;
		if(rfmax < marr[i].num) rfmax = marr[i].num;
	}
	if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin);
	else fprintf(stderr, " (%d-%d)) ", rfmin, rfmax);

	if(s.nz == -1)
	{
		Cvector zero(marr[0].num);
		zero.init(CC0);
		afl[0] = afr[0];
		afr[0] = zero;
		abl[0] = zero;
		abr[0] = zero;
		return;
	}

	Dvector *p;
	p = new Dvector[s.nz+2];
	for(int i=0; i<=s.nz+1; ++i)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=i-1)
		{
			if(marr[o].equal(marr[i]) == 1) 
			{
				p[i] = p[o];
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			ModeArray ma = marr[i];
			int n = ma.num;
			Dvector tp(n);
			for(int m=0; m<=n-1; ++m)
				tp(m) = sqrt(fabs(ma(m).power()));
			p[i] = tp;
		}
	}

	Cmatrix *dmat;
	Ivector filld(s.nz+2);
	Ivector idd(s.nz+2);

	dmat = new Cmatrix[s.nz+2];
	filld.init(0);

	for(int j=0; j<=s.nz+1; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o].equal(marr[j]) == 1 && filld(o) == 1)
			{
				idd(j) = o;
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			dmat[j] = orthomat_inv(marr[j]);
			fprintf(stderr, ".");
			idd(j) = j;
			filld(j) = 1;
		}
	}
	filld.strip();

	Cmatrix m_m, m_p;
	Cmatrix *mmff, *mmfb, *mmbf, *mmbb;
	Cmatrix *mpff, *mpfb, *mpbf, *mpbb;
	Ivector fill(s.nz+1);
	Ivector id(s.nz+1);

	mmff = new Cmatrix[s.nz+1];
	mmfb = new Cmatrix[s.nz+1];
	mmbf = new Cmatrix[s.nz+1];
	mmbb = new Cmatrix[s.nz+1];
	mpff = new Cmatrix[s.nz+1];
	mpfb = new Cmatrix[s.nz+1];
	mpbf = new Cmatrix[s.nz+1];
	mpbb = new Cmatrix[s.nz+1];
	fill.init(0);

	for(int j=0; j<=s.nz; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o  ].equal(marr[j  ]) == 1 &&
			   marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1)
			{
				id(j) = o+1;
				done = 1;
			}
			if(marr[o  ].equal(marr[j+1]) == 1 &&
			   marr[o+1].equal(marr[j  ]) == 1 && fill(o) == 1)
			{
				id(j) = -1*(o+1);
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			int nl, nr;
			nl = marr[j  ].num;
			nr = marr[j+1].num;
			if(tm_m != NULL && tm_p != NULL) 
			{
				m_m = tm_m[j];
				m_p = tm_p[j];
			}
			else 
			transfermats(marr[j], marr[j+1], 
			             dmat[idd(j)], dmat[idd(j+1)], m_m, m_p);
			fprintf(stderr, ".");
			mmff[j] = m_m.submatrix(nl, nr, 0,  0 );
			mmfb[j] = m_m.submatrix(nl, nr, 0,  nr);
			mmbf[j] = m_m.submatrix(nl, nr, nl, 0 );
			mmbb[j] = m_m.submatrix(nl, nr, nl, nr);
			m_m.strip();
			mpff[j] = m_p.submatrix(nr, nl, 0,  0 );
			mpfb[j] = m_p.submatrix(nr, nl, 0,  nl);
			mpbf[j] = m_p.submatrix(nr, nl, nr, 0 );
			mpbb[j] = m_p.submatrix(nr, nl, nr, nl);
			m_p.strip();
			id(j) = j+1;
			fill(j) = 1;
		}
	}
	fill.strip();
	idd.strip();
	delete[] dmat;

	fprintf(stderr, ":");

	Cvector *t = NULL;
	if(s.nz >= 1)
	{
		t = new Cvector[s.nz];
		for(int i=0; i<=s.nz-1; ++i)
		{
			ModeArray ma = marr[i+1];
			int ni = ma.num;
			double l = s.hz(i+1)-s.hz(i);
			Cvector tmpt(ni);
			for(int j=0; j<=ni-1; ++j)
			{
/*
				if(ma(j).ppt == PROPAG)
				   tmpt(j) = expi(-l*ma(j).beta);
				else
				   tmpt(j) = exp(-l*ma(j).beta);
*/
				tmpt(j) = ma(j).efac(FORW, l);
			}
			t[i] = tmpt;
		}
		fprintf(stderr, ".");
	}

	Cmatrix u;
	Cmatrix ri, rr;
	Cmatrix *r;
	r = new Cmatrix[s.nz+1];

	u.unit(marr[s.nz].num);
	if(id(s.nz) > 0) 
		ri = sub(u, mult(mmbf[ id(s.nz)-1], mpfb[ id(s.nz)-1]));
	else    
	        ri = sub(u, mult(mpbf[-id(s.nz)-1], mmfb[-id(s.nz)-1]));
	ri.inverse();
	if(id(s.nz) > 0) 
		r[s.nz] = mult(ri, mult(mmbf[ id(s.nz)-1], mpff[ id(s.nz)-1]));
	else 
		r[s.nz] = mult(ri, mult(mpbf[-id(s.nz)-1], mmff[-id(s.nz)-1]));
	fprintf(stderr, ".");

	for(int j=s.nz-1; j>=0; --j)
	{
		u.unit(marr[j].num);
		if(id(j) > 0) 
			rr = add(mmbf[ id(j)-1], mult(mmbb[ id(j)-1], 
		                 diagmult(t[j], diagmult(r[j+1],t[j])))); 
		else
			rr = add(mpbf[-id(j)-1], mult(mpbb[-id(j)-1], 
		                 diagmult(t[j], diagmult(r[j+1],t[j])))); 
		if(id(j) > 0) ri = sub(u, mult(rr, mpfb[ id(j)-1])); 
		else          ri = sub(u, mult(rr, mmfb[-id(j)-1])); 
		ri.inverse();
		if(id(j) > 0) r[j] = mult(ri, mult(rr, mpff[ id(j)-1]));
		else          r[j] = mult(ri, mult(rr, mmff[-id(j)-1]));
		rr.strip();
		ri.strip();
		fprintf(stderr, ".");
	}
	u.strip();

	delete[] mpbf;
	delete[] mpbb;
	delete[] mmbf;
	delete[] mmbb;

	for(int m=0; m<=marr[0].num-1; ++m) (afr[0])(m)/=(p[0])(m);

	abr[0] = mult(r[0], afr[0]);
	r[0].strip();

	for(int j=1; j<=s.nz; ++j)
	{
		if(id(j-1) > 0)
			afl[j] = add(mult(mpff[ id(j-1)-1], afr[j-1]), 
				     mult(mpfb[ id(j-1)-1], abr[j-1]));
		else
			afl[j] = add(mult(mmff[-id(j-1)-1], afr[j-1]), 
				     mult(mmfb[-id(j-1)-1], abr[j-1]));
		afr[j] = diagmult(t[j-1], afl[j]);
		abr[j] = mult(r[j],   afr[j]); 
		r[j].strip();
		abl[j] = diagmult(t[j-1], abr[j]); 
		fprintf(stderr, ".");
	}

	if(id(s.nz) > 0)
		afl[s.nz+1] = add(mult(mpff[ id(s.nz)-1], afr[s.nz]), 
				  mult(mpfb[ id(s.nz)-1], abr[s.nz]));
	else
		afl[s.nz+1] = add(mult(mmff[-id(s.nz)-1], afr[s.nz]), 
				  mult(mmfb[-id(s.nz)-1], abr[s.nz]));

	delete[] mpff;
	delete[] mpfb;
	delete[] mmff;
	delete[] mmfb;
	delete[] r;

	for(int m=0; m<=marr[0].num-1; ++m)
	{ 
		(afr[0])(m)*=(p[0])(m);
		(abr[0])(m)*=(p[0])(m);
	}
	for(int j=1; j<=s.nz; ++j)
	{
		for(int m=0; m<=marr[j].num-1; ++m) 
		{
			(afl[j])(m)*=(p[j])(m);
			(afr[j])(m)*=(p[j])(m);
			(abl[j])(m)*=(p[j])(m);
			(abr[j])(m)*=(p[j])(m);
		}
	}
	for(int m=0; m<=marr[s.nz+1].num-1; ++m) 
		(afl[s.nz+1])(m)*=(p[s.nz+1])(m);

	fprintf(stderr, " Ok.\n");

	Cvector zero;
	zero = Cvector(marr[0].num);
	zero.init(CC0);
	abl[0] = zero;
	afl[0] = zero;
	zero = Cvector(marr[s.nz+1].num);
	zero.init(CC0);
	abl[s.nz+1] = zero;
	afr[s.nz+1] = zero;
	abr[s.nz+1] = zero;

	delete[] p;
	if(t != NULL) delete[] t;

	return;
}

/* mode expansion solver:
   given the structure with the accompanying mode expansion,
   where input mode amplitudes are specified by setf(0, f),
   find the remaining mode coefficients.
   Assumption: no backwards traveling fields on the last segment */
void BepField::bepsim() 
{
	Cmatrix *tm_m;	
	Cmatrix *tm_p;	
	tm_m = NULL;
	tm_p = NULL; 
	bepsim(tm_m, tm_p);	
	return;
}

/* mode expansion solver:
   given the structure with the accompanying mode expansion and
   input mode amplitudes inp on segment 0, 
   find the amplitudes ref and trans of the reflected and transmitted fields; 
   fill the intermediate amplitudes.
   Assumption: no backwards traveling fields on the last segment */
void BepField::bepsim(Cvector& inp, Cvector& ref, Cvector& trans)
{
	setf(0, SBRIGHT, inp);
	bepsim();
	ref   = b(0, SBRIGHT);
	trans = f(s.nz+1, SBLEFT);
	return;
}

/* mode expansion solver:
   given the structure with the accompanying mode expansion
   find the mode reflectivity matrix ref
   and transmittivity matrix trans
   the mode amplitude vectors are not accessed
   Assumption: no backwards traveling fields on the last segment */
void BepField::bepsim(Cmatrix& ref, Cmatrix& trans)    
{
	fprintf(stderr, "BEP (");
	if(s.nz < -1 ) bepfielderror("bepsim: nz < -1"); 
	if(s.nz+2 != nums) bepfielderror("bepsim: structure mismatch"); 
	fprintf(stderr, "%d x", s.nz+2);
	int rfmin, rfmax;
	rfmin = rfmax = marr[0].num;
	for(int i=1; i<=s.nz+1; ++i)
	{
		if(rfmin > marr[i].num) rfmin = marr[i].num;
		if(rfmax < marr[i].num) rfmax = marr[i].num;
	}
	if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin);
	else fprintf(stderr, " (%d-%d)) ", rfmin, rfmax);

	if(s.nz == -1)
	{
		ref   = Cmatrix(marr[0].num, marr[0].num);
		ref.init(CC0);
		trans.unit(marr[0].num);
		return;
	}

	Cmatrix *dmat;
	Ivector filld(s.nz+2);
	Ivector idd(s.nz+2);

	dmat = new Cmatrix[s.nz+2];
	filld.init(0);

	for(int j=0; j<=s.nz+1; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o].equal(marr[j]) == 1 && filld(o) == 1)
			{
				idd(j) = o;
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			dmat[j] = orthomat_inv(marr[j]);
			fprintf(stderr, ".");
			idd(j) = j;
			filld(j) = 1;
		}
	}
	filld.strip();

	Cmatrix m_m, m_p;
	Cmatrix *mmff, *mmfb, *mmbf, *mmbb;
	Cmatrix *mpff, *mpfb, *mpbf, *mpbb;
	Ivector fill(s.nz+1);
	Ivector id(s.nz+1);

	mmff = new Cmatrix[s.nz+1];
	mmfb = new Cmatrix[s.nz+1];
	mmbf = new Cmatrix[s.nz+1];
	mmbb = new Cmatrix[s.nz+1];
	mpff = new Cmatrix[s.nz+1];
	mpfb = new Cmatrix[s.nz+1];
	mpbf = new Cmatrix[s.nz+1];
	mpbb = new Cmatrix[s.nz+1];
	fill.init(0);

	for(int j=0; j<=s.nz; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o  ].equal(marr[j  ]) == 1 &&
			   marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1)
			{
				id(j) = o+1;
				done = 1;
			}
			if(marr[o  ].equal(marr[j+1]) == 1 &&
			   marr[o+1].equal(marr[j  ]) == 1 && fill(o) == 1)
			{
				id(j) = -1*(o+1);
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			int nl, nr;
			nl = marr[j  ].num;
			nr = marr[j+1].num;
			transfermats(marr[j], marr[j+1], 
			             dmat[idd(j)], dmat[idd(j+1)], m_m, m_p);
			fprintf(stderr, ".");
			mmff[j] = m_m.submatrix(nl, nr, 0,  0 );
			mmfb[j] = m_m.submatrix(nl, nr, 0,  nr);
			mmbf[j] = m_m.submatrix(nl, nr, nl, 0 );
			mmbb[j] = m_m.submatrix(nl, nr, nl, nr);
			m_m.strip();
			mpff[j] = m_p.submatrix(nr, nl, 0,  0 );
			mpfb[j] = m_p.submatrix(nr, nl, 0,  nl);
			mpbf[j] = m_p.submatrix(nr, nl, nr, 0 );
			mpbb[j] = m_p.submatrix(nr, nl, nr, nl);
			m_p.strip();
			id(j) = j+1;
			fill(j) = 1;
		}
	}
	fill.strip();
	idd.strip();
	delete[] dmat;

	fprintf(stderr, ":");

	Cvector *t = NULL;
	if(s.nz >= 1)
	{
		t = new Cvector[s.nz];
		for(int i=0; i<=s.nz-1; ++i)
		{
			ModeArray ma = marr[i+1];
			int ni = ma.num;
			double l = s.hz(i+1)-s.hz(i);
			Cvector tmpt(ni);
			for(int j=0; j<=ni-1; ++j)
			{
/*
				if(ma(j).ppt == PROPAG)
				   tmpt(j) = expi(-l*ma(j).beta);
				else
				   tmpt(j) = exp(-l*ma(j).beta);
*/
				tmpt(j) = ma(j).efac(FORW, l);
			}
			t[i] = tmpt;
		}
		fprintf(stderr, ".");
	}

	Cmatrix u;
	Cmatrix ri, rr;
	Cmatrix *r;
	r = new Cmatrix[s.nz+1];

	u.unit(marr[s.nz].num);
	if(id(s.nz) > 0) 
		ri = sub(u, mult(mmbf[ id(s.nz)-1], mpfb[ id(s.nz)-1]));
	else    
	        ri = sub(u, mult(mpbf[-id(s.nz)-1], mmfb[-id(s.nz)-1]));
	ri.inverse();
	if(id(s.nz) > 0) 
		r[s.nz] = mult(ri, mult(mmbf[ id(s.nz)-1], mpff[ id(s.nz)-1]));
	else 
		r[s.nz] = mult(ri, mult(mpbf[-id(s.nz)-1], mmff[-id(s.nz)-1]));
	fprintf(stderr, ".");

	for(int j=s.nz-1; j>=0; --j)
	{
		u.unit(marr[j].num);
		if(id(j) > 0) 
			rr = add(mmbf[ id(j)-1], mult(mmbb[ id(j)-1], 
		                 diagmult(t[j], diagmult(r[j+1],t[j])))); 
		else
			rr = add(mpbf[-id(j)-1], mult(mpbb[-id(j)-1], 
		                 diagmult(t[j], diagmult(r[j+1],t[j])))); 
		if(id(j) > 0) ri = sub(u, mult(rr, mpfb[ id(j)-1])); 
		else          ri = sub(u, mult(rr, mmfb[-id(j)-1])); 
		ri.inverse();
		if(id(j) > 0) r[j] = mult(ri, mult(rr, mpff[ id(j)-1]));
		else          r[j] = mult(ri, mult(rr, mmff[-id(j)-1]));
		rr.strip();
		ri.strip();
		fprintf(stderr, ".");
	}
	u.strip();

	delete[] mpbf;
	delete[] mpbb;
	delete[] mmbf;
	delete[] mmbb;


	ref = r[0];
	Cmatrix tf, tb, tmp;
	
	tf.unit(marr[0].num);
	tb = r[0];
	r[0].strip();

	for(int j=1; j<=s.nz; ++j)
	{
		if(id(j-1) > 0)
			tmp = add(mult(mpff[ id(j-1)-1], tf), 
			          mult(mpfb[ id(j-1)-1], tb));
		else
			tmp = add(mult(mmff[-id(j-1)-1], tf), 
		  	          mult(mmfb[-id(j-1)-1], tb));
		tf = diagmult(t[j-1], tmp);
		tb = mult(r[j], tf); 
		r[j].strip();
		fprintf(stderr, ".");
	}

	if(id(s.nz) > 0)
		trans = add(mult(mpff[ id(s.nz)-1], tf), 
		            mult(mpfb[ id(s.nz)-1], tb));
	else
		trans = add(mult(mmff[-id(s.nz)-1], tf), 
		            mult(mmfb[-id(s.nz)-1], tb));

	tf.strip();
	tb.strip();
	delete[] mpff;
	delete[] mpfb;
	delete[] mmff;
	delete[] mmfb;
	delete[] r;
	fprintf(stderr, " Ok.\n");

	if(t != NULL) delete[] t;

	return;
}

/* mode expansion solver:
   given the structure with the accompanying mode expansion
   find mode reflectivity and transmissivity matrices kff, kfb, kbf, kbb, 
   that link the amplitudes b(0, SBRIGHT), f(s.nz+1, SBLEFT) of outgoing modes 
   to the incoming mode amplitudes f(0, SBRIGHT), b(s.nz+1, SBLEFT) as 
   ( f(s.nz+1, SBLEFT) )   ( kff  kfb ) ( f(0, SBRIGHT)     )
   (                   ) = (          ) (                   )
   ( b(0, SBRIGHT)     )   ( kbf  kbb ) ( b(s.nz+1, SBLEFT) ) ;
   the mode amplitude vectors are not accessed */
void BepField::bepsim(Cmatrix& kff, Cmatrix& kfb, Cmatrix& kbf, Cmatrix& kbb)
{
	fprintf(stderr, "BEP (");
	if(s.nz+2 != nums) bepfielderror("bepsim: structure mismatch"); 
	if(s.nz < 1 ) bepfielderror("bepsim: nz < 1"); 
	fprintf(stderr, "%d x", s.nz+2);
	int rfmin, rfmax;
	rfmin = rfmax = marr[0].num;
	for(int i=1; i<=s.nz; ++i)
	{
		if(rfmin > marr[i].num) rfmin = marr[i].num;
		if(rfmax < marr[i].num) rfmax = marr[i].num;
	}
	if(rfmin > marr[s.nz+1].num) rfmin = marr[s.nz+1].num;
	if(rfmax < marr[s.nz+1].num) rfmax = marr[s.nz+1].num;
	if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin);
	else fprintf(stderr, " (%d - %d)) ", rfmin, rfmax);

	Cmatrix *dmat;
	Ivector filld(s.nz+2);
	Ivector idd(s.nz+2);

	dmat = new Cmatrix[s.nz+2];
	filld.init(0);

	for(int j=0; j<=s.nz+1; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o].equal(marr[j]) == 1 && filld(o) == 1)
			{
				idd(j) = o;
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			dmat[j] = orthomat_inv(marr[j]);
			fprintf(stderr, ".");
			idd(j) = j;
			filld(j) = 1;
		}
	}
	filld.strip();

	Cmatrix m_m, m_p;
	Cmatrix *mmff, *mmfb, *mmbf, *mmbb;
	Cmatrix *mpff, *mpfb, *mpbf, *mpbb;
	Ivector fill(s.nz+1);
	Ivector id(s.nz+1);

	mmff = new Cmatrix[s.nz+1];
	mmfb = new Cmatrix[s.nz+1];
	mmbf = new Cmatrix[s.nz+1];
	mmbb = new Cmatrix[s.nz+1];
	mpff = new Cmatrix[s.nz+1];
	mpfb = new Cmatrix[s.nz+1];
	mpbf = new Cmatrix[s.nz+1];
	mpbb = new Cmatrix[s.nz+1];
	fill.init(0);

	for(int j=0; j<=s.nz; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o  ].equal(marr[j  ]) == 1 &&
			   marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1)
			{
				id(j) = o+1;
				done = 1;
			}
			if(marr[o  ].equal(marr[j+1]) == 1 &&
			   marr[o+1].equal(marr[j  ]) == 1 && fill(o) == 1)
			{
				id(j) = -1*(o+1);
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			int nl, nr;
			nl = marr[j  ].num;
			nr = marr[j+1].num;
			transfermats(marr[j], marr[j+1], 
			             dmat[idd(j)], dmat[idd(j+1)], m_m, m_p);
			fprintf(stderr, ".");
			mmff[j] = m_m.submatrix(nl, nr, 0,  0 );
			mmfb[j] = m_m.submatrix(nl, nr, 0,  nr);
			mmbf[j] = m_m.submatrix(nl, nr, nl, 0 );
			mmbb[j] = m_m.submatrix(nl, nr, nl, nr);
			m_m.strip();
			mpff[j] = m_p.submatrix(nr, nl, 0,  0 );
			mpfb[j] = m_p.submatrix(nr, nl, 0,  nl);
			mpbf[j] = m_p.submatrix(nr, nl, nr, 0 );
			mpbb[j] = m_p.submatrix(nr, nl, nr, nl);
			m_p.strip();
			id(j) = j+1;
			fill(j) = 1;
		}
	}
	fill.strip();
	idd.strip();
	delete[] dmat;

	fprintf(stderr, ":");

	Cvector *t = NULL;
	t = new Cvector[s.nz];
	for(int i=0; i<=s.nz-1; ++i)
	{
		ModeArray ma = marr[i+1];
		int ni = ma.num;
		double l = s.hz(i+1)-s.hz(i);
		Cvector tmpt(ni);
		for(int j=0; j<=ni-1; ++j)
		{
/*
			if(ma(j).ppt == PROPAG)
			   tmpt(j) = expi(-l*ma(j).beta);
			else
			   tmpt(j) = exp(-l*ma(j).beta);
*/
			tmpt(j) = ma(j).efac(FORW, l);
		}
		t[i] = tmpt;
	}
	fprintf(stderr, ".");

	Cmatrix u, ri;
	Cmatrix *jff, *jfb, *jbf, *jbb;
	jff = new Cmatrix[s.nz+1];
	jfb = new Cmatrix[s.nz+1];
	jbf = new Cmatrix[s.nz+1];
	jbb = new Cmatrix[s.nz+1];

	for(int j=0; j<=s.nz; ++j)
	{
		u.unit(marr[j+1].num);
		if(id(j) > 0) 
		      ri = sub(u, mult(mpfb[ id(j)-1], mmbf[ id(j)-1]));
		else    
		      ri = sub(u, mult(mmfb[-id(j)-1], mpbf[-id(j)-1]));
		ri.inverse();
		if(id(j) > 0) 
		      jff[j] = mult(ri, mpff[ id(j)-1]);
		else          
		      jff[j] = mult(ri, mmff[-id(j)-1]);
		if(id(j) > 0) 
		      jfb[j] = mult(ri, mult(mpfb[ id(j)-1], mmbb[ id(j)-1]));
		else
		      jfb[j] = mult(ri, mult(mmfb[-id(j)-1], mpbb[-id(j)-1]));

		u.unit(marr[j].num);
		if(id(j) > 0) 
		      ri = sub(u, mult(mmbf[ id(j)-1], mpfb[ id(j)-1]));
		else    
		      ri = sub(u, mult(mpbf[-id(j)-1], mmfb[-id(j)-1]));
		ri.inverse();
		if(id(j) > 0) 
		      jbf[j] = mult(ri, mult(mmbf[ id(j)-1], mpff[ id(j)-1]));
		else
		      jbf[j] = mult(ri, mult(mpbf[-id(j)-1], mmff[-id(j)-1]));
		if(id(j) > 0) 
		      jbb[j] = mult(ri, mmbb[ id(j)-1]);
		else          
		      jbb[j] = mult(ri, mpbb[-id(j)-1]);
		fprintf(stderr, ".");
	}

	delete[] mmff;
	delete[] mmfb;
	delete[] mmbf;
	delete[] mmbb;
	delete[] mpff;
	delete[] mpfb;
	delete[] mpbf;
	delete[] mpbb;

	Cmatrix *kkff, *kkfb, *kkbf, *kkbb;
	kkff = new Cmatrix[s.nz+1];
	kkfb = new Cmatrix[s.nz+1];
	kkbf = new Cmatrix[s.nz+1];
	kkbb = new Cmatrix[s.nz+1];

	Cmatrix *aff, *afb, *abf, *abb;
	aff = new Cmatrix[s.nz];
	afb = new Cmatrix[s.nz];
	abf = new Cmatrix[s.nz];
	abb = new Cmatrix[s.nz];

	kkff[0] = jff[0];
	kkfb[0] = jfb[0];
	kkbf[0] = jbf[0];
	kkbb[0] = jbb[0];

	for(int j=0; j<=s.nz-1; ++j)
	{
		u.unit(marr[j+1].num);
		ri = sub(u, mult(diagmult(kkfb[j  ], t[j+1-1]), 
		                 diagmult(jbf[j+1], t[j+1-1])));
		ri.inverse();
		aff[j] = mult(ri, kkff[j]);
		afb[j] = mult(ri, mult(kkfb[j], diagmult(t[j+1-1], jbb[j+1])));
		ri = sub(u, mult(diagmult(jbf[j+1], t[j+1-1]), 
		                 diagmult(kkfb[j  ], t[j+1-1])));
		ri.inverse();
		abf[j] = mult(ri, mult(jbf[j+1], diagmult(t[j+1-1], kkff[j])));
		abb[j] = mult(ri, jbb[j+1]);

		ri = diagmult(jff[j+1], t[j+1-1]);
		kkff[j+1] = mult(ri, aff[j]);
		kkfb[j+1] = add(mult(ri, afb[j]), jfb[j+1]);
		ri = diagmult(kkbb[j], t[j+1-1]);
		kkbf[j+1] = add(mult(ri, abf[j]), kkbf[j]);
		kkbb[j+1] = mult(ri, abb[j]);
		fprintf(stderr, ".");
	}

	kff = kkff[s.nz];
	kfb = kkfb[s.nz];
	kbf = kkbf[s.nz];
	kbb = kkbb[s.nz];

	delete[] kkff;
	delete[] kkfb;
	delete[] kkbf;
	delete[] kkbb;

	delete[] jff;
	delete[] jfb;
	delete[] jbf;
	delete[] jbb;

	delete[] aff;
	delete[] afb;
	delete[] abf;
	delete[] abb;

	delete[] t;

	fprintf(stderr, " Ok.\n");
	return;
}

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

/* mode expansion solver, for QUEP implementation
   given the structure with the accompanying mode expansion,
   initialize the BepField (amplitude arrays),
   setup - transfer matrices for the input and output junction
         - a matrix that connects amplitudes on the first and last inner segment
         - matrices that determine the amplitudes on all inner segments */
void BepField::quepsetup(Cmatrix& m0mbf, Cmatrix& m0mbb,
                         Cmatrix& m0pff, Cmatrix& m0pfb,
                         Cvector& t1, Cmatrix& d0, Cmatrix& d1,
                         Cmatrix& mNmbf, Cmatrix& mNmbb,
                         Cmatrix& mNpff, Cmatrix& mNpfb,
                         Cvector& tN, Cmatrix& dN, Cmatrix& dNp,
                         Cmatrix& kff, Cmatrix& kfb, 
                         Cmatrix& kbf, Cmatrix& kbb,
                         Cmatrix *iff, Cmatrix *ifb, 
                         Cmatrix *ibf, Cmatrix *ibb,
                         Cvector *tx)
{
	fprintf(stderr, "BEP (");
	if(s.nz < 3)          bepfielderror("quepsetup: nz < 3"); 
	if(s.nz+2 != nums) bepfielderror("quepsetup: structure mismatch"); 
	fprintf(stderr, "%d x", s.nz+2);
	int rfmin, rfmax;
	if(afr[0].nel != marr[0].num) 
		bepfielderror("quepsetup: dimensions 0"); 
	abr[0] = Cvector(marr[0].num);
	afl[0] = Cvector(marr[0].num);
	abl[0] = Cvector(marr[0].num);
	rfmin = rfmax = marr[0].num;
	for(int i=1; i<=s.nz; ++i)
	{
		afl[i] = Cvector(marr[i].num);
		afr[i] = Cvector(marr[i].num);
		abl[i] = Cvector(marr[i].num);
		abr[i] = Cvector(marr[i].num);
		if(rfmin > marr[i].num) rfmin = marr[i].num;
		if(rfmax < marr[i].num) rfmax = marr[i].num;
	}
	if(abl[s.nz+1].nel != marr[s.nz+1].num) 
		bepfielderror("quepsetup: dimensions nz+1"); 
	afl[s.nz+1] = Cvector(marr[s.nz+1].num);
	afr[s.nz+1] = Cvector(marr[s.nz+1].num);
	abr[s.nz+1] = Cvector(marr[s.nz+1].num);
	if(rfmin > marr[s.nz+1].num) rfmin = marr[s.nz+1].num;
	if(rfmax < marr[s.nz+1].num) rfmax = marr[s.nz+1].num;
	if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin);
	else fprintf(stderr, " (%d - %d)) ", rfmin, rfmax);

	Cmatrix *dmat;
	Ivector filld(s.nz+2);
	Ivector idd(s.nz+2);

	dmat = new Cmatrix[s.nz+2];
	filld.init(0);

	for(int j=0; j<=s.nz+1; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o].equal(marr[j]) == 1 && filld(o) == 1)
			{
				idd(j) = o;
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			dmat[j] = orthomat_inv(marr[j]);
			fprintf(stderr, ".");
			idd(j) = j;
			filld(j) = 1;
		}
	}
	filld.strip();

	d0  = dmat[idd(0)];
	d1  = dmat[idd(1)];
	dN  = dmat[idd(s.nz)];
	dNp = dmat[idd(s.nz+1)];

	Cmatrix m_m, m_p;
	Cmatrix *mmff, *mmfb, *mmbf, *mmbb;
	Cmatrix *mpff, *mpfb, *mpbf, *mpbb;
	Ivector fill(s.nz+1);
	Ivector id(s.nz+1);

	mmff = new Cmatrix[s.nz+1];
	mmfb = new Cmatrix[s.nz+1];
	mmbf = new Cmatrix[s.nz+1];
	mmbb = new Cmatrix[s.nz+1];
	mpff = new Cmatrix[s.nz+1];
	mpfb = new Cmatrix[s.nz+1];
	mpbf = new Cmatrix[s.nz+1];
	mpbb = new Cmatrix[s.nz+1];
	fill.init(0);

	for(int j=0; j<=s.nz; ++j)
	{
		int done = 0;
		int o = 0;
		while(done == 0 && o<=j-1)
		{
			if(marr[o  ].equal(marr[j  ]) == 1 &&
			   marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1)
			{
				id(j) = o+1;
				done = 1;
			}
			if(marr[o  ].equal(marr[j+1]) == 1 &&
			   marr[o+1].equal(marr[j  ]) == 1 && fill(o) == 1)
			{
				id(j) = -1*(o+1);
				done = 1;
			}
			++o;
		}
		if(done == 0)
		{
			int nl, nr;
			nl = marr[j  ].num;
			nr = marr[j+1].num;
			transfermats(marr[j], marr[j+1], 
			             dmat[idd(j)], dmat[idd(j+1)], m_m, m_p);
			fprintf(stderr, ".");
			mmff[j] = m_m.submatrix(nl, nr, 0,  0 );
			mmfb[j] = m_m.submatrix(nl, nr, 0,  nr);
			mmbf[j] = m_m.submatrix(nl, nr, nl, 0 );
			mmbb[j] = m_m.submatrix(nl, nr, nl, nr);
			m_m.strip();
			mpff[j] = m_p.submatrix(nr, nl, 0,  0 );
			mpfb[j] = m_p.submatrix(nr, nl, 0,  nl);
			mpbf[j] = m_p.submatrix(nr, nl, nr, 0 );
			mpbb[j] = m_p.submatrix(nr, nl, nr, nl);
			m_p.strip();
			id(j) = j+1;
			fill(j) = 1;
		}
	}
	fill.strip();
	idd.strip();
	delete[] dmat;

	if(id(0) > 0) 
	{
		m0mbf = mmbf[id(0)-1];
		m0mbb = mmbb[id(0)-1];
		m0pff = mpff[id(0)-1];
		m0pfb = mpfb[id(0)-1];
	}
	else    
	{
		m0mbf = mpbf[-id(0)-1];
		m0mbb = mpbb[-id(0)-1];
		m0pff = mmff[-id(0)-1];
		m0pfb = mmfb[-id(0)-1];
	}
	if(id(s.nz) > 0) 
	{
		mNmbf = mmbf[id(s.nz)-1];
		mNmbb = mmbb[id(s.nz)-1];
		mNpff = mpff[id(s.nz)-1];
		mNpfb = mpfb[id(s.nz)-1];
	}
	else    
	{
		mNmbf = mpbf[-id(s.nz)-1];
		mNmbb = mpbb[-id(s.nz)-1];
		mNpff = mmff[-id(s.nz)-1];
		mNpfb = mmfb[-id(s.nz)-1];
	}
	fprintf(stderr, ":");

	Cvector *t = NULL;
	t = new Cvector[s.nz];
	for(int i=0; i<=s.nz-1; ++i)
	{
		ModeArray ma = marr[i+1];
		int ni = ma.num;
		double l = s.hz(i+1)-s.hz(i);
		Cvector tmpt(ni);
		for(int j=0; j<=ni-1; ++j)
		{
/*
			if(ma(j).ppt == PROPAG)
			   tmpt(j) = expi(-l*ma(j).beta);
			else
			   tmpt(j) = exp(-l*ma(j).beta);
*/
			tmpt(j) = ma(j).efac(FORW, l);
		}
		t[i] = tmpt;
	}
	fprintf(stderr, ".");

	t1 = t[0];
	for(int i=1; i<=s.nz-2; ++i) tx[i+1] = t[i];
	tN = t[s.nz-1];

	Cmatrix u, ri;
	Cmatrix *jff, *jfb, *jbf, *jbb;
	jff = new Cmatrix[s.nz];
	jfb = new Cmatrix[s.nz];
	jbf = new Cmatrix[s.nz];
	jbb = new Cmatrix[s.nz];

	for(int j=1; j<=s.nz-1; ++j)
	{
		u.unit(marr[j+1].num);
		if(id(j) > 0) 
		      ri = sub(u, mult(mpfb[ id(j)-1], mmbf[ id(j)-1]));
		else    
		      ri = sub(u, mult(mmfb[-id(j)-1], mpbf[-id(j)-1]));
		ri.inverse();
		if(id(j) > 0) 
		      jff[j] = mult(ri, mpff[ id(j)-1]);
		else          
		      jff[j] = mult(ri, mmff[-id(j)-1]);
		if(id(j) > 0) 
		      jfb[j] = mult(ri, mult(mpfb[ id(j)-1], mmbb[ id(j)-1]));
		else
		      jfb[j] = mult(ri, mult(mmfb[-id(j)-1], mpbb[-id(j)-1]));

		u.unit(marr[j].num);
		if(id(j) > 0) 
		      ri = sub(u, mult(mmbf[ id(j)-1], mpfb[ id(j)-1]));
		else    
		      ri = sub(u, mult(mpbf[-id(j)-1], mmfb[-id(j)-1]));
		ri.inverse();
		if(id(j) > 0) 
		      jbf[j] = mult(ri, mult(mmbf[ id(j)-1], mpff[ id(j)-1]));
		else
		      jbf[j] = mult(ri, mult(mpbf[-id(j)-1], mmff[-id(j)-1]));
		if(id(j) > 0) 
		      jbb[j] = mult(ri, mmbb[ id(j)-1]);
		else          
		      jbb[j] = mult(ri, mpbb[-id(j)-1]);
		fprintf(stderr, ".");
	}

	delete[] mmff;
	delete[] mmfb;
	delete[] mmbf;
	delete[] mmbb;
	delete[] mpff;
	delete[] mpfb;
	delete[] mpbf;
	delete[] mpbb;

	Cmatrix *kkff, *kkfb, *kkbf, *kkbb;
	kkff = new Cmatrix[s.nz];
	kkfb = new Cmatrix[s.nz];
	kkbf = new Cmatrix[s.nz];
	kkbb = new Cmatrix[s.nz];

	Cmatrix *aff, *afb, *abf, *abb;
	aff = new Cmatrix[s.nz-1];
	afb = new Cmatrix[s.nz-1];
	abf = new Cmatrix[s.nz-1];
	abb = new Cmatrix[s.nz-1];

	kkff[1] = jff[1];
	kkfb[1] = jfb[1];
	kkbf[1] = jbf[1];
	kkbb[1] = jbb[1];

	for(int j=1; j<=s.nz-2; ++j)
	{
		u.unit(marr[j+1].num);
		ri = sub(u, mult(diagmult(kkfb[j  ], t[j+1-1]), 
		                 diagmult(jbf[j+1], t[j+1-1])));
		ri.inverse();
		aff[j] = mult(ri, kkff[j]);
		afb[j] = mult(ri, mult(kkfb[j], diagmult(t[j+1-1], jbb[j+1])));
		ri = sub(u, mult(diagmult(jbf[j+1], t[j+1-1]), 
		                 diagmult(kkfb[j  ], t[j+1-1])));
		ri.inverse();
		abf[j] = mult(ri, mult(jbf[j+1], diagmult(t[j+1-1], kkff[j])));
		abb[j] = mult(ri, jbb[j+1]);

		ri = diagmult(jff[j+1], t[j+1-1]);
		kkff[j+1] = mult(ri, aff[j]);
		kkfb[j+1] = add(mult(ri, afb[j]), jfb[j+1]);
		ri = diagmult(kkbb[j], t[j+1-1]);
		kkbf[j+1] = add(mult(ri, abf[j]), kkbf[j]);
		kkbb[j+1] = mult(ri, abb[j]);
		fprintf(stderr, ".");
	}

	kff = kkff[s.nz-1];
	kfb = kkfb[s.nz-1];
	kbf = kkbf[s.nz-1];
	kbb = kkbb[s.nz-1];

	delete[] kkff;
	delete[] kkfb;
	delete[] kkbf;
	delete[] kkbb;

	delete[] jff;
	delete[] jfb;
	delete[] jbf;
	delete[] jbb;

	iff[s.nz-1] = aff[s.nz-2];
	ifb[s.nz-1] = afb[s.nz-2];
	ibf[s.nz-1] = abf[s.nz-2];
	ibb[s.nz-1] = abb[s.nz-2];

	for(int j=s.nz-1; j>=3; --j)
	{
		ri = diagmult(afb[j-2], t[j-1]);
		iff[j-1] = add(aff[j-2], mult(ri, ibf[j]));
		ifb[j-1] = mult(ri, ibb[j]);
		ri = diagmult(abb[j-2], t[j-1]);
		ibf[j-1] = add(abf[j-2], mult(ri, ibf[j]));
		ibb[j-1] = mult(ri, ibb[j]);
	}

	delete[] aff;
	delete[] afb;
	delete[] abf;
	delete[] abb;

	delete[] t;

	fprintf(stderr, " Ok.\n");
	return;
}

/* for the QUEP implementation: overlaps of perpendiculary oriented mode 
   expansion fields */
Cmatrix BepField::crossoverlap(const ModeArray& mh, double zpos,
                               Cvector& t1, Cvector& tN, 
                               Cmatrix& kff, Cmatrix& kfb, 
                               Cmatrix& kbf, Cmatrix& kbb,
                               Cmatrix *iff, Cmatrix *ifb, 
                               Cmatrix *ibf, Cmatrix *ibb)
{

	Cmatrix g, vt;
	double p = zpos;

	Cmatrix vff, vfb, vbf, vbb;

	Cmatrix *gff, *gfb, *gbf, *gbb;	
	gff = new Cmatrix[s.nz+1];
	gfb = new Cmatrix[s.nz+1];
	gbf = new Cmatrix[s.nz+1];
	gbb = new Cmatrix[s.nz+1];

	for(int l=1; l<=s.nz; ++l)
	{
		mh.crossoverlap(marr[l], p, s.segment(l), 
		                gff[l], gfb[l], gbf[l], gbb[l]);
	}

	g = gff[1];
	vt = g;
	g = gfb[1];
	vt.addeq(mult(g, diagmult(kbf, t1)));
	for(int l=2; l<=s.nz-1; ++l)
	{
		g = gff[l];
		vt.addeq(mult(g, diagmult(iff[l], t1)));
		g = gfb[l]; 
		vt.addeq(mult(g, diagmult(ibf[l], t1)));
	}
	g = gff[s.nz]; 
	vt.addeq(mult(g, diagmult(kff, t1)));
	vff = vt;

	g = gbf[1];
	vt = g;
	g = gbb[1];
	vt.addeq(mult(g, diagmult(kbf, t1)));
	for(int l=2; l<=s.nz-1; ++l)
	{
		g = gbf[l];
		vt.addeq(mult(g, diagmult(iff[l], t1)));
		g = gbb[l];
		vt.addeq(mult(g, diagmult(ibf[l], t1)));
	}
	g = gbf[s.nz];
	vt.addeq(mult(g, diagmult(kff, t1)));
	vbf = vt;

	g = gfb[1];
	vt = mult(g, diagmult(kbb, tN));
	for(int l=2; l<=s.nz-1; ++l)
	{
		g = gff[l];
		vt.addeq(mult(g, diagmult(ifb[l], tN)));
		g = gfb[l];
		vt.addeq(mult(g, diagmult(ibb[l], tN)));
	}
	g = gff[s.nz];
	vt.addeq(mult(g, diagmult(kfb, tN)));
	g = gfb[s.nz];
	vt.addeq(g);
	vfb = vt;

	g = gbb[1];
	vt = mult(g, diagmult(kbb, tN));
	for(int l=2; l<=s.nz-1; ++l)
	{
		g = gbf[l];
		vt.addeq(mult(g, diagmult(ifb[l], tN)));
		g = gbb[l];
		vt.addeq(mult(g, diagmult(ibb[l], tN)));
	}
	g = gbf[s.nz];
	vt.addeq(mult(g, diagmult(kfb, tN)));
	g = gbb[s.nz];
	vt.addeq(g);
	vbb = vt;

	vt.strip();
	g.strip();

	for(int l=1; l<=s.nz; ++l)
	{
		gff[l].strip();
		gfb[l].strip();
		gbf[l].strip();
		gbb[l].strip();
	}
	delete[] gff;
	delete[] gfb;
	delete[] gbf;
	delete[] gbb;
	
	Cmatrix vmat(mh.num+mh.num, marr[1].num+marr[s.nz].num);
	vmat.setsubmatrix(vff, 0,      0     );
	vmat.setsubmatrix(vfb, 0,      marr[1].num);
	vmat.setsubmatrix(vbf, mh.num, 0     );
	vmat.setsubmatrix(vbb, mh.num, marr[1].num);    	
	vff.strip();
	vfb.strip();
	vbf.strip();
	vbb.strip();
	return vmat;
} 

/* ------------------------------------------------------------------- */
/* for the vBEP / vQUEP implementaion */

/* prepare for the vTE/vTM mode representation, for given ky */
void BepField::vectorify(Complex ky)
{
	for(int i=0; i<=s.nz+1; ++i) marr[i].vectorify(ky);
	return;
}

/* compute the discretized bidirectional mode spectra:
   vTE and vTM modes, uniform wavenumber ky, 
   2*Mx expansion terms per vertical slice on a vertical 
   computational window cw, 
   type = 0: boundary conditions LIMD for TE and LIMN for TM modes,
   type = 1: boundary conditions LIMN for TE and LIMD for TM modes */
void BepField::discretize(Complex ky, int type, Interval cw, int Mx, int quiet)
{
	Waveguide wg;
	ModeArray ma;
	if(quiet != 1) fprintf(stderr, "SpecDisc (vM):");
	for(int i=0; i<=s.nz+1; ++i)
	{	
		int done = 0; 
		int j=0;
		while(j<i && done==0)
		{
			if(s(i).equal(s(j))==1)
			{
				marr[i] = marr[j];
				done = 1;
			}
			++j;
		}
		if(done==0)
		{
			wg = s(i);
			switch(type)
			{
			case 0:
				modespectrum(wg, TE, cw, LIMD, Mx, ma, 1);
				marr[i] = ma;
				if(quiet != 1) fprintf(stderr, ".");
				modespectrum(wg, TM, cw, LIMN, Mx, ma, 1);
				marr[i].merge(ma);
				if(quiet != 1) fprintf(stderr, ".");
				break;
			case 1:
				modespectrum(wg, TE, cw, LIMN, Mx, ma, 1);
				marr[i] = ma;
				if(quiet != 1) fprintf(stderr, ".");
				modespectrum(wg, TM, cw, LIMD, Mx, ma, 1);
				marr[i].merge(ma);
				if(quiet != 1) fprintf(stderr, ".");
				break;	
			default: 
				bepfielderror("discretize, ky: type");
				break;
			}
			marr[i].vectorify(ky);
			if(quiet != 1) fprintf(stderr, ",");
		}
	}
	if(quiet != 1) fprintf(stderr, " Ok.\n");
	Cvector zero;
	zero = Cvector(marr[0].num); zero.init(CC0);
	setf(0, SBRIGHT, zero);
	zero = Cvector(marr[s.nz+1].num); zero.init(CC0);
	setb(s.nz+1, SBLEFT, zero);
	check();
	return;
}

/* full vBEP initialization: setup the container for the bidirectional 
   mode expansion, compute the spectral discretization 
   rc: the structure under investigation, 
   ky: fields of vTE and vTM polarization, with y-wavenumber ky
   wx: vertical computational window,
   Mx: number of expansion terms per vertical slice,
   suppress log output, if quiet == 1 */
BepField::BepField(Circuit rc, Complex ky, Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	if(quiet != 1)
	{
fprintf(stderr, "------------ vBEP -------------------------------------------------- '14 ---\n");
fprintf(stderr, "vTE, vTM, ky = (%.10g + i %.10g) k0\n",
                ky.re/val_k0(rc.lambda), ky.im/val_k0(rc.lambda));
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d\n", wx.x0, wx.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, "lambda: %.10g  k0: %g", rc.lambda, val_k0(rc.lambda));
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	(*this) = BepField(rc, 1);
	discretize(ky, 0, wx, Mx, quiet);
}
BepField::BepField(Circuit rc, Complex ky, Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc, ky, wx, Mx, 0);
}
BepField::BepField(SegWgStruct rc, Complex ky, 
                   Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), ky, wx, Mx, quiet);
}
BepField::BepField(SegWgStruct rc, Complex ky, 
                   Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), ky, wx, Mx, 0);
}
/* ... with ky specified by the input field: 
   guided mode of polarization pol, order ord, 
   coming in at border bdr, with amplitude amp, 
   at angle of incidence theta (^o) */
BepField::BepField(Circuit rc, SBorder bdr,
                   Polarization pol, int ord, Complex amp, double theta,
                   Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc, 1);
	Waveguide ewg;
	switch(bdr) 
	{
	case SBLEFT:
		ewg = s(0);
		break;
	case SBRIGHT:
		ewg = s(s.nz+1);
		break;
	}
	ModeArray ma;
	modeanalysis(ewg, pol, ma, 1);
	if(ord < 0 || ord >= ma.num) bepfielderror("vBEP, setup, mode order");
	double N = ma(ord).neff;
	double b = ma(ord).beta;
	if(theta <= -90.0 || theta >= 90.0) bepfielderror("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; amp = %g+i%g; angle = %g^o\n", 
                ma(ord).ids, N, b, amp.re, amp.im, theta);
	}
	(*this) = BepField(rc, ky, wx, Mx, quiet);
	if(pol == TE) input(bdr, amp, ord);
	else          input(bdr, amp, Mx+ord);
}
BepField::BepField(Circuit rc, SBorder bdr,
                   Polarization pol, int ord, Complex amp, double theta,
                   Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc, bdr, pol, ord, amp, theta, wx, Mx, 0);
}
BepField::BepField(SegWgStruct rc, SBorder bdr,
                   Polarization pol, int ord, Complex amp, double theta,
                   Interval wx, int Mx, int quiet)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), 
	                   bdr, pol, ord, amp, theta, wx, Mx, quiet);
}
BepField::BepField(SegWgStruct rc, SBorder bdr,
                   Polarization pol, int ord, Complex amp, double theta,
                   Interval wx, int Mx)
{
	marr = NULL; afl = afr = abl = abr = NULL;
	(*this) = BepField(rc.circuit(), 
	                   bdr, pol, ord, amp, theta, wx, Mx, 0);
}


/* ------------------------------------------------------------------- */
/* preparation of physical field plots:
   selection of the global phase / snapshot time 0 within one period,
   a uniform polarization is assumed */

/* multiply the entire field by a phase factor exp(i phi),
   the internal representation is modified accordingly */
void BepField::adjustphase(double phi)
{
	Complex e = expi(phi);
	for(int j=0; j<=nums-1; ++ j)
	{
		afl[j].multeq(e);
		afr[j].multeq(e);
		abl[j].multeq(e);
		abr[j].multeq(e);
	}
	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 BepField::adjustphase(double x, double z)
{
	int seg = s.segidx(z);
	if(marr[seg].num <= 0) return;
	Polarization p = (marr[seg])(0).pol;
	Fcomp cp = principalcomp(p);
	adjustphase(-field(cp, x, z).arg());
	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 BepField::adjustphase(int numx, int numz)
{
	if(marr[0].num <= 0) return;
	Polarization p = (marr[0])(0).pol;
	Fcomp cp = principalcomp(p);
	Interval wx, wz;
	s.defaultwindow(wx, wz);
	double x, z;
	double xs, zs;
	xs = wx.len()/((double) numx);
	zs = wz.len()/((double) numz);
	x = wx.x0;
	double xm = 0.0, zm = 0.0, maxf = 0.0;
	double f;
	for(int xi=0; xi <=numx; ++xi)
	{
		z = wz.x0;
		for(int zi=0; zi <=numz; ++zi)
		{
			f = field(cp, x, z).sqabs();
			if(f > maxf) { xm=x; zm=z; maxf=f; }
			z += zs;
		}
		x += xs;
	}
	adjustphase(-field(cp, xm, zm).arg());
	return;
}

/* ... as before, with (coarse) defaults for numx, numz */
void BepField::adjustphase()
{
	Interval wx, wz;
	s.defaultwindow(wx, wz);
	int numx, numz;	
	numx = (int) ceil(sqrt(2500.0*wx.len()/wz.len()));
	numz = numx*((int) ceil(wz.len()/wx.len()));
	adjustphase(numx, numz);
	return;
}

/* ------------------------------------------------------------------- */
/* visualization of the light propagation, output: MATLAB scripts */
	
/* mode interference plots, image of component cp 
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   foa: MOD, SQR, REP, IMP, ORG = REP
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void BepField::plot(Fcomp cp, Afo foa, 
                    double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.plot(cp, foa, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	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 BepField::fplot(Fcomp cp, 
                     double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.fplot(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	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 BepField::phasemap(Fcomp cp,
                        double xbeg, double xend, double zbeg, double zend,
                        int npx, int npz, char ext0, char ext1) const 
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.phasemap(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* image plot of the field "strength", id: one of mE, mH, qE, qH, mW  
   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 BepField::plot(FSid id,
                    double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.plot(id, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* electromagnetic energy density image
   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 BepField::plot(double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.plot(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	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 BepField::movie(double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, int ntfr, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.movie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1);
	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 BepField::fmovie(double xbeg, double xend, double zbeg, double zend,
                      int npx, int npz, int ntfr, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.fmovie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1);
	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 BepField::viewer(double xbeg, double xend, double zbeg, double zend,
                      int npx, int npz, char ext0, char ext1) const
{
	QuepField qfld;
	qfld.h = (*this);
	qfld.QF_pol    = BF_pol;
	qfld.QF_lambda = BF_lambda;
	qfld.QF_vM     = BF_vM;
	qfld.QF_ky     = BF_ky;
	qfld.viewer(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

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

/* orthogonality matrix for a mode superposition, inverted */
Cmatrix orthomat_inv(ModeArray& ma)
{
	Cmatrix dff, dfb, dbf, dbb;
	Cmatrix d;
	int n = ma.num;
	dff = Cmatrix(n, n);
	dfb = Cmatrix(n, n);
	dbf = Cmatrix(n, n);
	dbb = Cmatrix(n, n);
	for(int l=0; l<=n-1; ++l)
	{
		for(int m=0; m<=n-1; ++m)
		{
			bidir_overlaps(ma(l), ma(m), 
			               dff(l, m), dfb(l, m), 
			               dbf(l, m), dbb(l, m));
		}
	}
	d = Cmatrix(n+n, n+n);
	d.setsubmatrix(dff, 0, 0);
	d.setsubmatrix(dfb, 0, n);
	d.setsubmatrix(dbf, n, 0);
	d.setsubmatrix(dbb, n, n);
	dff.strip();
	dfb.strip();
	dbf.strip();
	dbb.strip();

	d.inverse();
	return d;
}


/* overlap matrix for two mode superpositions */
Cmatrix overlapmat(ModeArray& m2, ModeArray& m1)
{    
	Cmatrix off, ofb, obf, obb;
	Cmatrix o;

	off = Cmatrix(m2.num, m1.num);
	ofb = Cmatrix(m2.num, m1.num);
	obf = Cmatrix(m2.num, m1.num);
	obb = Cmatrix(m2.num, m1.num);
	for(int l=0; l<=m2.num-1; ++l)
	{
		for(int m=0; m<=m1.num-1; ++m)
		{
			bidir_overlaps(m2(l), m1(m), 
			               off(l, m), ofb(l, m), 
			               obf(l, m), obb(l, m));
		}
	}
	o = Cmatrix(m2.num+m2.num, m1.num+m1.num);
	o.setsubmatrix(off, 0,      0     );
	o.setsubmatrix(ofb, 0,      m1.num);
	o.setsubmatrix(obf, m2.num, 0     );
	o.setsubmatrix(obb, m2.num, m1.num);
	off.strip();
	ofb.strip();
	obf.strip();
	obb.strip();
	return o;
}

/* compute transfer matrices for a junction at z=0:
   m1: local modes for z<0 
   m2: local modes for z>0
   mat1, mat2: output, overlap integrals between forward and backward 
       propagating modes on both sides of the junction */
void transfermats(ModeArray& m1, ModeArray& m2, Cmatrix& mat1, Cmatrix& mat2)
{
	if(m1.num <= 0 || m2.num <= 0) return;
	Cmatrix o;
	Cmatrix d;
	o = overlapmat(m2, m1);
	d = orthomat_inv(m2);
	mat2 = mult(d, o);
	d.strip();
	d = orthomat_inv(m1);
	o.adjoint();
	mat1 = mult(d, o);
	d.strip();
	o.strip();

	return;
}

/* compute transfer matrices for a junction at z=0:
   m1: local modes for z<0 
   m2: local modes for z>0
   d1, d2: given orthogonality matrices for mode sets m1, m2
   mat1, mat2: output, overlap integrals between forward and backward 
       propagating modes on both sides of the junction */
void transfermats(ModeArray& m1, ModeArray& m2, 
                  Cmatrix& d1, Cmatrix& d2, Cmatrix& mat1, Cmatrix& mat2)
{
	if(m1.num <= 0 || m2.num <= 0) return;
	Cmatrix o;
	o = overlapmat(m2, m1);
	mat2 = mult(d2, o);
	o.adjoint();
	mat1 = mult(d1, o);
	o.strip();
	return;
}

/* compute transfer matrices for a junction at z=0:
   m1: local modes for z<0 
   m2: local modes for z>0
   m1ff - m2bb: output, overlap integrals between forward and backward 
       propagating modes on both sides of the junction */
void transfermats(ModeArray& m1, ModeArray& m2, 
                  Cmatrix& m1ff, Cmatrix& m1fb, Cmatrix& m1bf, Cmatrix& m1bb,
                  Cmatrix& m2ff, Cmatrix& m2fb, Cmatrix& m2bf, Cmatrix& m2bb)
{
	if(m1.num <= 0 || m2.num <= 0) return;

	Cmatrix mat1, mat2;
	transfermats(m1, m2, mat1, mat2);

	m2ff = mat2.submatrix(m2.num, m1.num, 0,      0     );
	m2fb = mat2.submatrix(m2.num, m1.num, 0,      m1.num);
	m2bf = mat2.submatrix(m2.num, m1.num, m2.num, 0     );
	m2bb = mat2.submatrix(m2.num, m1.num, m2.num, m1.num);
	mat2.strip();
	
	m1ff = mat1.submatrix(m1.num, m2.num, 0,      0     );
	m1fb = mat1.submatrix(m1.num, m2.num, 0,      m2.num);
	m1bf = mat1.submatrix(m1.num, m2.num, m1.num, 0     );
	m1bb = mat1.submatrix(m1.num, m2.num, m1.num, m2.num);
	mat1.strip();
	return;
}

/* - specific driver routines --------------------------------------------- */

/* a long three segment coupler:
   m1: local modes for z<0 
   m2: local modes for 0<z<l, modes are assumed to be ordered by (squared)
       propagation constants, fundamental guided mode first
   m3: local modes for l<z
   amp1_f (input):  amplitudes of the forward propagating modes at z=0-, 
   amp1_b (output): amplitudes of the backward propagating modes at z=0-,
                    mode reflectivities   
   amp3_f (output): amplitudes of the forward propagating modes at z=L+,
                    mode transmittivities 
   amp2_f_a, amp2_b_a (output): amplitudes of forward and backwards traveling 
                    modes at z=0+
   amp2_f_b, amp2_b_b (output): amplitudes of forward and backwards traveling 
                    modes at z=L- */
void long_3segment_coupler(ModeArray& m1, 
                           ModeArray& m2, 
                           ModeArray& m3, 
                           Cvector fpert, Cvector bpert,
	                   Cvector amp1_f, double l, 
                           Cvector& amp1_b, Cvector& amp3_f, 
                           Cvector& amp2_f_a, Cvector& amp2_b_a,
                           Cvector& amp2_f_b, Cvector& amp2_b_b)
{
	fprintf(stderr, "Long three segment coupler:\n");
	Cmatrix m_1, m_2, q_2, q_3;
	transfermats(m1, m2, m_1, m_2);
	transfermats(m2, m3, q_2, q_3);

	int ng, nr;
	ng = 0;	
	while(m2(ng).ppt == PROPAG && m2(ng).npcB > 0.0) ++ng;
	nr = m2.num-ng;
	fprintf(stderr, "Coupling segment: %d guided modes,\n", ng);
	fprintf(stderr, "                  %d modes in the radiation field.\n", 
		nr);
	
	Cmatrix m1ffg, m1ffr, m1fbg, m1fbr,
	        m1bfg, m1bfr, m1bbg, m1bbr; 

	fprintf(stderr, "    >>> ");
	m1ffg = m_1.submatrix(m1.num, ng, 0,      0        );
	m1ffr = m_1.submatrix(m1.num, nr, 0,      ng       );
	m1fbg = m_1.submatrix(m1.num, ng, 0,      m2.num   );
	m1fbr = m_1.submatrix(m1.num, nr, 0,      m2.num+ng);
	m1bfg = m_1.submatrix(m1.num, ng, m1.num, 0        );
	m1bfr = m_1.submatrix(m1.num, nr, m1.num, ng       );
	m1bbg = m_1.submatrix(m1.num, ng, m1.num, m2.num   );
	m1bbr = m_1.submatrix(m1.num, nr, m1.num, m2.num+ng);
	m_1.strip();
	
	Cmatrix m2fgf, m2fgb,
	        m2frf, m2frb,
	        m2bgf, m2bgb,
                m2brf, m2brb;

	m2fgf = m_2.submatrix(ng, m1.num, 0,         0     );
	m2fgb = m_2.submatrix(ng, m1.num, 0,         m1.num);
	m2frf = m_2.submatrix(nr, m1.num, ng,        0     );
	m2frb = m_2.submatrix(nr, m1.num, ng,        m1.num);
	m2bgf = m_2.submatrix(ng, m1.num, m2.num,    0     );
	m2bgb = m_2.submatrix(ng, m1.num, m2.num,    m1.num);
	m2brf = m_2.submatrix(nr, m1.num, m2.num+ng, 0     );
	m2brb = m_2.submatrix(nr, m1.num, m2.num+ng, m1.num);
	m_2.strip();
	
	Cmatrix q2fgf, q2fgb,
	        q2frf, q2frb,
	        q2bgf, q2bgb,
	        q2brf, q2brb;

	q2fgf = q_2.submatrix(ng, m3.num, 0,         0     );
	q2fgb = q_2.submatrix(ng, m3.num, 0,         m3.num);
	q2frf = q_2.submatrix(nr, m3.num, ng,        0     );
	q2frb = q_2.submatrix(nr, m3.num, ng,        m3.num);
	q2bgf = q_2.submatrix(ng, m3.num, m2.num,    0     );
	q2bgb = q_2.submatrix(ng, m3.num, m2.num,    m3.num);
	q2brf = q_2.submatrix(nr, m3.num, m2.num+ng, 0     );
	q2brb = q_2.submatrix(nr, m3.num, m2.num+ng, m3.num);
	q_2.strip();
	
	Cmatrix q3ffg, q3ffr, q3fbg, q3fbr,
	        q3bfg, q3bfr, q3bbg, q3bbr;
	
	q3ffg = q_3.submatrix(m3.num, ng, 0,      0        );
	q3ffr = q_3.submatrix(m3.num, nr, 0,      ng       );
	q3fbg = q_3.submatrix(m3.num, ng, 0,      m2.num   );
	q3fbr = q_3.submatrix(m3.num, nr, 0,      m2.num+ng);
	q3bfg = q_3.submatrix(m3.num, ng, m3.num, 0        );
	q3bfr = q_3.submatrix(m3.num, nr, m3.num, ng       );
	q3bbg = q_3.submatrix(m3.num, ng, m3.num, m2.num   );
	q3bbr = q_3.submatrix(m3.num, nr, m3.num, m2.num+ng);
	q_3.strip();

	fprintf(stderr, ".");
	Cmatrix tfp(ng, ng);
	Cmatrix tbp(ng, ng);
	Cmatrix tfm(ng, ng);
	Cmatrix tbm(ng, ng);
	Complex e;

	tfp.init(CC0);
	tbp.init(CC0);
	tfm.init(CC0);
	tbm.init(CC0);
	for(int j=0; j<=ng-1; ++j)
	{
		e = CCI*(Complex(m2(j).beta, 0.0)+fpert(j))*l*(-1.0);
		tfp(j, j) = exp(e);
		e = CCI*(Complex(m2(j).beta, 0.0)+bpert(j))*l;
		tbp(j, j) = exp(e);
		e = CCI*(Complex(m2(j).beta, 0.0)+fpert(j))*l;
		tfm(j, j) = exp(e);
		e = CCI*(Complex(m2(j).beta, 0.0)+bpert(j))*l*(-1.0);
		tbm(j, j) = exp(e);
	}
	
	fprintf(stderr, ".");
	Cmatrix u;
	Cmatrix ra, rb, rc, rd, re, rf;

	u.unit(m1.num);
	ra = add(u, mult(-1.0, mult(m1bfr, m2frb)));
	rb = add(mult(m1bfg, mult(tfm, q2fgf)), mult(m1bbg, mult(tbm, q2bgf))); 
	rc = mult(m1bfr, m2frf);
	u.unit(m3.num);
	rd = add(u, mult(-1.0, mult(q3fbr, q2brf)));
	re = add(mult(q3ffg, mult(tfp, m2fgf)), mult(q3fbg, mult(tbp, m2bgf))); 
	rf = add(mult(q3ffg, mult(tfp, m2fgb)), mult(q3fbg, mult(tbp, m2bgb))); 

	fprintf(stderr, ".");
	ra.inverse();
	Cmatrix r1, r2, r;
	fprintf(stderr, ".");
	r1 = add(rd, mult(-1.0, mult(rf, mult(ra, rb))));
	fprintf(stderr, ".");
	r1.inverse();
	fprintf(stderr, ".");
	r2 = add(re, mult(rf, mult(ra, rc)));
	r = mult(r1, r2);
	r1.strip();
	r2.strip();

	fprintf(stderr, ".");
	Dvector p1(m1.num);	
	Dvector p2(m2.num);	
	Dvector p3(m3.num);	
	for(int i=0; i<=m1.num-1; ++i) p1(i) = sqrt(fabs(m1(i).power()));
	for(int i=0; i<=m2.num-1; ++i) p2(i) = sqrt(fabs(m2(i).power()));
	for(int i=0; i<=m3.num-1; ++i) p3(i) = sqrt(fabs(m3(i).power()));
	Cvector f1(m1.num);
	for(int i=0; i<=m1.num-1; ++i) f1(i) = amp1_f(i)/p1(i);
	
	fprintf(stderr, ".");
	amp3_f = mult(r, f1);
	amp1_b = add(mult(mult(ra, rb), amp3_f), mult(mult(ra, rc), f1));

	fprintf(stderr, ".");
	Cvector v2(m2.num);
	v2.init(CC0);
	v2.setsubvector(add(mult(m2fgf, amp1_f), mult(m2fgb, amp1_b)), 0);
	v2.setsubvector(add(mult(m2frf, amp1_f), mult(m2frb, amp1_b)), ng);
	amp2_f_a = v2;
	v2.init(CC0);
	v2.setsubvector(add(mult(m2bgf, amp1_f), mult(m2bgb, amp1_b)), 0);
	amp2_b_a = v2;
	v2.init(CC0);
	v2.setsubvector(mult(q2fgf, amp3_f), 0);
	amp2_f_b = v2;
	v2.init(CC0);
	v2.setsubvector(mult(q2bgf, amp3_f), 0);
	v2.setsubvector(mult(q2brf, amp3_f), ng);
	amp2_b_b = v2;

	fprintf(stderr, ".");
	for(int i=0; i<=m3.num-1; ++i) amp3_f(i) = amp3_f(i)*p3(i);
	for(int i=0; i<=m1.num-1; ++i) amp1_b(i) = amp1_b(i)*p1(i);
	for(int i=0; i<=m2.num-1; ++i) 
	{
		amp2_f_a(i) = amp2_f_a(i)*p2(i);
		amp2_b_a(i) = amp2_b_a(i)*p2(i);
		amp2_f_b(i) = amp2_f_b(i)*p2(i);
		amp2_b_b(i) = amp2_b_b(i)*p2(i);
	}
	fprintf(stderr, " <<<\n");
/*
Cvector t, q;
Cvector vf1, vb1;
Cvector vf2ag, vf2ar, vb2ag, vb2ar;
Cvector vf2bg, vf2br, vb2bg, vb2br;
Cvector vf3, vb3;

vf1 = amp1_f;
vb1 = amp1_b;
vf3 = amp3_f;
vb3 = Cvector(m3.num); vb3.init(CC0);
vf2ag = amp2_f_a.subvector(ng, 0);
vf2ar = amp2_f_a.subvector(nr, ng);
vb2ag = amp2_b_a.subvector(ng, 0);
vb2ar = amp2_b_a.subvector(nr, ng);
vf2bg = amp2_f_b.subvector(ng, 0);
vf2br = amp2_f_b.subvector(nr, ng);
vb2bg = amp2_b_b.subvector(ng, 0);
vb2br = amp2_b_b.subvector(nr, ng);

t = add( add( mult(m1bfg, vf2ag), mult(m1bfr, vf2ar)), mult(m1bbg, vb2ag)); 
q = add( vb1, mult(-1.0, t));
fprintf(stderr, "(1): %g\n", q.norm());

t = add( mult(m2fgf, vf1), mult(m2fgb, vb1)); 
q = add( vf2ag, mult(-1.0, t));
fprintf(stderr, "(2): %g\n", q.norm());

t = add( mult(m2frf, vf1), mult(m2frb, vb1)); 
q = add( vf2ar, mult(-1.0, t));
fprintf(stderr, "(3): %g\n", q.norm());

t = add( mult(m2bgf, vf1), mult(m2bgb, vb1)); 
q = add( vb2ag, mult(-1.0, t));
fprintf(stderr, "(4): %g\n", q.norm());

t = mult(q2fgf, vf3); 
q = add( vf2bg, mult(-1.0, t));
fprintf(stderr, "(5): %g\n", q.norm());

t = mult(q2bgf, vf3); 
q = add( vb2bg, mult(-1.0, t));
fprintf(stderr, "(6): %g\n", q.norm());

t = mult(q2brf, vf3); 
q = add( vb2br, mult(-1.0, t));
fprintf(stderr, "(7): %g\n", q.norm());

t = add( add( mult(q3ffg, vf2bg), mult(q3fbg, vb2bg)), mult(q3fbr, vb2br)); 
q = add( vf3, mult(-1.0, t));
fprintf(stderr, "(8): %g\n", q.norm());

t = add( add( mult(m1ffg, vf2ag), mult(m1ffr, vf2ar)), mult(m1fbg, vb2ag)); 
q = add( vf1, mult(-1.0, t));
fprintf(stderr, "(9): %g\n", q.norm());

t = add( add( mult(q3bfg, vf2bg), mult(q3bbg, vb2bg)), mult(q3bbr, vb2br)); 
q = t;
fprintf(stderr, "(10): %g\n", q.norm());

t = add( mult(m2brf, vf1), mult(m2brb, vb1)); 
q = t;
fprintf(stderr, "(11): %g\n", q.norm());

t = add( mult(q2frf, vf3), mult(q2frb, vb3)); 
q = t;
fprintf(stderr, "(12): %g\n", q.norm());
*/

	fprintf(stderr, "Ok.\n");

	return;
}

/* image plot of the light propagation across the long three segment coupler,
   output: MATLAB script 
   m1: local modes for z<0
   m2: local modes for 0<z<l
   m3: local modes for l<z
   amp1_f: amplitudes of the forward propagating modes at z=0- 
   amp1_b: amplitudes of the backward propagating modes at z=0-
   amp3_f: amplitudes of the forward propagating modes at z=l+
   amp2_f_a, amp2_b_a: amplitudes of forward and backwards traveling 
                       modes at z=0+
   amp2_f_b, amp2_b_b: amplitudes of forward and backwards traveling 
                       modes at z=l-
   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 l3sc_plot(ModeArray& m1, ModeArray& m2, ModeArray& m3, 
               Cvector& amp1_f, Cvector& amp1_b, Cvector amp3_f, 
               Cvector& amp2_f_a, Cvector& amp2_b_a, 
               Cvector& amp2_f_b, Cvector& amp2_b_b, 
               double l,  
               double xbeg, double xend, double zbeg, double zend,
               int npx, int npz, char ext0, char ext1) 
{
	FILE *dat;
	char name[13] = "l3sc__I.m";
	double x, dx, z, dz;
	double f; 
	int i, j;
	name[4] = ext0;
	name[5] = ext1;
	fprintf(stderr, "l3sc_plot(z=%g -> %g) >> %s\n", zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_l3scgeoxz_2(dat, m1(m1.num-1).wg, 
	                       m2(m2.num-1).wg, 
	                       m3(m3.num-1).wg,
                               l, xbeg, xend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);

	Dmatrix fld(npx, npz);
	dx = (xend-xbeg)/(npx-1);
	dz = (zend-zbeg)/(npz-1);

	Cvector amp3_b(amp3_f.nel);
	amp3_b.init(CC0);
	for(j=0; j<=npz-1; ++j)
	{
		z = zbeg+j*dz;
		for(i=0; i<=npx-1; ++i) 
		{
			x = xbeg+i*dx;
			if(z<=0.0) 
			   f = m1.field(amp1_f, amp1_b, SZ, x, z).abs();
			else
			{
			   if(z <= l)
			   {
//			     if(z <= l/2.0)
//			        f = m2.field(amp2_f_a, amp2_b_a, 
//			                     SZ, x, z).abs();
//			     else
//			        f = m2.field(amp2_f_b, amp2_b_b, 
//			                     SZ, x, z-l).abs();
			        f = (l-z)/l*m2.field(amp2_f_a, amp2_b_a, 
			                             SZ, x, z).abs()
			              + z/l*m2.field(amp2_f_b, amp2_b_b, 
			                             SZ, x, z-l).abs();
			   }
			   else
			     f = m3.field(amp3_f, amp3_b, SZ, x, z-l).abs();
			}
			fld(i, j) = sqrt(f);
		}
	}
	mlout_fld(dat, npx, npz, SZ, fld);
	name[7] = 0; 
	mlout_1Dl3scimage(name, dat); 
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}