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

/*
 * quepfld.cpp
 * Quadridirectional mode expansion, crossed 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"matlvis.h"    
#include"quepfld.h"    

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

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

/* initialize */
QuepField::QuepField()
{
	h = BepField();
	v = BepField();
	cwx = Interval();
	cwz = Interval();
}
QuepField::QuepField(SegWgStruct hst, Interval wx, Interval wz)
{
	if(hst.nz < 1) 
	quepfielderror("at least one inner segment (hst.nz >=1) required");
	SegWgStruct vst;
	vst = hst.rotate();

	if(wx.x0 >= vst.hz(0)) 
		quepfielderror("illegal computational window (1)");
	if(wx.x1 <= vst.hz(vst.nz)) 
		quepfielderror("illegal computational window (2)");
	if(wz.x0 >= hst.hz(0)) 
		quepfielderror("illegal computational window (3)");
	if(wz.x1 <= hst.hz(hst.nz)) 
		quepfielderror("illegal computational window (4)");
	cwx = wx;
	cwz = wz;
	SegWgStruct hste, vste;
	hste = hst.expand(cwz);
	vste = vst.expand(cwx);
	h = BepField(hste);
	v = BepField(vste); 
}
QuepField::QuepField(Circuit rc, Interval wx, Interval wz)
{
	SegWgStruct hst(rc);
	(*this) = QuepField(hst, wx, wz);
}

/* full QUEP 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,
   wz: horizontal computational window,
   Mx: number of expansion terms per horizontal layer,
   suppress log output, if quiet == 1 */
QuepField::QuepField(Circuit rc, Polarization p, 
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{
	if(quiet != 1)
	{
fprintf(stderr, "\n------------- QUEP ------------------------------------------------- '04 ---\n");
fprintf(stderr, "T%c  ", polCHR(p));
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d,  ", wx.x0, wx.x1, Mx);
fprintf(stderr, "z in (%.10g, %.10g), Mz: %d\n", wz.x0, wz.x1, Mz);
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) = QuepField(rc, wx, wz);
	discretize(p, Mx, Mz, 0);
}
QuepField::QuepField(Circuit rc, Polarization p, 
                     Interval wx, int Mx, Interval wz, int Mz)
{
	(*this) = QuepField(rc, p, wx, Mx, wz, Mz, 0);
}
QuepField::QuepField(SegWgStruct rc, Polarization p, 
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{
	(*this) = QuepField(rc.circuit(), p, wx, Mx, wz, Mz, quiet);
}
QuepField::QuepField(SegWgStruct rc, Polarization p, 
                     Interval wx, int Mx, Interval wz, int Mz)
{
	(*this) = QuepField(rc.circuit(), p, wx, Mx, wz, Mz, 0);
}

/* destroy */
QuepField::~QuepField()
{
	h.clear();
	v.clear();
}

/* copy constructor */
QuepField::QuepField(const QuepField& q)
{
	h = q.h;
	v = q.v;
	cwx = q.cwx;
	cwz = q.cwz;
	QF_pol    = q.QF_pol;
	QF_lambda = q.QF_lambda;
	QF_vM     = q.QF_vM;
	QF_ky     = q.QF_ky;
}
	
/* assignment */
QuepField& QuepField::operator=(const QuepField& q)
{
	if(this != &q)
	{
		h = q.h; 
		v = q.v; 
		cwx = q.cwx;
		cwz = q.cwz;
		QF_pol    = q.QF_pol;
		QF_lambda = q.QF_lambda;
		QF_vM     = q.QF_vM;
		QF_ky     = q.QF_ky;
	}
	return *this;
}

/* delete all segment entries */
void QuepField::clear()
{
	h.clear();
	v.clear();
}

/* input from FILE dat */
void QuepField::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		quepfielderror("QuepField: read");
	}
	else
	{
		h.read(dat);
		v.read(dat);
		cwx.read(dat);
		cwz.read(dat);
		check();
	}
	return;
}

/* output to FILE dat */
void QuepField::write(FILE *dat)
{
	comment("QuepField", dat);
	comment("BepField, horizontal", dat);
	h.write(dat);
	comment("BepField, vertical", dat);
	v.write(dat);
	comment("CW, vertical", dat);
	cwx.write(dat);
	comment("CW, horizontal", dat);
	cwz.write(dat);
	return;
}

/* input from default file */
void QuepField::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] = 'q';
	name[7] = 'e';
	name[8] = 'p';
	name[9] = 0;

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

/* output to default file */
void QuepField::write(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] = 'q';
	name[7] = 'e';
	name[8] = 'p';
	name[9] = 0;

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

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

/* the waveguides assiciated with the borders of the circuit */
Waveguide QuepField::port(CBorder bdr) const
{
	switch(bdr)
	{
	case RIGHT:  return h.s(h.s.nz+1);
	case TOP:    return v.s(v.s.nz+1);
	case LEFT:   return h.s(0);
	case BOTTOM: return v.s(0);
	}
	quepfielderror("port, boundary type");
	return h.s(h.s.nz+1);
}

/* the arrays of in- and outgoing modes crossing the circuit borders */
ModeArray QuepField::modes(CBorder bdr) const
{
	switch(bdr)
	{
	case RIGHT:  return h(h.s.nz+1);
	case TOP:    return v(v.s.nz+1);
	case LEFT:   return h(0);
	case BOTTOM: return v(0);
	}
	quepfielderror("modes, boundary type");
	return h(h.s.nz+1);
}

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

/* mode m with amplitude A at border bdr */
void QuepField::input(CBorder bdr, Complex A, const Mode& m)
{
	switch(bdr)
	{
	case RIGHT:  h.input(SBRIGHT, A, m); return; 
	case TOP:    v.input(SBRIGHT, A, m); return; 
	case LEFT:   h.input(SBLEFT, A, m); return; 
	case BOTTOM: v.input(SBLEFT, A, m); return; 
	}
	return;
}
/* mode of order o with amplitude A at border bdr */
void QuepField::input(CBorder bdr, Complex A, int o)
{
	switch(bdr)
	{
	case RIGHT:  h.input(SBRIGHT, A, o); return; 
	case TOP:    v.input(SBRIGHT, A, o); return; 
	case LEFT:   h.input(SBLEFT, A, o); return; 
	case BOTTOM: v.input(SBLEFT, A, o); return; 
	}
	return;
}
/* tilted Gaussian beam gb, multiplied by amplitude A, at border bdr */
void QuepField::input(CBorder bdr, const Gaussianbeam& gb, Complex A)
{
	switch(bdr)
	{
	case RIGHT:  h.input(SBRIGHT, gb, A); return; 
	case TOP:    v.input(SBRIGHT, gb, A); return; 
	case LEFT:   h.input(SBLEFT, gb, A); return; 
	case BOTTOM: v.input(SBLEFT, gb, A); return; 
	}
	return;
}
/* ... a normalized beam, unspecified phase */
void QuepField::input(CBorder bdr, const Gaussianbeam& gb)
{
	input(bdr, gb, CC1);
	return;
}

/* --- output amplitudes & powers ------------------------------------ */
/*     & variants for vTE / vTM modes, specification of Polarization   */

/* amplitude of outgoing mode m at border bdr */
Complex QuepField::Aout(CBorder bdr, const Mode& m) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Aout(SBRIGHT, m);
	case TOP:    return v.Aout(SBRIGHT, m); 
	case LEFT:   return h.Aout(SBLEFT, m);
	case BOTTOM: return v.Aout(SBLEFT, m);
	}
	return CC0;
}

/* amplitude of outgoing mode of order o at border bdr */
Complex QuepField::Aout(CBorder bdr, int o) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Aout(SBRIGHT, o);
	case TOP:    return v.Aout(SBRIGHT, o); 
	case LEFT:   return h.Aout(SBLEFT, o);
	case BOTTOM: return v.Aout(SBLEFT, o);
	}
	return CC0;
}

/* vM: amplitude of outgoing mode of order o at border bdr, 
   with polarization pol */
Complex QuepField::Aout(CBorder bdr, Polarization pol, int o) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Aout(SBRIGHT, pol, o);
	case TOP:    return v.Aout(SBRIGHT, pol, o); 
	case LEFT:   return h.Aout(SBLEFT, pol, o);
	case BOTTOM: return v.Aout(SBLEFT, pol, o);
	}
	return CC0;
}

/* vM: amplitude of outgoing mode of order o passing border bdr, 
   with polarization pol, evaluated at propagation coordinate p */
Complex QuepField::Aout(CBorder bdr, Polarization pol, int o, double p) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Aout(SBRIGHT, pol, o, p);
	case TOP:    return v.Aout(SBRIGHT, pol, o, p); 
	case LEFT:   return h.Aout(SBLEFT, pol, o, p);
	case BOTTOM: return v.Aout(SBLEFT, pol, o, p);
	}
	return CC0;
}

/* vM: amplitude of incoming mode of order o passing border bdr, 
   with polarization pol, evaluated at propagation coordinate p */
Complex QuepField::Ain(CBorder bdr, Polarization pol, int o, double p) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Ain(SBRIGHT, pol, o, p);
	case TOP:    return v.Ain(SBRIGHT, pol, o, p); 
	case LEFT:   return h.Ain(SBLEFT, pol, o, p);
	case BOTTOM: return v.Ain(SBLEFT, pol, o, p);
	}
	return CC0;
}


/* output power associated with mode m at border bdr */
double QuepField::Pout(CBorder bdr, const Mode& m) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pout(SBRIGHT, m);
	case TOP:    return v.Pout(SBRIGHT, m); 
	case LEFT:   return h.Pout(SBLEFT, m);
	case BOTTOM: return v.Pout(SBLEFT, m);
	}
	return 0.0;
}

/* output power associated with mode of order o at border bdr */
double QuepField::Pout(CBorder bdr, int o) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pout(SBRIGHT, o);
	case TOP:    return v.Pout(SBRIGHT, o); 
	case LEFT:   return h.Pout(SBLEFT, o);
	case BOTTOM: return v.Pout(SBLEFT, o);
	}
	return 0.0;
}

/* output power associated with mode of order o at border bdr 
   with polarization pol (vM) */
double QuepField::Pout(CBorder bdr, Polarization pol, int o) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pout(SBRIGHT, pol, o);
	case TOP:    return v.Pout(SBRIGHT, pol, o); 
	case LEFT:   return h.Pout(SBLEFT, pol, o);
	case BOTTOM: return v.Pout(SBLEFT, pol, o);
	}
	return 0.0;
}

/* guided output power crossing border bdr */
double QuepField::Pgout(CBorder bdr) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pgout(SBRIGHT);
	case TOP:    return v.Pgout(SBRIGHT); 
	case LEFT:   return h.Pgout(SBLEFT);
	case BOTTOM: return v.Pgout(SBLEFT);
	}
	return 0.0;
}

/* guided output power crossing border bdr, 
   associated with polarization pol (vM) */
double QuepField::Pgout(CBorder bdr, Polarization pol) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pgout(SBRIGHT, pol);
	case TOP:    return v.Pgout(SBRIGHT, pol); 
	case LEFT:   return h.Pgout(SBLEFT, pol);
	case BOTTOM: return v.Pgout(SBLEFT, pol);
	}
	return 0.0;
}


/* total output power crossing border bdr (guided and nonguided fields) */
double QuepField::Pout(CBorder bdr) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pout(SBRIGHT);
	case TOP:    return v.Pout(SBRIGHT); 
	case LEFT:   return h.Pout(SBLEFT);
	case BOTTOM: return v.Pout(SBLEFT);
	}
	return 0.0;
}

/* total output power crossing border bdr (guided and nonguided fields)
   associated with polarization pol (vM) */
double QuepField::Pout(CBorder bdr, Polarization pol) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pout(SBRIGHT, pol);
	case TOP:    return v.Pout(SBRIGHT, pol); 
	case LEFT:   return h.Pout(SBLEFT, pol);
	case BOTTOM: return v.Pout(SBLEFT, pol);
	}
	return 0.0;
}

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

/* total input power crossing border bdr (guided and nonguided fields) */
double QuepField::Pin(CBorder bdr) const 
{
	switch(bdr)
	{
	case RIGHT:  return h.Pin(SBRIGHT);
	case TOP:    return v.Pin(SBRIGHT); 
	case LEFT:   return h.Pin(SBLEFT);
	case BOTTOM: return v.Pin(SBLEFT);
	}
	return 0.0;
}

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

/* inspect modal expansions & amplitudes at the border regions */
void QuepField::inspect(CBorder bdr) const
{
	ModeArray ma = modes(bdr);
	Cvector ain, aout;
	char bid = ' ';
	switch(bdr)
	{
	case LEFT: 
		ain  = h.f(0, SBRIGHT);
		aout = h.b(0, SBRIGHT);
		bid = 'L';
		break;
	case RIGHT: 
		aout = h.f(h.s.nz+1, SBLEFT);
		ain  = h.b(h.s.nz+1, SBLEFT);
		bid = 'R';
		break;
	case BOTTOM: 
		ain  = v.f(0, SBRIGHT);
		aout = v.b(0, SBRIGHT);
		bid = 'B';
		break;
	case TOP: 
		aout = v.f(v.s.nz+1, SBLEFT);
		ain  = v.b(v.s.nz+1, SBLEFT);
		bid = 'T';
		break;
	}
	for(int j=0; j<=ma.num-1; ++j)
        fprintf(stderr, "%c[%d] %s bq = %g  |A_in| = %g  |A_out| = %g\n",
	        bid, j, ma(j).ids, ma(j).betaq, ain(j).abs(), aout(j).abs());
	return;
}

void QuepField::inspect() const
{
	inspect(LEFT); fprintf(stderr, "\n");
	inspect(TOP); fprintf(stderr, "\n");
	inspect(RIGHT); fprintf(stderr, "\n");
	inspect(BOTTOM); fprintf(stderr, "\n");
	return;
}


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

/* check the QuepField for consistency 
   & set values for polarization, wavelength, vectorfields  */
void QuepField::check()
{
	QF_vM = (h(0))(0).vM;
	QF_ky = (h(0))(0).ky;
	QF_pol = (h(0))(0).pol;
	QF_lambda = (h(0))(0).wg.lambda;
	ModeArray mai;
	for(int i=0; i<=h.s.nz+1; ++i)
	{
		mai = h(i);
		int nmi = mai.num;
		for(int j=0; j<=nmi-1; ++j)
		{
			Mode mij = mai(j);
			if(mij.wg.lambda != QF_lambda) 
				quepfielderror("check: wavelengths mixed");
//			if(mij.bdt_b != LIMD && mij.special!=2) 
//				quepfielderror("check: boundary type (_b)");
//			if(mij.bdt_t != LIMD && mij.special!=2) 
//				quepfielderror("check: boundary type (_t)");
			if(fabs(mij.bp_b-cwx.x0)> COMPTOL_HX && mij.special!=2)
			      quepfielderror("check: mode boundary mismatch");
			if(fabs(mij.bp_t-cwx.x1)> COMPTOL_HX && mij.special!=2) 
			      quepfielderror("check: mode boundary mismatch");
		}
	}
	for(int i=0; i<=v.s.nz+1; ++i)
	{
		mai = v(i);
		int nmi = mai.num;
		for(int j=0; j<=nmi-1; ++j)
		{
			Mode mij = mai(j);
			if(mij.wg.lambda != QF_lambda) 
				quepfielderror("check: wavelengths mixed");
//			if(mij.bdt_b != LIMD && mij.special!=2) 
//				quepfielderror("check: boundary type (_b)");
//			if(mij.bdt_t != LIMD && mij.special!=2) 
//				quepfielderror("check: boundary type (_t)");
			if(fabs(mij.bp_b-cwz.x0)>COMPTOL_HX && mij.special!=2) 
				quepfielderror("check: mode boundary mismatch");
			if(fabs(mij.bp_t-cwz.x1)>COMPTOL_HX && mij.special!=2) 
				quepfielderror("check: mode boundary mismatch");
		}
	}
	if(fabs(h.s.hz(0)-cwz.x0)>COMPTOL_HX) 
		quepfielderror("check: segment-cw-boundary mismatch");
	if(fabs(h.s.hz(h.s.nz)-cwz.x1)>COMPTOL_HX) 
		quepfielderror("check: segment-cw-boundary mismatch");
	if(fabs(v.s.hz(0)-cwx.x0)>COMPTOL_HX) 
		quepfielderror("check: segment-cw-boundary mismatch");
	if(fabs(v.s.hz(v.s.nz)-cwx.x1)>COMPTOL_HX) 
		quepfielderror("check: segment-cw-boundary mismatch");

	if(QF_vM)
	{
		for(int i=0; i<=h.s.nz+1; ++i)
		{
			mai = h(i);
			int nmi = mai.num;
			for(int j=0; j<=nmi-1; ++j)
			{
				Mode mij = mai(j);
				if(mij.vM != QF_vM) 
					quepfielderror("check: vM not unique");
				if((mij.ky-QF_ky).sqabs() > SLAMS_BQTOL) 
					quepfielderror("check: ky not unique");
			}
		}
		for(int i=0; i<=v.s.nz+1; ++i)
		{
			mai = v(i);
			int nmi = mai.num;
			for(int j=0; j<=nmi-1; ++j)
			{
				Mode mij = mai(j);
				if(mij.vM != QF_vM) 
					quepfielderror("check: vM not unique");
				if((mij.ky-QF_ky).sqabs() > SLAMS_BQTOL) 
					quepfielderror("check: ky not unique");
			}
		}
		return;
	}
	for(int i=0; i<=h.s.nz+1; ++i)
	{
		mai = h(i);
		int nmi = mai.num;
		for(int j=0; j<=nmi-1; ++j)
		{
			Mode mij = mai(j);
			if(mij.pol != QF_pol) 
				quepfielderror("check: polarizations mixed");
		}
	}
	for(int i=0; i<=v.s.nz+1; ++i)
	{
		mai = v(i);
		int nmi = mai.num;
		for(int j=0; j<=nmi-1; ++j)
		{
			Mode mij = mai(j);
			if(mij.pol != QF_pol) 
				quepfielderror("check: polarizations mixed");
		}
	}
	return;
}

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

/* compute the discretized quadridirectional mode spectra:
   modes of polarization p,
   Mx expansion terms per vertical slice,
   Mz terms per horizontal layer */
void QuepField::discretize(Polarization p, int Mx, int Mz, int quiet)
{
	if(quiet != 1) fprintf(stderr, "Horiz");
	h.discretize(p, cwx, Mx, quiet);
	if(quiet != 1) fprintf(stderr, "Vertc");
	v.discretize(p, cwz, Mz, quiet);
	check();
	return;
}

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

/* quadridirectional mode interference evaluation */

/* field superposition at point x, z 
   cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex QuepField::field(Fcomp cp, double x, double z) const
{
	Complex ex, ey, ez, hx, hy, hz, f;
	if(QF_vM)
	{
		switch(cp)
		{
		case EX: return h.field(EX, x, z)-v.field(EZ, z, x); 
		case EY: return h.field(EY, x, z)-v.field(EY, z, x); 
		case EZ: return h.field(EZ, x, z)-v.field(EX, z, x); 
		case HX: return h.field(HX, x, z)+v.field(HZ, z, x); 
		case HY: return h.field(HY, x, z)+v.field(HY, z, x); 
		case HZ: return h.field(HZ, x, z)+v.field(HX, z, x); 
		case SX: 
			 ey = field(EY, x, z);
			 hz = field(HZ, x, z);
			 ez = field(EZ, x, z);
			 hy = field(HY, x, z);
			 f = ey*hz.conj()-ez*hy.conj();
			 return Complex(0.5*f.re, 0.0); 
			 break;
		case SY: 
			 ez = field(EZ, x, z);
			 hx = field(HX, x, z);
			 ex = field(EX, x, z);
			 hz = field(HZ, x, z);
			 f = ez*hx.conj()-ex*hz.conj();
			 return Complex(0.5*f.re, 0.0); 
			 break;
		case SZ: 
			 ey = field(EY, x, z);
			 hx = field(HX, x, z);
			 ex = field(EX, x, z);
			 hy = field(HY, x, z);
			 f = ex*hy.conj()-ey*hx.conj();
			 return Complex(0.5*f.re, 0.0); 
			 break;
		default: return CC0;
			 break;
		}
	}
	switch(cp)
	{
	case EX: return h.field(EX, x, z)-v.field(EZ, z, x); 
	case EY: return h.field(EY, x, z)+v.field(EY, z, x); 
	case EZ: return h.field(EZ, x, z)-v.field(EX, z, x); 
	case HX: return h.field(HX, x, z)-v.field(HZ, z, x); 
	case HY: return h.field(HY, x, z)+v.field(HY, z, x); 
	case HZ: return h.field(HZ, x, z)-v.field(HX, z, x); 
	case SX: 
	         ey = field(EY, x, z);
	         hz = field(HZ, x, z);
	         ez = field(EZ, x, z);
	         hy = field(HY, x, z);
	         f = ey*hz.conj()-ez*hy.conj();
	         return Complex(0.5*f.re, 0.0); 
	         break;
	case SY: 
	         ez = field(EZ, x, z);
	         hx = field(HX, x, z);
	         ex = field(EX, x, z);
	         hz = field(HZ, x, z);
	         f = ez*hx.conj()-ex*hz.conj();
	         return Complex(0.5*f.re, 0.0); 
	         break;
	case SZ: 
	         ey = field(EY, x, z);
	         hx = field(HX, x, z);
	         ex = field(EX, x, z);
	         hy = field(HY, x, z);
	         f = ex*hy.conj()-ey*hx.conj();
	         return Complex(0.5*f.re, 0.0); 
	         break;
	default: return CC0;
	         break;
	}
	return CC0;
}

/* local field "strength",  id: one of mE, mH, qE, qH, mW */
double QuepField::lfs(FSid id, double x, double z) const
{
	Complex ex, ey, ez, hx, hy, hz;
	double eps_xz;
	switch(id)
	{
	case mE:
		ex = field(EX, x, z);
		ey = field(EY, x, z);
		ez = field(EZ, x, z);
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
	case mH:
		hx = field(HX, x, z);
		hy = field(HY, x, z);
		hz = field(HZ, x, z);
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
	case qE:
		ex = field(EX, x, z);
		ey = field(EY, x, z);
		ez = field(EZ, x, z);
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
	case qH:
		hx = field(HX, x, z);
		hy = field(HY, x, z);
		hz = field(HZ, x, z);
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
	case mW:
		ex = field(EX, x, z);
		ey = field(EY, x, z);
		ez = field(EZ, x, z);
		hx = field(HX, x, z);
		hy = field(HY, x, z);
		hz = field(HZ, x, z);
		eps_xz = h.s(h.s.segidx(z)).eps(x);
		return  0.25*(  MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			+EP0*eps_xz*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
	default:
		quepfielderror("lfs, FSid?");
                break;
	}
	return 0.0;
}

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

/* the six electric and magnetic components at x, z */
void QuepField::emfield(double x, double z,
                        Complex& ex, Complex& ey, Complex& ez,
                        Complex& hx, Complex& hy, Complex& hz) const
{
	if(QF_vM) quepfielderror("emfield: vM not yet implemented");
	Complex hex, hey, hez, hhx, hhy, hhz;
	Complex vex, vey, vez, vhx, vhy, vhz;
	h.emfield(x, z, hex, hey, hez, hhx, hhy, hhz);
	v.emfield(z, x, vex, vey, vez, vhx, vhy, vhz);
	ex = hex-vez; 
	ey = hey+vey; 
	ez = hez-vex; 
	hx = hhx-vhz; 
	hy = hhy+vhy; 
	hz = hhz-vhx; 
	return;
}


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

/* quadridirectional mode expansion solver:
   given the structure with the accompanying mode expansions,
   where input mode amplitudes are specified by
   h.setf(0, f), h.setb(h.s.nz+1, b),
   v.setf(0, u), v.setb(v.s.nz+1, d), or the input methods,
   find the remaining mode coefficients. */
void QuepField::quepsim() 
{
	if(h.s.nz < 3) quepfielderror("quepsim: #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[h.s.nz];
        hIfb = new Cmatrix[h.s.nz];
        hIbf = new Cmatrix[h.s.nz];
        hIbb = new Cmatrix[h.s.nz];        
        hTX  = new Cvector[h.s.nz];        
	Cvector hT1, hTN;
	Cmatrix hD0, hD1, hDN, hDNp;
	if(v.s.nz < 3) quepfielderror("quepsim: #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[v.s.nz];
        vIfb = new Cmatrix[v.s.nz];
        vIbf = new Cmatrix[v.s.nz];
        vIbb = new Cmatrix[v.s.nz];        
        vTX  = new Cvector[v.s.nz];        
	Cvector vT1, vTN;
	Cmatrix vD0, vD1, vDN, vDNp;
	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);     
	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", h.s.nz+2);
	int rfmin, rfmax, rfn;
	rfmin = rfmax = h(0).num;
	for(int i=1; i<=h.s.nz+1; ++i)
	{
		rfn = 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", v.s.nz+2);
	rfmin = rfmax = v(0).num;
	for(int i=1; i<=v.s.nz+1; ++i)
	{
		rfn = 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 = v(1).num;
	dvb = v(v.s.nz).num;

	Cmatrix hR0mff, hR0mfb, hR0mbf, hR0mbb; 
	dh  = h(0).num;
	vmat = v.crossoverlap(h(0), 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(h(1).equal(h(0)) != 1)
	{
		dh  = h(1).num;
		vmat = v.crossoverlap(h(1), 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  = h(h.s.nz).num;
	vmat = v.crossoverlap(h(h.s.nz), h.s.hz(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(h(h.s.nz+1).equal(h(h.s.nz)) != 1)
	{
		dh  = h(h.s.nz+1).num;
		vmat = v.crossoverlap(h(h.s.nz+1), h.s.hz(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 = h(1).num;
	dhb = h(h.s.nz).num;

	Cmatrix vR0mff, vR0mfb, vR0mbf, vR0mbb; 
	dv  = v(0).num;
	vmat = h.crossoverlap(v(0), 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(v(1).equal(v(0)) != 1)
	{
		dv  = v(1).num;
		vmat = h.crossoverlap(v(1), 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  = v(v.s.nz).num;
	vmat = h.crossoverlap(v(v.s.nz), v.s.hz(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(v(v.s.nz+1).equal(v(v.s.nz)) != 1)
	{
		dv  = v(v.s.nz+1).num;
		vmat = h.crossoverlap(v(v.s.nz+1), v.s.hz(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 = h(1).num;
	dbn = h(h.s.nz).num;
	du1 = v(1).num;
	ddn = v(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, ".:");

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

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

	fprintf(stderr, ".");
	q.strip();
	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(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));
	h.setb(0, SBRIGHT, sol);

	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));
	h.setf(h.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));
	v.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));
	v.setf(v.s.nz+1, SBLEFT, sol);

	fprintf(stderr, ".");

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

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

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

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

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

/* QUEP simulations as before, with a test whether the total relative 
   power balance is violated by not more than pviol,
   return value 1: ok, 0: no success. */
int QuepField::quepsim(double pviol)
{
	quepsim();
	double pin  =  Pin(LEFT)+ Pin(RIGHT)+ Pin(BOTTOM)+ Pin(TOP);
	double pout = Pout(LEFT)+Pout(RIGHT)+Pout(BOTTOM)+Pout(TOP);
	fprintf(stderr, "P_in=%g, P_out=%g, ", pin, pout);
	if(pout>pin*(1.0-pviol) && pout<pin*(1.0+0.5*pviol))
	{
		fprintf(stderr, "Acc.\n");
		return 1;
	}
	fprintf(stderr, "Rej.\n");
	return 0;
}

/* QUEP simulations as before, but with a movable computational window:
   try at most numtr times to obtain a result with a violation of the 
   total relative power balance of not more than pviol,
   by repeating the calculation with computational window boundary
   positions shifted by small integer multiples of +/- bshift times the 
   window width, return value 1: ok, 0: no success. */
int QuepField::quepsim_s(int numtr, double pviol, double bshift)
{
	int done = 0;
	int tr = 0;
	bshift = fabs(bshift);
	Interval cwx0 = cwx;
	Interval cwz0 = cwz;
	double dbx0 = bshift*cwx0.len();
	double dbz0 = bshift*cwz0.len();
	double dbx = dbx0;
	double dbz = dbz0;
	while(done == 0 && tr <= numtr-1)
	{
		quepsim();
		double pin  =  Pin(LEFT)+ Pin(RIGHT)+ Pin(BOTTOM)+ Pin(TOP);
		double pout = Pout(LEFT)+Pout(RIGHT)+Pout(BOTTOM)+Pout(TOP);
		fprintf(stderr, "P_in=%g, P_out=%g, ", pin, pout);
		if(pout>pin*(1.0-pviol) && pout<pin*(1.0+0.5*pviol))
		{
			fprintf(stderr, "Acc(%d).\n", tr);
			done = 1;
		}
		else
		{
		    fprintf(stderr, "Rej(%d).\n", tr);
		    if(tr <= numtr-2)
		    {
			Interval ncwx(cwx0.x0-dbx, cwx0.x1+dbx);
			Interval ncwz(cwz0.x0-dbz, cwz0.x1+dbz);

			Interval ix(v.s.hz(1), v.s.hz(v.s.nz-1));
			Interval iz(h.s.hz(1), h.s.hz(h.s.nz-1));
			double l0, l;
			l0 =  cwx0.len();
			l  = ncwx.len();
			if(ncwx.x0 >= ix.x0) 
				ncwx.x0 = ix.x0-fabs(ix.x0-cwx0.x0)*l/l0;
			if(ncwx.x1 <= ix.x1) 
				ncwx.x1 = ix.x1+fabs(cwx0.x1-ix.x1)*l/l0;
			l0 =  cwz0.len();
			l  = ncwz.len();
			if(ncwz.x0 >= iz.x0) 
				ncwz.x0 = iz.x0-fabs(iz.x0-cwz0.x0)*l/l0;
			if(ncwz.x1 <= iz.x1) 
				ncwz.x1 = iz.x1+fabs(cwz0.x1-iz.x1)*l/l0;
			fprintf(stderr, 
				"x in (%.10g, %.10g), z in (%.10g, %.10g)\n",
				ncwx.x0, ncwx.x1, ncwz.x0, ncwz.x1);
			SegWgStruct st = h.s.shrink();
			QuepField q(st, ncwx, ncwz);;
			if(QF_vM)
			{
				int Mx = h(0).num/2;
				int Mz = v(0).num/2;
				q.discretize(QF_ky, Mx, Mz, 0);
			}
			else
			{
				int Mx = h(0).num;
				int Mz = v(0).num;
				Polarization p = (h(0))(0).pol;
				q.discretize(p, Mx, Mz, 0);
			}
			Cvector a;

			a = h.f(0, SBRIGHT); 
			for(int i=0; i<=h(0).num-1; ++i) 
				a(i) *= ((h(0))(i)).efac(FORW, ncwz.x0-cwz.x0);
			q.h.setf(0, SBRIGHT, a);

			a = h.b(h.s.nz+1, SBLEFT); 
			for(int i=0; i<=h(h.s.nz+1).num-1; ++i) 
				a(i) *= ((h(h.s.nz+1))(i)).efac(BACK, ncwz.x1-cwz.x1);
			q.h.setb(h.s.nz+1, SBLEFT, a);

			a = v.f(0, SBRIGHT); 
			for(int i=0; i<=v(0).num-1; ++i) 
				a(i) *= ((v(0))(i)).efac(FORW, ncwx.x0-cwx.x0);
			q.v.setf(0, SBRIGHT, a);

			a = v.b(v.s.nz+1, SBLEFT); 
			for(int i=0; i<=v(v.s.nz+1).num-1; ++i) 
				a(i) *= ((v(v.s.nz+1))(i)).efac(BACK, ncwx.x1-cwx.x1);
			q.v.setb(v.s.nz+1, SBLEFT, a);

			(*this) = q;
			if(dbx > 0.0) dbx = -dbx; 
			else          dbx = -dbx+dbx0;
			if(dbz > 0.0) dbz = -dbz; 
			else          dbz = -dbz+dbz0;
		    }
		}
		++tr;
	}
	return done;
}

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

/* prepare for the vTE/vTM mode representation, for given ky */
void QuepField::vectorify(Complex ky)
{
	h.vectorify(ky);	
	v.vectorify(ky);
	return;
}

/* vQUEP: compute the discretized quadridirectional mode spectra:
   vTE and vTM modes, uniform wavenumber ky,
   2*Mx expansion terms per vertical slice,
   2*Mz terms per horizontal layer */
void QuepField::discretize(Complex ky, int Mx, int Mz, int quiet)
{
	if(quiet != 1) fprintf(stderr, "Horiz");
	h.discretize(ky, 0, cwx, Mx, quiet);
	if(quiet != 1) fprintf(stderr, "Vertc");
	v.discretize(ky, 1, cwz, Mz, quiet);
	check();
	return;
}

/* full vQUEP initialization: setup the container for the quadridirectional 
   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,
   wz: horizontal computational window,
   Mx: number of expansion terms per horizontal layer,
   suppress log output, if quiet == 1 */
QuepField::QuepField(Circuit rc, Complex ky, 
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{
	if(quiet != 1)
	{
fprintf(stderr, "------------- vQUEP ------------------------------------------------ '14 ---\n");
fprintf(stderr, "vTE, vTM, ky = (%g + i %g) k0\n", 
                ky.re/val_k0(rc.lambda), ky.im/val_k0(rc.lambda));
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d,  ", wx.x0, wx.x1, Mx);
fprintf(stderr, "z in (%.10g, %.10g), Mz: %d\n", wz.x0, wz.x1, Mz);
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) = QuepField(rc, wx, wz);
	discretize(ky, Mx, Mz, 0);
}
QuepField::QuepField(Circuit rc, Complex ky, 
                     Interval wx, int Mx, Interval wz, int Mz)
{
	(*this) = QuepField(rc, ky, wx, Mx, wz, Mz, 0);
}
QuepField::QuepField(SegWgStruct rc, Complex ky, 
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{
	(*this) = QuepField(rc.circuit(), ky, wx, Mx, wz, Mz, quiet);
}
QuepField::QuepField(SegWgStruct rc, Complex ky, 
                     Interval wx, int Mx, Interval wz, int Mz)
{
	(*this) = QuepField(rc.circuit(), ky, wx, Mx, wz, Mz, 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) */
QuepField::QuepField(Circuit rc, CBorder bdr,
                     Polarization pol, int ord, Complex amp, double theta,
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{
	(*this) = QuepField(rc, wx, wz);
	Waveguide ewg = port(bdr);
	ModeArray ma;
	modeanalysis(ewg, pol, ma, 1);
	if(ord < 0 || ord >= ma.num) quepfielderror("vQUEP, setup, mode order");
	double N = ma(ord).neff;
	double b = ma(ord).beta;
	if(theta <= -90.0 || theta >= 90.0) quepfielderror("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; amp = %g+i%g; angle = %g^o\n", 
                ma(ord).ids, N, b, amp.re, amp.im, theta);
	}
	(*this) = QuepField(rc, ky, wx, Mx, wz, Mz, quiet);
	switch(bdr)
	{
	case LEFT: case RIGHT:
		if(pol == TE) input(bdr, amp, ord);
		else          input(bdr, amp, Mx+ord);
		break;
	case TOP: case BOTTOM:
		if(pol == TE) input(bdr, amp, ord);
		else          input(bdr, amp, Mz+ord);
		break;
	}
}
QuepField::QuepField(Circuit rc, CBorder bdr,
                     Polarization pol, int ord, Complex amp, double theta,
                     Interval wx, int Mx, Interval wz, int Mz)
{	
	(*this) = QuepField(rc, bdr, pol, ord, amp, theta,
                            wx, Mx, wz, Mz, 0);
}
QuepField::QuepField(SegWgStruct rc, CBorder bdr,
                     Polarization pol, int ord, Complex amp, double theta,
                     Interval wx, int Mx, Interval wz, int Mz, int quiet)
{	
	(*this) = QuepField(rc.circuit(), bdr, pol, ord, amp, theta,
                            wx, Mx, wz, Mz, quiet);
}
QuepField::QuepField(SegWgStruct rc, CBorder bdr,
                     Polarization pol, int ord, Complex amp, double theta,
                     Interval wx, int Mx, Interval wz, int Mz)
{	
	(*this) = QuepField(rc.circuit(), bdr, pol, ord, amp, theta,
                            wx, Mx, wz, Mz, 0);
}


	


/* ------------------------------------------------------------------- */
/* 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),
   the internal representaion is modified accordingly */
void QuepField::adjustphase(double phi)
{
	h.adjustphase(phi);
	v.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 QuepField::adjustphase(double x, double z)
{
	if(QF_vM) quepfielderror("adjustphase, vM: specify a component");
	Fcomp cp = principalcomp(QF_pol);
	adjustphase(cp, 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 the computational window */
void QuepField::adjustphase(int numx, int numz)
{
	if(QF_vM) quepfielderror("adjustphase, vM: specify a component");
	Fcomp cp = principalcomp(QF_pol);
	adjustphase(cp, numx, numz);
	return;
}

/* ... as before, with (coarse) defaults for numx, numz */
void QuepField::adjustphase()
{
	if(QF_vM) quepfielderror("adjustphase, vM: specify a component");
	Fcomp cp = principalcomp(QF_pol);
	adjustphase(cp);
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   component cp exhibits the maximum field at position (x, z) */
void QuepField::adjustphase(Fcomp cp, double x, double z)
{
	adjustphase(-field(cp, x, z).arg());
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   component cp 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 the computational window */
void QuepField::adjustphase(Fcomp cp, int numx, int numz)
{
	double x, z;
	double xs, zs;
	xs = cwx.len()/((double) numx);
	zs = cwz.len()/((double) numz);
	x = cwx.x0;
	double xm = 0.0, zm = 0.0, maxf = 0.0;
	double f;
	for(int xi=0; xi <=numx; ++xi)
	{
		z = cwz.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 QuepField::adjustphase(Fcomp cp)
{
	int numx, numz;	
	numx = (int) ceil(sqrt(2500.0*cwx.len()/cwz.len()));
	numz = numx*((int) ceil(cwz.len()/cwx.len()));
	adjustphase(cp, numx, numz);
	return;
}

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

/* the optical field, discretized on a regular rectangular mesh,  
   npx x npz points on a display window dwx x dwz
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   foa: ORG, MOD, SQR, REP, IMP (ORG = REP) */
Cmatrix QuepField::field(Fcomp cp, Interval dwx, int npx, 
                                   Interval dwz, int npz) const
{
	if(npx <= 1 || npz <= 1) quepfielderror("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 QuepField::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 QuepField::field(Fcomp cp, double x0, double z0, 
                         double x1, double z1, int np) const
{
	if(np <= 1) quepfielderror("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 QuepField::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);
}

/* ------------------------------------------------------------------- */
/* 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 QuepField::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);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, h.s, 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_lowlevcolormap(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 QuepField::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);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Phasemap");
	mlout_gengeoxz(dat, h.s, 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;
}

/* 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 QuepField::plot(FSid id,
                     double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	if(npx <= 1 || npz <= 1) 
		quepfielderror("plot, lfs: npx, npz"); 
	FILE *dat;
	char name[13] = "pl____I.m";
	name[2] = mqchr(id);
	name[3] = idchr(id);
	name[4] = ext0;
	name[5] = ext1;
	fprintf(stderr, "plot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, h.s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Dmatrix fld(npx, npz);
	double x, dx, z, dz;
	dx = (xend-xbeg)/(npx-1);
	dz = (zend-zbeg)/(npz-1);
	z = zbeg;
	for(int j=0; j<=npz-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) 
		{
			fld(i, j) = lfs(id, x, z);
			x += dx;
		}
		z += dz;
	}
	mlout_fld(dat, npx, npz, SZ, fld);
	name[7] = 0; 
	char desc[10];
	iddesc(id, desc);
	mlout_genimage(SZ, MOD, name, desc, dat); 
	mlout_lowlevcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* electromagnetic energy density image
   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 QuepField::plot(double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	return plot(mW, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
}

/* 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 QuepField::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);
	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 st = h.s.circuit();
	int nsm = 2*(st.nz+2)*(st.nx+1)+2*(st.nx+2)*(st.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*(st.nz+2)*npx+2*(st.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<=st.nx+1; ++l)
	{
		if(l==0)
		{
			if(xbeg < st.hx(0)) x0 = xbeg; 	
			else                x0 = st.hx(0)-COMPTOL_HX;
		}
                else x0 = st.hx(l-1);
		if(l==st.nx+1)
		{
			if(xend > st.hx(st.nx)) x1 = xend; 	
			else                    x1 = st.hx(st.nx)+COMPTOL_HX; 	
		}
                else x1 = st.hx(l);
		Interval r(x0, x1);
		int o = dwx.clip(r);
		if(o)
		{
		    for(int s=0; s<=st.nz; ++s)
		    {
			if(fabs(st.n(l,s)-st.n(l,s+1)) > COMPTOL_N)
			{
			    zp = st.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<=st.nz+1; ++s)
	{
		if(s==0)
		{
			if(zbeg < st.hz(0)) z0 = zbeg; 	
			else                z0 = st.hz(0)-COMPTOL_HX;
		}
                else z0 = st.hz(s-1);
		if(s==st.nz+1)
		{
			if(zend > st.hz(st.nz)) z1 = zend; 	
			else                    z1 = st.hz(st.nz)+COMPTOL_HX; 	
		}
                else z1 = st.hz(s);
		Interval r(z0, z1);
		int o = dwz.clip(r);
		if(o)
		{
		    for(int l=0; l<=st.nx; ++l)
		    {
			if(fabs(st.n(l,s)-st.n(l+1,s)) > COMPTOL_N)
			{
			    xp = st.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 images of component cp, or EY (TE) | HY (TM), resp.,
   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 QuepField::movie(Fcomp cp, 
                      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;
        name[2] = fldchr(cp);
        name[3] = cpchr(cp);      
	name[4] = ext0;
	name[5] = ext1;
	tmax = val_period_T(QF_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);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, h.s, 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(QF_lambda), famp); 
	fclose(dat);
	return;
}
void QuepField::movie(double xbeg, double xend, double zbeg, double zend,
                      int npx, int npz, int ntfr, char ext0, char ext1) const
{
	if(QF_vM) quepfielderror("movie, vM: specify a component");
	Fcomp cp = principalcomp(QF_pol);
	movie(cp, xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1);
	return;
}
		
/* animation of the light propagation, fancy style ;-)
   the frames show images of component cp, or EY (TE) | HY (TM), resp.,
   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 QuepField::fmovie(Fcomp cp, 
                       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;
        name[2] = fldchr(cp);
        name[3] = cpchr(cp);      
	name[4] = ext0;
	name[5] = ext1;
	tmax = val_period_T(QF_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);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, h.s, 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 st = h.s.circuit();
	int nsm = 2*(st.nz+2)*(st.nx+1)+2*(st.nx+2)*(st.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*(st.nz+2)*npx+2*(st.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<=st.nx+1; ++l)
	{
		if(l==0)
		{
			if(xbeg < st.hx(0)) x0 = xbeg; 	
			else                x0 = st.hx(0)-COMPTOL_HX;
		}
                else x0 = st.hx(l-1);
		if(l==st.nx+1)
		{
			if(xend > st.hx(st.nx)) x1 = xend; 	
			else                    x1 = st.hx(st.nx)+COMPTOL_HX; 	
		}
                else x1 = st.hx(l);
		Interval r(x0, x1);
		int o = dwx.clip(r);
		if(o)
		{
		    for(int s=0; s<=st.nz; ++s)
		    {
			if(fabs(st.n(l,s)-st.n(l,s+1)) > COMPTOL_N)
			{
			    zp = st.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<=st.nz+1; ++s)
	{
		if(s==0)
		{
			if(zbeg < st.hz(0)) z0 = zbeg; 	
			else                z0 = st.hz(0)-COMPTOL_HX;
		}
                else z0 = st.hz(s-1);
		if(s==st.nz+1)
		{
			if(zend > st.hz(st.nz)) z1 = zend; 	
			else                    z1 = st.hz(st.nz)+COMPTOL_HX; 	
		}
                else z1 = st.hz(s);
		Interval r(z0, z1);
		int o = dwz.clip(r);
		if(o)
		{
		    for(int l=0; l<=st.nx; ++l)
		    {
			if(fabs(st.n(l,s)-st.n(l+1,s)) > COMPTOL_N)
			{
			    xp = st.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(QF_lambda), famp);
	fclose(dat);
	return;
}
void QuepField::fmovie(double xbeg, double xend, double zbeg, double zend,
                       int npx, int npz, int ntfr, char ext0, char ext1) const
{
	if(QF_vM) quepfielderror("fmovie, vM: specify a component");
	Fcomp cp = principalcomp(QF_pol);
	fmovie(cp, 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 QuepField::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";
	name[3] = ext0;
	name[4] = ext1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   xbeg, npx, xend, zbeg, npz, zend, name);
	dat = fopen(name, "w+");
	name[6] = 0; 
	mlout_viewertop(dat, name, QF_pol, QF_lambda, QF_vM);
	mlout_gengeoxz(dat, h.s, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix f;
	Fcomp cp;
	if(QF_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(QF_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) = h.s.n(x, z); x += dx; }
		z += dz;
	}
	mlout_n(dat, npx, npz, n);
	mlout_fldviewer(dat, name);
	fclose(dat);
	return;
}