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

/*
 * cmt.cpp
 * Procedures related to the variants of coupled mode theory
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"cylfunc.h"
#include"slamode.h"
#include"slamarr.h"
#include"matlvis.h"
#include"cylfunc.h"
#include"wgm.h"
#include"bend.h"
#include"fim.h"
#include"cmt.h"

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

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

/*
 * hybrid CMT:
 * analytical modal basis fields,
 * combined with FEM-discretized amplitude functions
 */

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

/* component identifiers for use by the HCMT solver, translation */
Fcomp fkcp2fcomp(FKcp cp)
{
	switch(cp)
	{
	case fcEx: return EX; break;		
	case fcEy: return EY; break;		
	case fcEz: return EZ; break;		
	case fcHx: return HX; break;		
	case fcHy: return HY; break;		
	case fcHz: return HZ; break;		
	default: ;
	}
	cmterror("fkcp2fcomp: cp");
	return EY;
}

FKcp fcomp2fkcp(Fcomp cp)
{
	switch(cp)
	{
	case EX: return fcEx; break;		
	case EY: return fcEy; break;		
	case EZ: return fcEz; break;		
	case HX: return fcHx; break;		
	case HY: return fcHy; break;		
	case HZ: return fcHz; break;		
	default: ;
	}
	cmterror("fcomp2fkcp: cp");
	return fcEy;
}

/* ... for use by the vectorial (3D) HCMT solver, translation */
Fcomp c3cp2fcomp(C3cp cp)
{
	switch(cp)
	{
	case f3Ex: return EX; break;		
	case f3Ey: return EY; break;		
	case f3Ez: return EZ; break;		
	case f3Hx: return HX; break;		
	case f3Hy: return HY; break;		
	case f3Hz: return HZ; break;		
	default: ;
	}
	cmterror("c3cp2fcomp: cp");
	return EY;
}

C3cp fcomp2c3cp(Fcomp cp)
{
	switch(cp)
	{
	case EX: return f3Ex; break;		
	case EY: return f3Ey; break;		
	case EZ: return f3Ez; break;		
	case HX: return f3Hx; break;		
	case HY: return f3Hy; break;		
	case HZ: return f3Hz; break;		
	default: ;
	}
	cmterror("fcomp2c3cp: cp");
	return f3Ey;
}


/* HCMT basic field types -> characters */
void hcmtbftchr(HcmtBft t, char& c0, char& c1, char &c2)
{
	switch(t)
	{
	case MHF: c0='M'; c1='H'; c2='F'; break;
	case MHB: c0='M'; c1='H'; c2='B'; break;
	case MVU: c0='M'; c1='V'; c2='U'; break;
	case MVD: c0='M'; c1='V'; c2='D'; break;
	case SRC: c0='S'; c1='R'; c2='C'; break;
	case GBM: c0='G'; c1='B'; c2='M'; break;
	case RES: c0='R'; c1='E'; c2='S'; break;
	case CCC: c0='C'; c1='C'; c2='C'; break;
	case CCA: c0='C'; c1='C'; c2='A'; break;
	case WGM: c0='W'; c1='G'; c2='M'; break;
	case MAT: c0='M'; c1='A'; c2='T'; break;
	case FIM: c0='F'; c1='I'; c2='M'; break;
	default: break;
	}
	return;
}

/* for evaluation of Hcmt integrals & mode overlaps:
   parts of guided mode profiles will be neglected where the field strength 
   is below MCROP times the mode field profile maximum */
double MCROP = 0.001;

/* for evaluation of Hcmt integrals & mode overlaps:
   contributions of bend mode profiles will be neglected 
   for radii below  BM_R_CROP_IN and radii beyond BM_R_CROP_OUT */
double BM_R_CROP_IN  = 0.001;
double BM_R_CROP_OUT = 100.0;

/* for evaluation of Hcmt integrals & mode overlaps:
   contributions of WGM profiles will be neglected 
   for radii below  WM_R_CROP_IN and radii beyond WM_R_CROP_OUT */
double WM_R_CROP_IN  = 0.001;
double WM_R_CROP_OUT = 100.0;


/* ------------------------------------------------------------------------ 
   segments of the x-z-plane, for piecewise integration <-> FEM matrices 
   ------------------------------------------------------------------------ */

/* part of the initialization */
void Isegment::set_za_zb_al(Interval cwx, Interval cwz)
{
	cwx0 = cwx.x0;
	cwx1 = cwx.x1;
	cwz0 = cwz.x0;
	cwz1 = cwz.x1;
	cal = cos(alpha);
	sal = sin(alpha);
	if(fabs(cal) < 1.0e-8) 
	{
		za = cwz0;
		zb = cwz1;
		return;
	}
	double za0, za1, zb0, zb1;
	za0 = z0+sa*cal-sal/cal*(cwx0-x0-sa*sal);
	za1 = z0+sa*cal-sal/cal*(cwx1-x0-sa*sal);
	zb0 = z0+sb*cal-sal/cal*(cwx0-x0-sb*sal);
	zb1 = z0+sb*cal-sal/cal*(cwx1-x0-sb*sal);
	za = za0;	
	if(za1<za) za = za1;
	if(zb0<za) za = zb0;
	if(zb1<za) za = zb1;
	zb = za0;	
	if(za1>zb) zb = za1;
	if(zb0>zb) zb = zb0;
	if(zb1>zb) zb = zb1;
	if(za < cwz0) za = cwz0;
	if(zb > cwz1) zb = cwz1;
	return;
}

int Isegment::set_LPC_slopes(double &s1, double &s2)
{
	double r = sqrt( (cwx1-cwx0)*(cwx1-cwx0)+(cwz1-cwz0)*(cwz1-cwz0));
	double x1, z1, x2, z2, z;
	x1 = x0 + r*cos(s1);
	z1 = z0 + r*sin(s1);
	x2 = x0 + r*cos(s2);
	z2 = z0 + r*sin(s2);
	if(fabs(z1-z0)<1.0e-6)
	{
		if(x1>x0)
		{
			ubt = CST; ubcx = cwx1;
			lbsg = (x2-x0)/(z2-z0);
			lbt = LPC;
			z = z0; if(z2 < z) z = z2;
			if(z<cwz0) z = cwz0;
			za = z;
			z = z0; if(z2 > z) z = z2;
			if(z>cwz1) z = cwz1;
			zb = z;
			return 0;
		}
		else
		{
			lbt = CST; lbcx = cwx0;
			ubsg = (x2-x0)/(z2-z0);
			ubt = LPC;
			z = z0; if(z2 < z) z = z2;
			if(z<cwz0) z = cwz0;
			za = z;
			z = z0; if(z2 > z) z = z2;
			if(z>cwz1) z = cwz1;
			zb = z;
			return 0;
		}
		
	}
	if(fabs(z2-z0)<1.0e-6)
	{
		if(x2>x0)
		{
			ubt = CST; ubcx = cwx1;
			lbsg = (x1-x0)/(z1-z0);
			lbt = LPC;
			z = z0; if(z1 < z) z = z1;
			if(z<cwz0) z = cwz0;
			za = z;
			z = z0; if(z1 > z) z = z1;
			if(z>cwz1) z = cwz1;
			zb = z;
			return 0;
		}
		else
		{
			lbt = CST; lbcx = cwx0;
			ubsg = (x1-x0)/(z1-z0);
			ubt = LPC;
			z = z0; if(z1 < z) z = z1;
			if(z<cwz0) z = cwz0;
			za = z;
			z = z0; if(z1 > z) z = z1;
			if(z>cwz1) z = cwz1;
			zb = z;
			return 0;
		}
	}
	if((z1<z0 && z2>z0) || (z2<z0 && z1>z0))
	{	
		if(x1>x0 || x2>x0)
		{
			double a, b;
			a = s1; b = 0.0;
			set_LPC_slopes(a, b);
			s1 = 0.0; 
			return 1;
		}
		if(x1<x0 || x2<x0)
		{
			double a, b;
			a = s1; b = PI;
			set_LPC_slopes(a, b);
			s1 = PI; 
			return 1;
		}
		cmterror("Isegment::set_LPC_slopes: confused");
		return 0;
	}
	if(x1>x2)
	{
		ubsg = (x1-x0)/(z1-z0);
		ubt = LPC;
		lbsg = (x2-x0)/(z2-z0);
		lbt = LPC;
		z = z0; if(z1 < z) z = z1; if(z2 < z) z = z2;
                if(z<cwz0) z = cwz0;
                za = z;
                z = z0; if(z1 > z) z = z1; if(z2 > z) z = z2;
                if(z>cwz1) z = cwz1;
                zb = z;
		return 0;
	}
	else
	{
		lbsg = (x1-x0)/(z1-z0);
		lbt = LPC;
		ubsg = (x2-x0)/(z2-z0);
		ubt = LPC;
		z = z0; if(z1 < z) z = z1; if(z2 < z) z = z2;
                if(z<cwz0) z = cwz0;
                za = z;
                z = z0; if(z1 > z) z = z1; if(z2 > z) z = z2;
                if(z>cwz1) z = cwz1;
                zb = z;
		return 0;
	}
	return 0;
}


/* the lower boundary */
double Isegment::lbxpos(double z) const 
{
	double dz, dx;
	double xa, xb, x;
	switch(lbt)
	{
	case CST:
		return lbcx;
	case CRV:
		dz = z-lbcz;
		dx = sqrt(fabs(lbr*lbr-dz*dz));
		return lbcx+lbsg*dx;
	case STL:
		if(fabs(sal) < 1.0e-8) return cwx0;
		xa = x0+sa*sal+cal/sal*(z0+sa*cal-z);
		xb = x0+sb*sal+cal/sal*(z0+sb*cal-z);
		x = xa; if(xb<x) x = xb;
		if(x < cwx0) x = cwx0;
		if(x > cwx1) x = cwx1;
		return x;
	case LPC:
		x = x0+lbsg*(z-z0);
		if(x < cwx0) x = cwx0;
		if(x > cwx1) x = cwx1;
		return x;
	}
	return 0.0;
}

/* the upper boundary */
double Isegment::ubxpos(double z) const
{
	double dz, dx;
	double xa, xb, x;
	switch(ubt)
	{
	case CST:
		return ubcx;
	case CRV:
		dz = z-ubcz;
		dx = sqrt(fabs(ubr*ubr-dz*dz));
		return ubcx+ubsg*dx;
	case STL:
		if(fabs(sal) < 1.0e-8) return cwx1;
		xa = x0+sa*sal+cal/sal*(z0+sa*cal-z);
		xb = x0+sb*sal+cal/sal*(z0+sb*cal-z);
		x = xa; if(xb>x) x = xb;
		if(x > cwx1) x = cwx1;
		if(x < cwx0) x = cwx0;
		return x;
	case LPC:
		x = x0+ubsg*(z-z0);
		if(x < cwx0) x = cwx0;
		if(x > cwx1) x = cwx1;
		return x;
	}
	return 0.0;
}
	
/* for inspection purposes: save segment information to a file */
void Isegment::write(FILE *dat)
{
	fprintf(dat, "[za, zb] = [%g, %g]\n", za, zb);
	switch(lbt)
	{
	case CST:
		fprintf(dat, "lb: CST, x = %g\n", lbcx);
		break;
	case CRV:
		fprintf(dat, "lb: CRV, r = %g, sg = %g\n", 
			      lbr, lbsg);
		fprintf(dat, "    xs = %g, zs = %g\n", 
			      lbcx, lbcz);
		break;
	case STL:
		fprintf(dat, "lb: STL, sa = %g, sb = %g\n", 
			      sa, sb);
		fprintf(dat, "    x0 = %g, z0 = %g, al = %g\n", 
			      x0, z0, alpha);
	case LPC:
		fprintf(dat, "lb: LPC, slope = %g\n", lbsg);
		break;
	default: break;
	}
	switch(ubt)
	{
	case CST:
		fprintf(dat, "ub: CST, x = %g\n", ubcx);
		break;
	case CRV:
		fprintf(dat, "ub: CRV, r = %g, sg = %g\n", 
			      ubr, ubsg);
		fprintf(dat, "    xs = %g, zs = %g\n", 
			      ubcx, ubcz);
		break;
	case STL:
		fprintf(dat, "ub: STL\n");
		break;
	case LPC:
		fprintf(dat, "ub: LPC, slope = %g\n", ubsg);
		break;
	default: break;
	}
	return;
}

/* for inspection purposes: MATLAB plot commands for the segment outline */
void Isegment::plot(FILE *dat, char *s)
{
	Dvector zp(203), xp(203);
	double z;
	int zi=0;
	for(int i=0; i<=100; ++i)
	{
		z = za+((double)i)*(zb-za)/100.0;
		zp(zi) = z;
		xp(zi) = lbxpos(z);
		++zi;
	}
	for(int i=100; i>=0; --i)
	{
		z = za+((double)i)*(zb-za)/100.0;
		zp(zi) = z;
		xp(zi) = ubxpos(z);
		++zi;
	}
	zp(zi) = z;
	xp(zi) = lbxpos(z);
	++zi;
	if(zi != 203) cmterror("Isegment:plot");
	fprintf(dat, "plot([ ");
	for(int i=0; i<=202; ++i) fprintf(dat, "%g ", zp(i));
	fprintf(dat, " ], [ ");
	for(int i=0; i<=202; ++i) fprintf(dat, "%g ", xp(i));
	fprintf(dat, " ], '%s');\n", s);
	return;
}

/* ------------------------------------------------------------------------ 
   basis elements for the hybrid CMT, one of the interacting fields 
   ------------------------------------------------------------------------ */

/* initialization: generic */
HcmtElement::HcmtElement()
{
	vM = 0;
	ky = 0.0;
	type = MHF;
	pol = TE;	
	k0 = val_k0(1.0);
	invommu0 = val_invommu0(1.0);
	invomep0 = val_invomep0(1.0);
	dir = FORW;
	s0 = -1.0; 
	sN = 1.0;
	ds_hom = 0.01;
	ds_disc = 0.01;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
	pramp.strip();
	prpos.strip();
}

/* destroy */
HcmtElement::~HcmtElement()
{
	if(seg != NULL)
	{

		for(int k=0; k<=Ns; ++k) 
		{
			if(seg[k] != NULL) 
			{
				delete[] seg[k];
				seg[k] = NULL;
			}
		}	
		delete[] seg;
	}
	seg = NULL;
	nseg.strip();
}

/* copy constructor */
HcmtElement::HcmtElement(const HcmtElement& he)
{
	vM = he.vM;
	ky = he.ky;
	pol = he.pol;
	k0 = he.k0;
	invommu0 = he.invommu0;
	invomep0 = he.invomep0;
	type = he.type;
	ain = he.ain;
	aout = he.aout;
	s0 = he.s0;
	sN = he.sN;
	ds_hom = he.ds_hom;
	ds_disc = he.ds_disc;
	u = he.u;
	sbeg = he.sbeg;
	scen = he.scen;
	send = he.send;
	sdsL = he.sdsL;
	sdsR = he.sdsR;
	given = he.given;
	recEli = he.recEli;
	mode = he.mode;
	rmode = he.rmode;
	bmode = he.bmode;
	wgmode = he.wgmode;
	fimode = he.fimode;
	mext = he.mext;
	dir = he.dir;
	d_beta = he.d_beta;
	d_alpha = he.d_alpha;
	d_kappa = he.d_kappa;
	d_gamma = he.d_gamma;
	d_kz = he.d_kz;
	v_z0 = he.v_z0;
	ord = he.ord;
	ordd = he.ordd;
	orda = he.orda;
	x0 = he.x0;
	xn = he.xn;
	numx = he.numx;
	dx = he.dx;
	z0 = he.z0;
	zn = he.zn;
	numz = he.numz;
	dz = he.dz;
	baB = he.baB;
	width = he.width;
	alpha = he.alpha;
	sal = he.sal;
	cal = he.cal;
	nb = he.nb;
	kn = he.kn;
	kqnq = he.kqnq;
	FDfEy = he.FDfEy;
	FDfHx = he.FDfHx;
	FDfHz = he.FDfHz;
	FDfdzEy = he.FDfdzEy;
	FDfdxEy = he.FDfdxEy;
	FDfdzHxmdxHz = he.FDfdzHxmdxHz;
	Ns = he.Ns;
	nseg = he.nseg;
	seg = NULL;
	if(he.seg != NULL)
	{
		seg = new Isegment*[Ns+1];
		for(int k=0; k<=Ns; ++k)
		{
			seg[k] = new Isegment[MaxNumIsegs];
			for(int j=0; j<=nseg(k)-1; ++j)
				seg[k][j] = he.seg[k][j];
		}
	}
	pramp = he.pramp;
	prpos = he.prpos;
}

/* assignment */
HcmtElement& HcmtElement::operator=(const HcmtElement& he)
{
        if(this != &he)
        {
		if(seg != NULL)
		{
			for(int k=0; k<=Ns; ++k) 
			{
				if(seg[k] != NULL) 
				{
					delete[] seg[k];
					seg[k] = NULL;
				}
			}	
			delete[] seg;
		}

		vM = he.vM;
		ky = he.ky;
		pol = he.pol;
		k0 = he.k0;
		invommu0 = he.invommu0;
		invomep0 = he.invomep0;
		type = he.type;
		ain = he.ain;
		aout = he.aout;
		s0 = he.s0;
		sN = he.sN;
		ds_hom = he.ds_hom;
		ds_disc = he.ds_disc;
		u = he.u;
		sbeg = he.sbeg;
		scen = he.scen;
		send = he.send;
		sdsL = he.sdsL;
		sdsR = he.sdsR;
		given = he.given;
		recEli = he.recEli;
		mode = he.mode;
		rmode = he.rmode;
		bmode = he.bmode;
		wgmode = he.wgmode;
		fimode = he.fimode;
		mext = he.mext;
		dir = he.dir;
		d_beta = he.d_beta;
		d_alpha = he.d_alpha;
		d_kappa = he.d_kappa;
		d_gamma = he.d_gamma;
		d_kz = he.d_kz;
		v_z0 = he.v_z0;
		ord = he.ord;
		ordd = he.ordd;
		orda = he.orda;
		x0 = he.x0;
		xn = he.xn;
		numx = he.numx;
		dx = he.dx;
		z0 = he.z0;
		zn = he.zn;
		numz = he.numz;
		dz = he.dz;
		baB = he.baB;
		width = he.width;
		alpha = he.alpha;
		sal = he.sal;
		cal = he.cal;
		nb = he.nb;
		kn = he.kn;
		kqnq = he.kqnq;
		FDfEy = he.FDfEy;
		FDfHx = he.FDfHx;
		FDfHz = he.FDfHz;
		FDfdzEy = he.FDfdzEy;
		FDfdxEy = he.FDfdxEy;
		FDfdzHxmdxHz = he.FDfdzHxmdxHz;
		Ns = he.Ns;
		nseg = he.nseg;
		seg = NULL;
		if(he.seg != NULL)
		{
			seg = new Isegment*[Ns+1];
			for(int k=0; k<=Ns; ++k)
			{
				seg[k] = new Isegment[MaxNumIsegs];
				for(int j=0; j<=nseg(k)-1; ++j) 
				{
					seg[k][j] = he.seg[k][j];
				}
			}
		}
		pramp = he.pramp;
		prpos = he.prpos;
        }
        return *this;
}

/* write (most) element information to file */
void HcmtElement::write(FILE *dat)
{
	fprintf(dat, "   // Element //\n");
	fprintf(dat, "vM: %d\n", vM);
	fprintf(dat, "ky: %.16g\n", ky);
	fprintf(dat, "pol: %d\n", pol);
	fprintf(dat, "k0: %.16g\n", k0);
	fprintf(dat, "invommu0: %.16g\n", invommu0);
	fprintf(dat, "invomep0: %.16g\n", invomep0);
	fprintf(dat, "type: %d\n", type);
	fprintf(dat, "ain: %.16g+i%.16g\n", ain.re, ain.im);
	fprintf(dat, "aout: %.16g+i%.16g\n", aout.re, aout.im);
	fprintf(dat, "s0: %.16g\n", s0);
	fprintf(dat, "sN: %.16g\n", sN);
	fprintf(dat, "ds_hom: %.16g\n", ds_hom);
	fprintf(dat, "ds_disc: %.16g\n", ds_disc);
	fprintf(dat, "u:\n"); u.write(dat);
	fprintf(dat, "sbeg:\n"); sbeg.write(dat);
	fprintf(dat, "scen:\n"); scen.write(dat);
	fprintf(dat, "send:\n"); send.write(dat);
	fprintf(dat, "sdsL:\n"); sdsL.write(dat);
	fprintf(dat, "sdsR:\n"); sdsR.write(dat);
	fprintf(dat, "given:\n"); given.write(dat);
	fprintf(dat, "recEli: %d\n", recEli);
	fprintf(dat, "mode: ???\n"); // mode.write(dat);
	fprintf(dat, "rmode: ???\n"); // rmode.write(dat);
	fprintf(dat, "bmode: ???\n"); // bmode.write(dat);
	fprintf(dat, "wgmode: ???\n"); // wgmode.write(dat);
	fprintf(dat, "fimode: ???\n"); // fimode.write(dat);
	fprintf(dat, "mext:\n"); mext.write(dat);
	fprintf(dat, "dir: %d\n", dir);
	fprintf(dat, "d_beta: %.16g\n", d_beta);
	fprintf(dat, "d_alpha: %.16g\n", d_alpha);
	fprintf(dat, "d_kappa: %.16g\n", d_kappa);
	fprintf(dat, "d_gamma: %.16g+i%.16g\n", d_gamma.re, d_gamma.im);
	fprintf(dat, "d_kz: %.16g+i%.16g\n", d_kz.re, d_kz.im);
	fprintf(dat, "v_z0: %.16g\n", v_z0);
	fprintf(dat, "ord: %d\n", ord);
	fprintf(dat, "ordd: %.16g\n", ordd);
	fprintf(dat, "orda: %d\n", orda);
	fprintf(dat, "x0: %.16g\n", x0);
	fprintf(dat, "xn: %.16g\n", xn);
	fprintf(dat, "numx: %d\n", numx);
	fprintf(dat, "dx: %.16g\n", dx);
	fprintf(dat, "z0: %.16g\n", z0);
	fprintf(dat, "zn: %.16g\n", zn);
	fprintf(dat, "numz: %d\n", numz);
	fprintf(dat, "dz: %.16g\n", dz);
	fprintf(dat, "baB: %.16g\n", baB);
	fprintf(dat, "width: %.16g\n", width);
	fprintf(dat, "alpha: %.16g\n", alpha);
	fprintf(dat, "sal: %.16g\n", sal);
	fprintf(dat, "cal: %.16g\n", cal);
	fprintf(dat, "nb: %.16g\n", nb);
	fprintf(dat, "kn: %.16g\n", kn);
	fprintf(dat, "kqnq: %.16g\n", kqnq);
	fprintf(dat, "FDfEy:\n"); FDfEy.write(dat);
	fprintf(dat, "FDfHx:\n"); FDfHx.write(dat);
	fprintf(dat, "FDfHz:\n"); FDfHz.write(dat);
	fprintf(dat, "FDfdzEy:\n"); FDfdzEy.write(dat);
	fprintf(dat, "FDfdxEy:\n"); FDfdxEy.write(dat);
	fprintf(dat, "FDfdzHxmdxHz:\n"); FDfdzHxmdxHz.write(dat);
	fprintf(dat, "Ns: %d\n", Ns);
	fprintf(dat, "nseg:\n"); nseg.write(dat);
	for(int k=0; k<=Ns; ++k)
	{
		for(int j=0; j<=nseg(k)-1; ++j) (seg[k][j]).write(dat);
	}
	fprintf(dat, "pramp:\n"); pramp.write(dat);
	fprintf(dat, "prpos:\n"); prpos.write(dat);
	return;
}

/* vHCMT, evanescent slab modes:
   amplitude discretization will be discontinued, if the absolute 
   slab mode level drops below EvSM_DISC_CROP */
double EvSM_DISC_CROP = 1.0e-8;

/* initialization, guided modes
   t: MHF or MHB for modes that travel forward / backward along z
      MVU or MVD for modes that travel upward / downward along x
   m: the mode object
   FEM discretization of the amplitude as a function of
   the x- or z-coordinate in [s_0, s_n], stepsize d_s, regular */
HcmtElement::HcmtElement(HcmtBft t, Mode m, double s_0, double s_n, double d_s)
{
	switch(t)
	{
	case MHF: case MHB: case MVU: case MVD: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, mode");
		break;
	}	
	if(m.vM == 1) 
	{
		vM = 1;
		if(m.ky.im != 0.0) cmterror("HcmtElement constructor, mode, vM, Im ky");
		ky = m.ky.re;
	}
	else 
	{
		vM = 0;
		ky = 0.0;
	}
	pol = m.pol;	
	k0 = m.k0;
	invommu0 = m.invommu0;
	invomep0 = m.invomep0;
	mode = m;
	dir = FORW;
	if(type == MHB || type == MVD) dir = BACK;
	d_beta = (dir == FORW) ? mode.beta : -mode.beta;
	if(m.vM == 1) 
	{
		if(dir == FORW)
		{
			d_kz = mode.kz;
			v_z0 = 0.5*(s_0+s_n);
			d_beta  =  mode.kz.re;
			d_alpha = -mode.kz.im;
/*
			double sb = s_0;
			double se = s_n;
			double r;
			r = expi(-d_kz*(se-v_z0)).abs()/expi(-d_kz*(sb-v_z0)).abs();
			while(r < EvSM_DISC_CROP && fabs(se-sb) > 3.0*d_s)
			{
				sb += d_s;
				se -= d_s;
				r = expi(-d_kz*(se-v_z0)).abs()/expi(-d_kz*(sb-v_z0)).abs();
			}
			s_0 = sb;
			s_n = se;
*/
			
		}
		else 
		{
			d_kz = -mode.kz;
			v_z0 = 0.5*(s_0+s_n);
			d_beta  = -mode.kz.re;
			d_alpha =  mode.kz.im;
/*
			double sb = s_0;
			double se = s_n;
			double r;
			r = expi(-d_kz*(sb-v_z0)).abs()/expi(-d_kz*(se-v_z0)).abs();
			while(r < EvSM_DISC_CROP && fabs(se-sb) > 3.0*d_s)
			{
				sb += d_s;
				se -= d_s;
				r = expi(-d_kz*(sb-v_z0)).abs()/expi(-d_kz*(se-v_z0)).abs();
			}
			s_0 = sb;
			s_n = se;
*/
		}
	}
	mext = mode.extent(MCROP);
	s0 = s_0;
	sN = s_n;
	ds_hom = d_s;
	ds_disc = d_s;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
}

/* initialization, guided modes
   t: MHF or MHB for modes that travel forward / backward along z
      MVU or MVD for modes that travel upward / downward along x
   m: the mode object
   FEM discretization of the amplitude as a function of
   the x- or z-coordinate in [s_0, s_n], irregular,
   elements of maximum length d_s_disc symmetrically positioned at each
   respective interface in str, in-between a regular discretization
   with maximum stepsize d_s_hom */
HcmtElement::HcmtElement(HcmtBft t, Mode m, Circuit str, 
                         double s_0, double s_n, 
                         double d_s_hom, double d_s_disc)
{
	switch(t)
	{
	case MHF: case MHB: case MVU: case MVD: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, mode");
		break;
	}	
	if(m.vM == 1) 
	{
		vM = 1;
		if(m.ky.im != 0.0) cmterror("HcmtElement constructor, mode, vec, Im ky");
		ky = m.ky.re;
	}
	else 
	{
		vM = 0;
		ky = 0.0;
	}
	pol = m.pol;	
	k0 = m.k0;
	invommu0 = m.invommu0;
	invomep0 = m.invomep0;
	mode = m;
	dir = FORW;
	if(type == MHB || type == MVD) dir = BACK;
	d_beta = (dir == FORW) ? mode.beta : -mode.beta;
	mext = mode.extent(MCROP);
	s0 = s_0;
	sN = s_n;
	ds_hom = d_s_hom;
	ds_disc = d_s_disc;
	Dvector p;
	switch(t)
	{
	case MHF: case MHB: 
		p = str.hz;
		break;
	case MVU: case MVD: 
		p= str.hx;
		break;
	default: 
		cmterror("HcmtElement constructor, mode");
		break;
	}	
	seg = NULL;
	initFEMstr(p);
}

HcmtElement::HcmtElement(HcmtBft t, Mode m, SegWgStruct str, 
                         double s_0, double s_n, 
                         double d_s_hom, double d_s_disc)
{
	Circuit cstr = str.circuit(); 
	switch(t)
	{
	case MHF: case MHB: case MVU: case MVD: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, mode");
		break;
	}	
	if(m.vM == 1) 
	{
		vM = 1;
		if(m.ky.im != 0.0) cmterror("HcmtElement constructor, mode, vM, Im ky");
		ky = m.ky.re;
	}
	else 
	{
		vM = 0;
		ky = 0.0;
	}
	pol = m.pol;	
	k0 = m.k0;
	invommu0 = m.invommu0;
	invomep0 = m.invomep0;
	mode = m;
	dir = FORW;
	if(type == MHB || type == MVD) dir = BACK;
	d_beta = (dir == FORW) ? mode.beta : -mode.beta;
	mext = mode.extent(MCROP);
	s0 = s_0;
	sN = s_n;
	ds_hom = d_s_hom;
	ds_disc = d_s_disc;
	Dvector p;
	switch(t)
	{
	case MHF: case MHB: 
		p = cstr.hz;
		break;
	case MVU: case MVD: 
		p= cstr.hx;
		break;
	default: 
		cmterror("HcmtElement constructor, mode");
		break;
	}	
	seg = NULL;
	initFEMstr(p);
}

/* initialization, localized source of polarization p, wavelength wl,
   outgoing Hankel functions H_2 of order o, localized at x_0, z_0,
   embedded in a medium with refractive index n,
   t: SRC,
   FEM discretization of the amplitude as a function of
   the radial coordinate in [s_0, s_n], stepsize d_s */
HcmtElement::HcmtElement(HcmtBft t, Polarization p, double wl,
                         int o, double x_0, double z_0, double n,
                         double s_0, double s_n, double d_s)
{
	if(p != TE) cmterror("TM polarization: not yet implemented");
	switch(t)
	{
	case SRC: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, source");
		break;
	}	
	vM = 0;
	ky = 0.0;
	pol = p;	
	k0 = val_k0(wl);
	invommu0 = val_invommu0(wl);
	invomep0 = val_invomep0(wl);
	ord = o;
	ordd = ((double)ord);
	orda = abs(ord);
	x0 = x_0;
	z0 = z_0;
	nb = n;
	kn = k0*nb; 
	kqnq = kn*kn; 
	s0 = s_0;
	sN = s_n;
	ds_hom = d_s;
	ds_disc = d_s;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
}

/* initialization, Gaussian beam with polarization p, wavelength wl,
   width w at the waist, beam waist localized at x_0, z_0, a beam that
   propagates at an angle a with respect to the pos. z-axis in a
   medium with refractive index n,
   t: GBM,
   FEM discretization of the amplitude as a function of
   the coordinate along the beam axis in [s_0, s_n], stepsize d_s */
HcmtElement::HcmtElement(HcmtBft t, Polarization p, double wl,
                         double w, double x_0, double z_0, double a, double n,
                         double s_0, double s_n, double d_s)
{
	switch(t)
	{
	case GBM: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, Gaussian beam");
		break;
	}	
	vM = 0;
	ky = 0.0;
	pol = p;	
	k0 = val_k0(wl);
	invommu0 = val_invommu0(wl);
	invomep0 = val_invomep0(wl);
	width = w,
	x0 = x_0;
	z0 = z_0;
	alpha = a;
	nb = n;
	kn = k0*nb; 
	kqnq = kn*kn; 
	s0 = s_0;
	sN = s_n;
	ds_hom = d_s;
	ds_disc = d_s;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
	if(pol == TE)
	{
		double p0 = invommu0/8.0*sqrt(2.0*PI)
			    *(2.0*width*width*kqnq-1.0)/(kn*width);
		baB = 1.0/sqrt(p0);
	}
	else
	{
		double p0 = invomep0/8.0/(nb*nb)*sqrt(2.0*PI)
			    *(2.0*width*width*kqnq-1.0)/(kn*width);
		baB = 1.0/sqrt(p0);
	}
	sal = sin(alpha);
	cal = cos(alpha);
}

/* initialization, bimodal resonance of a dieletric square,
   established by the cores of modes m1, m2, center location at x_0, z_0,
   a single resonance amplitude is introduced */
HcmtElement::HcmtElement(HcmtBft t, Mode m1, Mode m2, double x_0, double z_0)
{
	switch(t)
	{
	case RES: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, resonance");
	}
	if(m1.pol != TE) cmterror("TM polarization: not yet implemented");
	if(m1.pol != m2.pol) cmterror("RES element: polarization mismatch");
	if(m1.wg.equal(m2.wg) != 1) cmterror("RES element: waveguide mismatch");
	if(m1.wg.nx != 1) cmterror("RES element: only one core layer allowed");
	if(m1.wg.checksymmetry() != 1) cmterror("RES element: symmetry");
	if(m1.vM != 0 || m2.vM != 0) cmterror("RES element: vM not implemented");
	 mode = m1;
	rmode = m2;
	vM = 0;
	ky = 0.0;
	pol = mode.pol;	
	k0 = mode.k0;
	invommu0 = mode.invommu0;
	invomep0 = mode.invomep0;
	 mode.translate(-( mode.wg.hx(0)+ mode.wg.hx(1))/2.0);
	rmode.translate(-(rmode.wg.hx(0)+rmode.wg.hx(1))/2.0);
	x0 = x_0;
	z0 = z_0;
	width = (mode.wg.hx(1)-mode.wg.hx(0))/2.0;
	Interval e, re;
	 e =  mode.extent(MCROP);
	re = rmode.extent(MCROP);
	mext.x0 = e.x0; if(re.x0 < mext.x0) mext.x0=re.x0;
	mext.x1 = e.x1; if(re.x1 > mext.x1) mext.x1=re.x1;
	if(fabs(mext.x0+mext.x1)>1.0e-10) cmterror("RES element: geometry");
	seg = NULL;
	Ns = 0;
	u = Cvector(Ns+1);
	u.init(CC0);
}

/* initialization, bend mode traveling round a circular cavity
   t: CCC or CCA for modes that travel clockwise / anticlockwise
   m: the bend mode object
   FEM discretization of the amplitude as a function of
   the angular coordinate in [p_0, p_n], stepsize d_p, regular
   kR: prescribed approximate integer angular wavenumber (optional) */
HcmtElement::HcmtElement(HcmtBft t, BDMode m, 
                         double x_0, double z_0, double d_p, int kR)
{
	switch(t)
	{
	case CCC: case CCA:
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, circular cavity");
	}

	bmode = m;
	vM = 0;
	ky = 0.0;
	pol = bmode.pol;	
	k0 = bmode.k0;
	invommu0 = val_invommu0(bmode.wg.lambda);
	invomep0 = val_invomep0(bmode.wg.lambda);
	x0 = x_0;
	z0 = z_0;
	dir = FORW;
	if(type == CCA) dir = BACK;
	d_gamma = (dir == FORW) ? bmode.c_gamma() : bmode.c_gamma()*(-1.0);
	d_beta =  (dir == FORW) ? bmode.beta  : -bmode.beta;
	d_alpha = (dir == FORW) ? bmode.alpha : -bmode.alpha;
	double k;
	k = ((double)kR)/bmode.R;
	d_kappa = (dir == FORW) ? k : -k;
	s0 = 0.0;
	sN = 2.0*PI;
	ds_hom = d_p;
	ds_disc = d_p;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
	u.strip();
	Ns -= 1;
	u = Cvector(Ns+1);
	u.init(CC0);
}
HcmtElement::HcmtElement(HcmtBft t, BDMode m, 
                         double x_0, double z_0, double d_p)
{
	switch(t)
	{
	case CCC: case CCA:
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, circular cavity");
	}
	int kR;
	kR = (int)(floor(m.beta*m.R+0.5));
	bmode = m;
	vM = 0;
	ky = 0.0;
	pol = bmode.pol;	
	k0 = bmode.k0;
	invommu0 = val_invommu0(bmode.wg.lambda);
	invomep0 = val_invomep0(bmode.wg.lambda);
	x0 = x_0;
	z0 = z_0;
	dir = FORW;
	if(type == CCA) dir = BACK;
	d_gamma = (dir == FORW) ? bmode.c_gamma() : bmode.c_gamma()*(-1.0);
	d_beta =  (dir == FORW) ? bmode.beta  : -bmode.beta;
	d_alpha = (dir == FORW) ? bmode.alpha : -bmode.alpha;
	double k;
	k = ((double)kR)/bmode.R;
	d_kappa = (dir == FORW) ? k : -k;
	s0 = 0.0;
	sN = 2.0*PI;
	ds_hom = d_p;
	ds_disc = d_p;
	Ns = -1;
	nseg.strip();
	seg = NULL;
	initFEMreg();
	u.strip();
	Ns -= 1;
	u = Cvector(Ns+1);
	u.init(CC0);
}

/* initialization, whispering gallery mode supported by a
   circular cavity
   t:   WGM 
   wgm: the corresponding whispering gallery mode object */
HcmtElement::HcmtElement(HcmtBft t, WGMode wgm)
{
	switch(t)
	{
	case WGM: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, whispering gallery mode");
	}

	wgmode = wgm;
	vM = 0;
	ky = 0.0;
	pol = wgm.pol;	
	k0 = val_k0(val_lambda(wgm.omega));
	invommu0 = val_invommu0(val_lambda(wgm.omega));
	invomep0 = val_invomep0(val_lambda(wgm.omega));
	x0 = wgm.x0;
	z0 = wgm.z0;
	u.strip();
	seg = NULL;
	Ns = 0;
	u = Cvector(Ns+1);
	u.init(CC0);
}

/* initialization, resonant field with complex frequency omega [rad/fs],
   discretized on a rectangular FD grid
   x_0 -> (numx) -> x_n, z_0 -> (numz) -> z_n, numx, numz: #intervals,
   file input, t: MAT */
HcmtElement::HcmtElement(HcmtBft t, Polarization p, Complex omega,
                         const char *fname, 
                         double x_0, double x_n, int num_x,
                         double z_0, double z_n, int num_z)
{
	if(t!=MAT)
		cmterror("HcmtElement constructor, discretized field");
	type = t;
	if(p!=TE)
		cmterror("HcmtElement constructor, MAT, TM not implemented");
	vM = 0;
	ky = 0.0;
	pol = p;
	double l = 1.883651567/omega.re;
	k0 = val_k0(l);
	invommu0 = val_invommu0(l);
	invomep0 = val_invomep0(l);
	x0 = x_0; xn = x_n; numx = num_x; dx = (xn-x0)/((double)numx);
	z0 = z_0; zn = z_n; numz = num_z; dz = (zn-z0)/((double)numz);
	FDfEy = Cmatrix(numx+1, numz+1);
	FDfHx = Cmatrix(numx+1, numz+1);
	FDfHz = Cmatrix(numx+1, numz+1);
	FDfdzEy = Cmatrix(numx+1, numz+1);
	FDfdxEy = Cmatrix(numx+1, numz+1);
	FDfdzHxmdxHz = Cmatrix(numx+1, numz+1);
	FDfEy.init(CC0);
	FDfHx.init(CC0);
	FDfHz.init(CC0);
	FDfdzEy.init(CC0);
	FDfdxEy.init(CC0);
	FDfdzHxmdxHz.init(CC0);
	FILE *f;
        fprintf(stderr, "HCMT_EL(MAT) << %s", fname);
        f = fopen(fname, "r");
	Complex c;
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfEy(i, j) = c;
		}
	}
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfHx(i, j) = c;
		}
	}
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfHz(i, j) = c;
		}
	}
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfdzEy(i, j) = c;
		}
	}
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfdxEy(i, j) = c;
		}
	}
	for(int i=0; i<=numx; ++i)
	{
		for(int j=0; j<=numz; ++j)
		{
			c.re = fgetdouble(f); c.im = fgetdouble(f);
			FDfdzHxmdxHz(i, j) = c;
		}
	}
	fclose(f);
        fprintf(stderr, ".\n");

// reciprocal time dependence (!)
	FDfEy.conj();
	FDfHx.conj();
	FDfHz.conj();
	FDfdzEy.conj();
	FDfdxEy.conj();
	FDfdzHxmdxHz.conj();

	u.strip();
	seg = NULL;
	Ns = 0;
	u = Cvector(Ns+1);
	u.init(CC0);
}

/* FEM discretization of the amplitude function: initialization, regular */
void HcmtElement::initFEMreg()
{
	if(s0 >= sN) cmterror("HcmtElement::initFEMreg: s0, sN");
	Ns = (int) floor((sN-s0)/ds_hom+0.5);
	if(Ns < 1) Ns = 1;
	ds_hom = (sN-s0)/((double) Ns);
	ds_disc = ds_hom;
	sbeg = Dvector(Ns+1);
	scen = Dvector(Ns+1);
	send = Dvector(Ns+1);
	sdsL = Dvector(Ns+1);
	sdsR = Dvector(Ns+1);
        for(int k=0; k<=Ns; ++k)
        {
                scen(k) = s0+((double) k)*ds_hom;
                sbeg(k) = scen(k)-ds_hom;
                send(k) = scen(k)+ds_hom;
		sdsL(k) = ds_hom;
		sdsR(k) = ds_hom;
        }
	u = Cvector(Ns+1);
	u.init(CC0);
	return;
}

/* FEM discretization of the amplitude function: initialization,
   irregular, structure-adapted */
void HcmtElement::initFEMstr(Dvector d)
{
	if(s0 >= sN) cmterror("HcmtElement::initFEMstr: s0, sN");
	if(ds_hom <= ds_disc) cmterror("HcmtElement::initFEMstr: d_s");
	int maxnp = 10000;
	Dvector p(maxnp+2);
	double nd, dsd;
	dsd = ds_disc;
	for(int i=1; i<=d.nel-1; ++i)
	{
		if(dsd > (d(i)-d(i-1))/3.0) dsd = (d(i)-d(i-1))/3.0;
	}
	int np = 0;
	p(np) = s0; ++np; 
	for(int i=0; i<=d.nel-1; ++i)
	{
		nd = d(i)-dsd;
		if(nd > s0+dsd-HDIST && nd < sN-dsd+HDIST)
		{
			p(np) = nd; ++np; 
			if(np > maxnp) 
				cmterror("HcmtElement::initFEMstr: maxnp");
		}
		nd = d(i);
		if(nd > s0+dsd-HDIST && nd < sN-dsd+HDIST)
		{
			p(np) = nd; ++np; 
			if(np > maxnp) 
				cmterror("HcmtElement::initFEMstr: maxnp");
		}
		nd = d(i)+dsd;
		if(nd > s0+dsd-HDIST && nd < sN-dsd+HDIST)
		{
			p(np) = nd; ++np; 
			if(np > maxnp) 
				cmterror("HcmtElement::initFEMstr: maxnp");
		}
	}
	p(np) = sN; ++np; 
	Dvector fp(maxnp+2);
	int nfp = 0;
	fp(nfp) = p(0); ++nfp;
	for(int i=1; i<=np-1; ++i)
	{
		if(p(i) - fp(nfp-1) > ds_hom)
		{
			double s, sr, ds;
			int ns;
			sr = p(i);	
			s  = fp(nfp-1);	
			ns = (int) ceil((sr-s)/ds_hom);
			if(ns < 1) ns = 1;
			ds = (sr-s)/((double) ns);
			for(int k=1; k<=ns; ++k)
			{
				fp(nfp) = s+((double) k)*ds; ++nfp; 
				if(nfp > maxnp) 
				   cmterror("HcmtElement::initFEMstr: maxnp");
			}
		}
		else 
		{
			fp(nfp) = p(i); ++nfp;
			if(nfp > maxnp) 
				cmterror("HcmtElement::initFEMstr: maxnp");
		}
	}
	Ns = nfp-1;
	sbeg = Dvector(Ns+1);
	scen = Dvector(Ns+1);
	send = Dvector(Ns+1);
	sdsL = Dvector(Ns+1);
	sdsR = Dvector(Ns+1);
	for(int k=0; k<=Ns; ++k)
	{
		scen(k) = fp(k);
		if(k > 0) sbeg(k)  = fp(k-1);
		else      sbeg(k)  = -AWAY;
		if(k < Ns) send(k) = fp(k+1);
		else       send(k) = AWAY;
		sdsL(k) = scen(k)-sbeg(k);
		sdsR(k) = send(k)-scen(k);
	}
// fprintf(stderr, "\n # = %d\n", Ns+1);
// for(int k=0; k<=Ns; ++k)
// fprintf(stderr, "(%d) [ %g : %g : %g : %g : %g ]\n", 
// k, sbeg(k), sdsL(k), scen(k), sdsR(k), send(k));
	u = Cvector(Ns+1);
	u.init(CC0);
	return;
}

/* linear 1-D elements, triangle functions */ 
double HcmtElement::elf(int k, double s) const
{
	if(k>=1 && k<=Ns-1)
	{
		if(s <= sbeg(k)) return 0.0;
		if(s >= send(k)) return 0.0;
		if(s <= scen(k)) return (s-sbeg(k))/sdsL(k);
		return (send(k)-s)/sdsR(k);
	}
	if(k==0)
	{
		if(s <= scen(0)) return 1.0;
		if(s >= send(0)) return 0.0;
		return (send(0)-s)/sdsR(k);
	}
	if(k==Ns)
	{
		if(s >= scen(Ns)) return 1.0;
		if(s <= sbeg(Ns)) return 0.0;
		return (s-sbeg(Ns))/sdsL(k);
	}
	return 0.0;
}

/* ... and their derivatives */ 
double HcmtElement::d_elf(int k, double s) const
{
	if(k>=1 && k<=Ns-1)
	{
		if(s <= sbeg(k)) return 0.0;
		if(s >= send(k)) return 0.0;
		if(s <= scen(k)) return 1.0/sdsL(k);
		return -1.0/sdsR(k);
	}
	if(k==0)
	{
		if(s <= scen(0)) return 0.0;
		if(s >= send(0)) return 0.0;
		return -1.0/sdsR(k);
	}
	if(k==Ns)
	{
		if(s >= scen(Ns)) return 0.0;
		if(s <= sbeg(Ns)) return 0.0;
		return 1.0/sdsL(k);
	}
	return 0.0;
}


/* linear 1-D elements, triangle functions, end elements are identified */ 
double HcmtElement::celf(int k, double s) const
{
	if(k>=1 && k<=Ns)
	{
		if(s < sbeg(k)) return 0.0;
		if(s > send(k)) return 0.0;
		if(s <= scen(k)) return (s-sbeg(k))/sdsL(k);
		return (send(k)-s)/sdsR(k);
	}
	if(k==0)
	{
		if(s <  scen(0)) return 0.0;
		if(s <= send(0)) return (send(0)-s)/sdsR(0);
		if(s < sbeg(Ns+1)) return 0.0;
		return (s-sbeg(Ns+1))/sdsL(Ns+1);
	}
	if(k>Ns)
	{
		cmterror("HcmtElement::celf: illegal index");
		return 0.0;
	}
	return 0.0;
}

/* ... and their derivatives */ 
double HcmtElement::d_celf(int k, double s) const
{
	if(k>=1 && k<=Ns)
	{
		if(s <= sbeg(k)) return 0.0;
		if(s >= send(k)) return 0.0;
		if(s <= scen(k)) return 1.0/sdsL(k);
		return -1.0/sdsR(k);
	}
	if(k==0)
	{
		if(s <= scen(0)) return 0.0;
		if(s <= send(0)) return -1.0/sdsR(0);
		if(s <= sbeg(Ns+1)) return 0.0;
		return 1.0/sdsL(Ns+1);
	}
	if(k>Ns)
	{
		cmterror("HcmtElement::d_celf: illegal index");
		return 0.0;
	}
	return 0.0;
}

/* the amplitude function, constant input / output levels */
Complex HcmtElement::amp(double s) const
{
	Complex a = CC0;
	for(int k=0; k<=Ns; ++k) a += u(k)*elf(k, s);	
	return a;
}

/* circular amplitude function, end elements are identified */
Complex HcmtElement::camp(double s) const
{
	Complex a = CC0;
	for(int k=0; k<=Ns; ++k) a += u(k)*celf(k, s);	
	return a;
}

/* amplitude function -> .xyf data file 
   foa: MOD, SQR, REP, IMP
   disp: output window on the element coordinate axis 
   np: number of output points
   ext0, ext1: filename id characters */
void HcmtElement::amp2file(Afo foa, Interval disp, int np, 
                           char ext0, char ext1) const
{
	if(np <= 0) cmterror("amp2file: np <= 0");
	char name[15] = "a______";
	char c0, c1, c2;
	hcmtbftchr(type, c0, c1, c2);
	name[1] = c0;
	name[2] = c1;
	name[3] = c2;
	name[4] = ext0;
	name[5] = ext1;
	name[6] = afochr(foa);
	double sa = disp.x0;
	double sb = disp.x1;

	switch(foa)
	{
	case ORG:
		fprintf(stderr, "   a_%c%c%c%c%c [%g, %g] ", 
		                c0, c1, c2, ext0, ext1, sa, sb);
		break;
	case MOD:
		fprintf(stderr, "   |a_%c%c%c%c%c| [%g, %g] ", 
		                c0, c1, c2, ext0, ext1, sa, sb);
		break;
	case SQR:
		fprintf(stderr, "   |a_%c%c%c%c%c|^2 [%g, %g] ", 
		                c0, c1, c2, ext0, ext1, sa, sb);
		break;
	case REP:
		fprintf(stderr, "   Re a_%c%c%c%c%c [%g, %g] ", 
		                c0, c1, c2, ext0, ext1, sa, sb);
		break;
	case IMP:
		fprintf(stderr, "   Im a_%c%c%c%c%c [%g, %g] ", 
		                c0, c1, c2, ext0, ext1, sa, sb);
		break;
	}

	Dvector f(np);
	Dvector s(np);
	switch(type)
	{
	case CCC: case CCA:
		if(np == 1) f(0) = realize(camp(0.5*(sa+sb)), foa); 
		else 
		{
			double d = (sb-sa)/((double) (np-1));
			for(int j=0; j<=np-1; ++j)  
			{
				s(j) = sa+((double) j)*d;
				f(j) = realize(camp(s(j)), foa);
			}
		}
		break;
	default:
		if(np == 1) f(0) = realize(amp(0.5*(sa+sb)), foa); 
		else 
		{
			double d = (sb-sa)/((double) (np-1));
			for(int j=0; j<=np-1; ++j)  
			{
				s(j) = sa+((double) j)*d;
				f(j) = realize(amp(s(j)), foa);
			}
		}
		break;
	}
	toxyf(name, s, f);
	return;
}
/* as before, with default parameters, absolute value only */
void HcmtElement::amp2file(char ext0, char ext1) const
{
	Interval disp;
	disp.x0 = s0-(sN-s0)/4.0;
	disp.x1 = sN+(sN-s0)/4.0;
	amp2file(MOD, disp, 1000, ext0, ext1);
	return;
}
/* as before, with default parameters, REP, IMP, MOD */
void HcmtElement::amp2files(char ext0, char ext1) const
{
	Interval disp;
	if(type == CCC || type == CCA)
	{
		disp.x0 = s0;
		disp.x1 = sN;
	}
	else
	{
		disp.x0 = s0-(sN-s0)/4.0;
		disp.x1 = sN+(sN-s0)/4.0;
	}
	int np = 1000;
	char name[15] = "a______";
	char c0, c1, c2;
	hcmtbftchr(type, c0, c1, c2);
	name[1] = c0;
	name[2] = c1;
	name[3] = c2;
	name[4] = ext0;
	name[5] = ext1;
	double sa = disp.x0;
	double sb = disp.x1;
	fprintf(stderr, "a_%c%c%c%c%c [%g, %g]:\n", 
		                c0, c1, c2, ext0, ext1, sa, sb);
	Cvector f(np);
	Dvector s(np);
	switch(type)
	{
	case CCC: case CCA:
		if(np == 1) f(0) = camp(0.5*(sa+sb));
		else 
		{
			double d = (sb-sa)/((double) (np-1));
			for(int j=0; j<=np-1; ++j)  
			{
				s(j) = sa+((double) j)*d;
				f(j) = camp(s(j));
			}
		}
		break;
	default:
		if(np == 1) f(0) = amp(0.5*(sa+sb));
		else 
		{
			double d = (sb-sa)/((double) (np-1));
			for(int j=0; j<=np-1; ++j)  
			{
				s(j) = sa+((double) j)*d;
				f(j) = amp(s(j));
			}
		}
		break;
	}
	fprintf(stderr, "   Re: ");
	name[6] = afochr(REP); toxyf(name, s, realize(f, REP));
	fprintf(stderr, "   Im: ");
	name[6] = afochr(IMP); toxyf(name, s, realize(f, IMP));
	fprintf(stderr, "  |.|: ");
	name[6] = afochr(MOD); toxyf(name, s, realize(f, MOD));
	return;
}

/* basis fields, elementwise, for use in the HCMT solver */
Complex HcmtElement::el_field(FKcp cp, int k, double x, double z)
{
	double a, dxa, dza, dra, dsa;
	Complex cdsa;
	Complex e, f, fhx, fhz, fey, fex, fez, fhy;
	Complex fhr, fhp;
	Complex fer, fep;
	Complex drfey, drfhr, drfhp;
	Complex drfhy, drfer, drfep;
	double r, rq, knr, p, cosp, sinp;
	double dxr, dzr, dxp, dzp;
	Complex feyr, dfeyr, ddfeyr, fhyr, dfhyr, ddfhyr;
	double s, xs, zs;
	Complex fEy0, fHx0, fHz0, dzEy0, dxEy0, dzHx0mdxHz0;
	Complex fHy0, fEx0, fEz0, dzHy0, dxHy0, dzEx0mdxEz0;
	double eps;
	double f1, f2;
	Complex fa, fb, fc, fd;
	int i, j;
	double ra, rb, rc, rd;

	if(pol != TE) goto el_field_TM;

	switch(type)
	{
	case MHF: case MHB:
		a = elf(k, z);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_beta*z);
		switch(cp)
		{
		case fcEy:	
			f = mode.cfield(EY, dir, x);
			return a*f*e;
			break;
		case fcHx:
			f = mode.cfield(HX, dir, x);
			return a*f*e;
			break;
		case fcHz:
			f = mode.cfield(HZ, dir, x);
			return a*f*e;
			break;
		case dzEy:
			f = mode.cfield(EY, dir, x);
			dza = d_elf(k, z);
			return (CCmI*d_beta*a+dza)*f*e;
			break;
		case dxEy:
			f = mode.cfield(HZ, dir, x);
			return CCmI*a*f*e/invommu0;
			break;
		case dzHxmdxHz:
			fhx = mode.cfield(HX, dir, x);
			fey = mode.cfield(EY, dir, x);
			dza = d_elf(k, z);
			return (CCmI*d_beta*a+dza)*fhx*e
			       -CCI*invommu0*a*(d_beta*d_beta
			                       -k0*k0*mode.wg.eps(x))
			       *fey*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case MVU: case MVD:
		a = elf(k, x);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_beta*x);
		switch(cp)
		{
		case fcEy:	
			f = mode.vertfield(EY, dir, z);
			return a*f*e;
			break;
		case fcHx:
			f = mode.vertfield(HX, dir, z);
			return a*f*e;
			break;
		case fcHz:
			f = mode.vertfield(HZ, dir, z);
			return a*f*e;
			break;
		case dzEy:
			f = mode.vertfield(HX, dir, z);
			return CCI*a*f*e/invommu0;
			break;
		case dxEy:
			f = mode.vertfield(EY, dir, z);
			dxa = d_elf(k, x);
			return (CCmI*d_beta*a+dxa)*f*e;
			break;
		case dzHxmdxHz:
			fey = mode.vertfield(EY, dir, z);
			fhz = mode.vertfield(HZ, dir, z);
			dxa = d_elf(k, x);
			return CCmI*invommu0*a*(d_beta*d_beta
			                       -k0*k0*mode.wg.eps(z))
			       *fey*e-(CCmI*d_beta*a+dxa)*fhz*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case SRC: 
		dx = x-x0; dz = z-z0;
		rq = dx*dx+dz*dz; r = sqrt(rq);
		a = elf(k, r);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-ordd*Complex(dx, dz).arg());
		knr = kn*r; if(knr <= 0.01) return CC0;
		  feyr = hankel_H2(orda, knr);
		 dfeyr = d_hankel_H2(orda, knr)*kn;
		ddfeyr = (ordd*ordd/rq-kqnq)*feyr-dfeyr/r;
		switch(cp)
		{
		case fcEy:	
			return a*feyr*e;
			break;
		case fcHx:
			return a*CCmI*invommu0
			       *(dfeyr*dz/r-CCI*ordd*feyr*dx/rq)*e;
			break;
		case fcHz:
			return a*CCI*invommu0
			       *(dfeyr*dx/r+CCI*ordd*feyr*dz/rq)*e;
			break;
		case dzEy:
			dra = d_elf(k, r);
			return ((dra*feyr+a*dfeyr)*dz/r
			        -CCI*ordd*a*feyr*dx/rq)*e;
			break;
		case dxEy:
			dra = d_elf(k, r);
			return ((dra*feyr+a*dfeyr)*dx/r
			        +CCI*ordd*a*feyr*dz/rq)*e;
			break;
		case dzHxmdxHz:
			dra = d_elf(k, r);
			return CCmI*invommu0*(dra*dfeyr
			           +a*(ddfeyr+dfeyr/r-ordd*ordd/rq*feyr))*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case GBM: 
		s = sal*(x-x0)+cal*(z-z0);
		a = elf(k, s);
		if(fabs(a) < 1.0e-10) return CC0; 
		dsa = d_elf(k, s);
		xs = cal*(x-x0)-sal*(z-z0);
		zs = s;
		switch(cp)
		{
		case fcEy:	
			fEy0 = gbm_phi(kn, width, xs, zs);
			return a*baB*fEy0;
			break;
		case fcHx: 
			fHx0 = CCmI*invommu0
			       *(gbm_dxphi(kn, width, xs, zs)*(-sal)
				+gbm_dzphi(kn, width, xs, zs)*cal);
			return a*baB*fHx0;
			break;
		case fcHz:
			fHz0 = CCI*invommu0
			       *(gbm_dxphi(kn, width, xs, zs)*cal
				+gbm_dzphi(kn, width, xs, zs)*sal);
			return a*baB*fHz0;
			break;
		case dzEy:
			fEy0 = gbm_phi(kn, width, xs, zs);
			dzEy0 = gbm_dxphi(kn, width, xs, zs)*(-sal)
			       +gbm_dzphi(kn, width, xs, zs)*cal;
			return baB*(dsa*cal*fEy0+a*dzEy0);
			break;
		case dxEy:
			fEy0 = gbm_phi(kn, width, xs, zs);
			dxEy0 = gbm_dxphi(kn, width, xs, zs)*cal
			       +gbm_dzphi(kn, width, xs, zs)*sal;
			return baB*(dsa*sal*fEy0+a*dxEy0);
			break;
		case dzHxmdxHz:
			fEy0 = gbm_phi(kn, width, xs, zs);
			fHx0 = CCmI*invommu0
			       *(gbm_dxphi(kn, width, xs, zs)*(-sal)
				+gbm_dzphi(kn, width, xs, zs)*cal);
			fHz0 = CCI*invommu0
			       *(gbm_dxphi(kn, width, xs, zs)*cal
				+gbm_dzphi(kn, width, xs, zs)*sal);
			return baB*(dsa*(cal*fHx0-sal*fHz0)
				  +a*CCI/invomep0*nb*nb*fEy0);
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case RES: 
		xs = x-x0;
		zs = z-z0;
		switch(cp)
		{
		case fcEy:	
			f =    mode.field(EY, xs)*rmode.field(EY, zs)
			    - rmode.field(EY, xs)* mode.field(EY, zs);
			return f;
			break;
		case fcHx: 
			f =    mode.field(EY, xs)*rmode.field(HZ, zs)
			    - rmode.field(EY, xs)* mode.field(HZ, zs);
			return CCmI*f;
			break;
		case fcHz:
			f =    mode.field(HZ, xs)*rmode.field(EY, zs)
			    - rmode.field(HZ, xs)* mode.field(EY, zs);
			return CCI*f;
			break;
		case dzEy:
			f =    mode.field(EY, xs)*rmode.field(HZ, zs)
			    - rmode.field(EY, xs)* mode.field(HZ, zs);
			return f/invommu0;
			break;
		case dxEy:
			f =    mode.field(HZ, xs)*rmode.field(EY, zs)
			    - rmode.field(HZ, xs)* mode.field(EY, zs);
			return f/invommu0;
			break;
		case dzHxmdxHz:
			f = CC0;
			f1 = mode.field(EY, xs);
				eps = rmode.wg.eps(zs); 
			f2 = (rmode.betaq-k0*k0*eps)*rmode.field(EY, zs);
			f += f1*f2;
			f1 = rmode.field(EY, xs);
				eps = mode.wg.eps(zs); 
			f2 = (mode.betaq-k0*k0*eps)*mode.field(EY, zs);
			f -= f1*f2;
				eps = mode.wg.eps(xs); 
			f1 = (mode.betaq-k0*k0*eps)*mode.field(EY, xs);
			f2 = rmode.field(EY, zs);
			f += f1*f2;
				eps = rmode.wg.eps(xs); 
			f1 = (rmode.betaq-k0*k0*eps)*rmode.field(EY, xs);
			f2 = mode.field(EY, zs);
			f -= f1*f2;
			return CCmI*invommu0*f;
			break;

		default: return CC0;
		}
		return CC0;
		break;

	case CCC: case CCA:
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<BM_R_CROP_IN ) return CC0;
		if(r>BM_R_CROP_OUT) return CC0;
		if(r<1.0e-8) return CC0;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;

		a = celf(k, p);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_kappa*bmode.R*p);
		switch(cp)
		{
		case fcEy:	
			f = bmode.field(EY, r);
			return a*f*e;
			break;
		case fcHx:
			if(dir == FORW)
			f =  bmode.field(HR, r)*cosp-bmode.field(HP, r)*sinp;
			else 
			f = -bmode.field(HR, r)*cosp-bmode.field(HP, r)*sinp;
			return a*f*e;
			break;
		case fcHz:
			if(dir == FORW)
			f =  bmode.field(HR, r)*sinp+bmode.field(HP, r)*cosp;
			else 
			f = -bmode.field(HR, r)*sinp+bmode.field(HP, r)*cosp;
			return a*f*e;
			break;
		case dzEy:
			dzr = dz/r;
			if(    (    PI/4.0<p && p<3.0*PI/4.0) 
			    || (5.0*PI/4.0<p && p<7.0*PI/4.0))
				dzp = cosp*sinp/dz;	
			else    dzp = cosp*cosp/dx;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			fey = bmode.field(EY, r);
			drfey = CCmI*bmode.field(HP, r)/invommu0;
			return (cdsa*dzp*fey + a*drfey*dzr)*e;
			break;
		case dxEy:
			dxr = dx/r;
			if(    (                p<    PI/4.0) 
			    || (3.0*PI/4.0<p && p<5.0*PI/4.0) 
			    || (7.0*PI/4.0<p                ))
				dxp = -cosp*sinp/dx;	
			else    dxp = -sinp*sinp/dz;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			fey = bmode.field(EY, r);
			drfey = CCmI*bmode.field(HP, r)/invommu0;
			return (cdsa*dxp*fey + a*drfey*dxr)*e;
			break;
		case dzHxmdxHz:
			dzr = dz/r;
			dxr = dx/r;
			if(    (    PI/4.0<p && p<3.0*PI/4.0) 
			    || (5.0*PI/4.0<p && p<7.0*PI/4.0))
				dzp = cosp*sinp/dz;	
			else    dzp = cosp*cosp/dx;
			if(    (                p<    PI/4.0) 
			    || (3.0*PI/4.0<p && p<5.0*PI/4.0) 
			    || (7.0*PI/4.0<p                ))
				dxp = -cosp*sinp/dx;	
			else    dxp = -sinp*sinp/dz;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			if(dir == FORW) fhr =  bmode.field(HR, r);
			else            fhr = -bmode.field(HR, r);
			fhp = bmode.field(HP, r);
			fhx = fhr*cosp-fhp*sinp;
			fhz = fhr*sinp+fhp*cosp;
			fey = bmode.field(EY, r);
			drfhr = (CCI*d_gamma*bmode.R*fhp-fhr)/r;
			drfhp = fhp/(-r)
			      - CCI*invommu0*(k0*k0*bmode.wg.eps(r-bmode.R)
			      -d_gamma*d_gamma*bmode.R*bmode.R/r/r)
			      *fey;
			return 
			((dzp*cdsa-dxp*a)*fhx -(dzp*a+dxp*cdsa)*fhz
			                     +a*drfhr*(dzr*cosp-dxr*sinp)
			                     -a*drfhp*(dzr*sinp+dxr*cosp))*e;	
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case WGM: 
		s = 1.0; if(wgmode.dir == BACK) s = -1.0;
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<WM_R_CROP_IN ) return CC0;
		if(r>WM_R_CROP_OUT) return CC0;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;

		e = expi(-((double) wgmode.m)*p*s);
		switch(cp)
		{
		case fcEy:	
			f = wgmode.profile(r);
			return f*e;
			break;
		case fcHx:
			return CCmI*(wgmode.d_profile(r)*sinp 
		                     - CCI*((double) wgmode.m)*s/r*cosp*wgmode.profile(r))*e*wgmode.invommu0();
			break;
		case fcHz:
			return CCI*(wgmode.d_profile(r)*cosp 
		                     + CCI*((double) wgmode.m)*s/r*sinp*wgmode.profile(r))*e*wgmode.invommu0();
			break;
		case dzEy:
			return (wgmode.d_profile(r)*sinp 
		                - CCI*((double) wgmode.m)*s/r*cosp*wgmode.profile(r))*e;
			break;
		case dxEy:
			return (wgmode.d_profile(r)*cosp 
		                + CCI*((double) wgmode.m)*s/r*sinp*wgmode.profile(r))*e;
			break;
		case dzHxmdxHz:
			eps = wgmode.cav.eps(x, z);
			return CCI*wgmode.c_omega()*EP0*eps*wgmode.profile(r)*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case MAT: 
		i = (int)(floor((x-x0)/dx));
		if(i == numx) i=numx-1;
		j = (int)(floor((z-z0)/dz));
		if(j == numz) j=numz-1;
		ra = (x0+((double)(i+1))*dx-x)/dx*(z0+((double)(j+1))*dz-z)/dz;
		rb = (x0+((double)(i+1))*dx-x)/dx*(z-((double)(j  ))*dz-z0)/dz;
		rc = (x-((double)(i  ))*dx-x0)/dx*(z0+((double)(j+1))*dz-z)/dz;
		rd = (x-((double)(i  ))*dx-x0)/dx*(z-((double)(j  ))*dz-z0)/dz;
		switch(cp)
		{
		case fcEy:	
			fa = FDfEy(i,   j  );
			fb = FDfEy(i,   j+1);
			fc = FDfEy(i+1, j  );
			fd = FDfEy(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		case fcHx:
			fa = FDfHx(i,   j  );
			fb = FDfHx(i,   j+1);
			fc = FDfHx(i+1, j  );
			fd = FDfHx(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		case fcHz:
			fa = FDfHz(i,   j  );
			fb = FDfHz(i,   j+1);
			fc = FDfHz(i+1, j  );
			fd = FDfHz(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		case dzEy:
			fa = FDfdzEy(i,   j  );
			fb = FDfdzEy(i,   j+1);
			fc = FDfdzEy(i+1, j  );
			fd = FDfdzEy(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		case dxEy:
			fa = FDfdxEy(i,   j  );
			fb = FDfdxEy(i,   j+1);
			fc = FDfdxEy(i+1, j  );
			fd = FDfdxEy(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		case dzHxmdxHz:
			fa = FDfdzHxmdxHz(i,   j  );
			fb = FDfdzHxmdxHz(i,   j+1);
			fc = FDfdzHxmdxHz(i+1, j  );
			fd = FDfdzHxmdxHz(i+1, j+1);
			return fa*ra+fb*rb+fc*rc+fd*rd;
			break;
		default: return CC0;
		}
		return CC0;
		break;
	case FIM: 
		cmterror("HcmtElement::el_field, TE: FIM not implemented");
		return CC0;
		break;
	
	}
	return CC0;

	el_field_TM: ;;

	switch(type)
	{
	case MHF: case MHB:
		a = elf(k, z);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_beta*z);
		switch(cp)
		{
		case fcHy:	
			f = mode.cfield(HY, dir, x);
			return a*f*e;
			break;
		case fcEx:
			f = mode.cfield(EX, dir, x);
			return a*f*e;
			break;
		case fcEz:
			f = mode.cfield(EZ, dir, x);
			return a*f*e;
			break;
		case dzHy:
			f = mode.cfield(HY, dir, x);
			dza = d_elf(k, z);
			return (CCmI*d_beta*a+dza)*f*e;
			break;
		case dxHy:
			f = mode.cfield(EZ, dir, x);
			return CCI*a*f*e/invomep0*mode.wg.eps(x);
			break;
		case dzExmdxEz:
			fex = mode.cfield(EX, dir, x);
			fhy = mode.cfield(HY, dir, x);
			dza = d_elf(k, z);
			eps = mode.wg.eps(x);
			return (CCmI*d_beta*a+dza)*fex*e
			       +CCI*invomep0*a/eps
			        *(d_beta*d_beta-k0*k0*eps)*fhy*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case MVU: case MVD:
		a = elf(k, x);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_beta*x);
		switch(cp)
		{
		case fcHy:	
			f = mode.vertfield(HY, dir, z);
			return a*f*e;
			break;
		case fcEx:
			f = mode.vertfield(EX, dir, z);
			return a*f*e;
			break;
		case fcEz:
			f = mode.vertfield(EZ, dir, z);
			return a*f*e;
			break;
		case dzHy:
			f = mode.vertfield(EX, dir, z);
			return CCmI*a*f*e/invomep0*mode.wg.eps(z);
			break;
		case dxHy:
			f = mode.vertfield(HY, dir, z);
			dxa = d_elf(k, x);
			return (CCmI*d_beta*a+dxa)*f*e;
			break;
		case dzExmdxEz:
			fhy = mode.vertfield(HY, dir, z);
			fez = mode.vertfield(EZ, dir, z);
			dxa = d_elf(k, x);
			eps = mode.wg.eps(z);
			return CCI*invomep0*a/eps
			       *(d_beta*d_beta-k0*k0*eps)
			       *fhy*e-(CCmI*d_beta*a+dxa)*fez*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case SRC: 
		cmterror("HcmtElement::el_field: SRC, TM not yet implemented");
		return CC0;
		break;

	case GBM: 
		s = sal*(x-x0)+cal*(z-z0);
		a = elf(k, s);
		if(fabs(a) < 1.0e-10) return CC0; 
		dsa = d_elf(k, s);
		xs = cal*(x-x0)-sal*(z-z0);
		zs = s;
		switch(cp)
		{
		case fcHy:	
			fHy0 = gbm_phi(kn, width, xs, zs);
			return a*baB*fHy0;
			break;
		case fcEx: 
			fEx0 = CCI*invomep0/(nb*nb)
			       *(gbm_dxphi(kn, width, xs, zs)*(-sal)
				+gbm_dzphi(kn, width, xs, zs)*cal);
			return a*baB*fEx0;
			break;
		case fcEz:
			fEz0 = CCmI*invomep0/(nb*nb)
			       *(gbm_dxphi(kn, width, xs, zs)*cal
				+gbm_dzphi(kn, width, xs, zs)*sal);
			return a*baB*fEz0;
			break;
		case dzHy:
			fHy0 = gbm_phi(kn, width, xs, zs);
			dzHy0 = gbm_dxphi(kn, width, xs, zs)*(-sal)
			       +gbm_dzphi(kn, width, xs, zs)*cal;
			return baB*(dsa*cal*fHy0+a*dzHy0);
			break;
		case dxHy:
			fHy0 = gbm_phi(kn, width, xs, zs);
			dxHy0 = gbm_dxphi(kn, width, xs, zs)*cal
			       +gbm_dzphi(kn, width, xs, zs)*sal;
			return baB*(dsa*sal*fHy0+a*dxHy0);
			break;
		case dzExmdxEz:
			fHy0 = gbm_phi(kn, width, xs, zs);
			fEx0 = CCI*invomep0/(nb*nb)
			       *(gbm_dxphi(kn, width, xs, zs)*(-sal)
				+gbm_dzphi(kn, width, xs, zs)*cal);
			fEz0 = CCmI*invomep0/(nb*nb)
			       *(gbm_dxphi(kn, width, xs, zs)*cal
				+gbm_dzphi(kn, width, xs, zs)*sal);
			return baB*(dsa*(cal*fEx0-sal*fEz0)
				  +a*CCmI/invommu0*fHy0);
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case RES: 
		cmterror("HcmtElement::el_field: RES, TM not yet implemented");
		return CC0;
		break;

	case CCC: case CCA:
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<BM_R_CROP_IN ) return CC0;
		if(r>BM_R_CROP_OUT) return CC0;
		if(r<1.0e-8) return CC0;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;

		a = celf(k, p);
		if(fabs(a) < 1.0e-10) return CC0; 
		e = expi(-d_kappa*bmode.R*p);
		switch(cp)
		{
		case fcHy:	
			f = bmode.field(HY, r);
			return a*f*e;
			break;
		case fcEx:
			if(dir == FORW)
			f =  bmode.field(ER, r)*cosp-bmode.field(EP, r)*sinp;
			else 
			f = -bmode.field(ER, r)*cosp-bmode.field(EP, r)*sinp;
			return a*f*e;
			break;
		case fcEz:
			if(dir == FORW)
			f =  bmode.field(ER, r)*sinp+bmode.field(EP, r)*cosp;
			else 
			f = -bmode.field(ER, r)*sinp+bmode.field(EP, r)*cosp;
			return a*f*e;
			break;
		case dzHy:
			dzr = dz/r;
			if(    (    PI/4.0<p && p<3.0*PI/4.0) 
			    || (5.0*PI/4.0<p && p<7.0*PI/4.0))
				dzp = cosp*sinp/dz;	
			else    dzp = cosp*cosp/dx;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			fhy = bmode.field(HY, r);
			drfhy = CCI*bmode.field(EP, r)*bmode.wg.eps(r-bmode.R)/invomep0;
			return (cdsa*dzp*fhy + a*drfhy*dzr)*e;
			break;
		case dxHy:
			dxr = dx/r;
			if(    (                p<    PI/4.0) 
			    || (3.0*PI/4.0<p && p<5.0*PI/4.0) 
			    || (7.0*PI/4.0<p                ))
				dxp = -cosp*sinp/dx;	
			else    dxp = -sinp*sinp/dz;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			fhy = bmode.field(HY, r);
			drfhy = CCI*bmode.field(EP, r)*bmode.wg.eps(r-bmode.R)/invomep0;
			return (cdsa*dxp*fhy + a*drfhy*dxr)*e;
			break;
		case dzExmdxEz:
			dzr = dz/r;
			dxr = dx/r;
			if(    (    PI/4.0<p && p<3.0*PI/4.0) 
			    || (5.0*PI/4.0<p && p<7.0*PI/4.0))
				dzp = cosp*sinp/dz;	
			else    dzp = cosp*cosp/dx;
			if(    (                p<    PI/4.0) 
			    || (3.0*PI/4.0<p && p<5.0*PI/4.0) 
			    || (7.0*PI/4.0<p                ))
				dxp = -cosp*sinp/dx;	
			else    dxp = -sinp*sinp/dz;
			dsa = d_celf(k, p);
			cdsa = CCmI*d_kappa*bmode.R*a+dsa;
			if(dir == FORW) fer =  bmode.field(ER, r);
			else            fer = -bmode.field(ER, r);
			fep = bmode.field(EP, r);
			fex = fer*cosp-fep*sinp;
			fez = fer*sinp+fep*cosp;
			fhy = bmode.field(HY, r);
			drfer = (CCI*d_gamma*bmode.R*fep-fer)/r;
			drfep = fep/(-r)
			      + CCI*invomep0*(k0*k0
			      -d_gamma*d_gamma*bmode.R*bmode.R/r/r/bmode.wg.eps(r-bmode.R))
			      *fhy;
			return 
			((dzp*cdsa-dxp*a)*fex -(dzp*a+dxp*cdsa)*fez
			                     +a*drfer*(dzr*cosp-dxr*sinp)
			                     -a*drfep*(dzr*sinp+dxr*cosp))*e;	
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case WGM:
		s = 1.0; if(wgmode.dir == BACK) s = -1.0;
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<WM_R_CROP_IN ) return CC0;
		if(r>WM_R_CROP_OUT) return CC0;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;

		e = expi(-((double) wgmode.m)*p*s);
		switch(cp)
		{
		case fcHy:	
			f = wgmode.profile(r);
			return f*e;
			break;
		case fcEx:
			eps = wgmode.cav.eps(x, z);
			return CCI*(wgmode.d_profile(r)*sinp 
		                     - CCI*((double) wgmode.m)*s/r*cosp*wgmode.profile(r))*e*wgmode.invomep0()/eps;
			break;
		case fcEz:
			eps = wgmode.cav.eps(x, z);
			return CCmI*(wgmode.d_profile(r)*cosp 
		                     + CCI*((double) wgmode.m)*s/r*sinp*wgmode.profile(r))*e*wgmode.invomep0()/eps;
			break;
		case dzHy:
			return (wgmode.d_profile(r)*sinp 
		                - CCI*((double) wgmode.m)*s/r*cosp*wgmode.profile(r))*e;
			break;
		case dxHy:
			return (wgmode.d_profile(r)*cosp 
		                + CCI*((double) wgmode.m)*s/r*sinp*wgmode.profile(r))*e;
			break;
		case dzExmdxEz:
			return CCmI*wgmode.c_omega()*MU0*wgmode.profile(r)*e;
			break;
		default: return CC0;
		}
		return CC0;
		break;

	case MAT: 
		cmterror("HcmtElement::el_field: MAT, TM not yet implemented");
		return CC0;
		break;
	case FIM: 
		cmterror("HcmtElement::el_field, TM: FIM not implemented");
		return CC0;
		break;
	}
	return CC0;
}

/* basis fields, elementwise, for use in the vectorial HCMT solver */
Complex HcmtElement::el_field_vec(C3cp cp, int k, double x, double z)
{
	double a, dxa, dza;
	Complex e, f, fhx, fhy, fhz, fex, fey, fez;
	double r, dx, dz;

	if(vM != 1) { cmterror("HcmtElement::el_field_vec: vM?"); return CC0; }

	switch(type)
	{
	case MHF: case MHB:
		a = elf(k, z);
		if(fabs(a) < 1.0e-10) return CC0; 
		// dza = d_elf(k, z);
		dza = d_elf(k, z)+d_alpha*a;
		// e = expi(-d_kz*(z-v_z0));
		e = expi(-d_beta*(z-v_z0));
		switch(cp)
		{
			case f3Ex:	
				f = mode.cfield(EX, dir, x);
				return a*f*e;
				break;
			case f3Ey:	
				f = mode.cfield(EY, dir, x);
				return a*f*e;
				break;
			case f3Ez:	
				f = mode.cfield(EZ, dir, x);
				return a*f*e;
				break;
			case f3Hx:
				f = mode.cfield(HX, dir, x);
				return a*f*e;
				break;
			case f3Hy:
				f = mode.cfield(HY, dir, x);
				return a*f*e;
				break;
			case f3Hz:
				f = mode.cfield(HZ, dir, x);
				return a*f*e;
				break;
			case c3Ex:
				fhx = mode.cfield(HX, dir, x);
				fey = mode.cfield(EY, dir, x);
				return (CCmI/invommu0*a*fhx-dza*fey)*e; 
				break;
			case c3Ey:
				fhy = mode.cfield(HY, dir, x);
				fex = mode.cfield(EX, dir, x);
				return (CCmI/invommu0*a*fhy+dza*fex)*e; 
				break;
			case c3Ez:
				fhz = mode.cfield(HZ, dir, x);
				return CCmI/invommu0*a*fhz*e; 
				break;
			case c3Hx:
				fex = mode.cfield(EX, dir, x);
				fhy = mode.cfield(HY, dir, x);
				return (CCI/invomep0*mode.wg.eps(x)*a*fex-dza*fhy)*e; 
				break;
			case c3Hy:
				fey = mode.cfield(EY, dir, x);
				fhx = mode.cfield(HX, dir, x);
				return (CCI/invomep0*mode.wg.eps(x)*a*fey+dza*fhx)*e; 
				break;
			case c3Hz:
				fez = mode.cfield(EZ, dir, x);
				return CCI/invomep0*mode.wg.eps(x)*a*fez*e; 
				break;
			default: return CC0;
		}
		return CC0;
		break;

	case MVU: case MVD:
		a = elf(k, x);
		if(fabs(a) < 1.0e-10) return CC0; 
		// dxa = d_elf(k, x);
		dxa = d_elf(k, x)+d_alpha*a;
		// e = expi(-d_kz*(x-v_z0));
		e = expi(-d_beta*(x-v_z0));
		switch(cp)
		{
			case f3Ex:	
				f = mode.vertfield(EX, dir, z);
				return a*f*e;
				break;
			case f3Ey:	
				f = mode.vertfield(EY, dir, z);
				return a*f*e;
				break;
			case f3Ez:	
				f = mode.vertfield(EZ, dir, z);
				return a*f*e;
				break;
			case f3Hx:
				f = mode.vertfield(HX, dir, z);
				return a*f*e;
				break;
			case f3Hy:
				f = mode.vertfield(HY, dir, z);
				return a*f*e;
				break;
			case f3Hz:
				f = mode.vertfield(HZ, dir, z);
				return a*f*e;
				break;
			case c3Ex:
				fhx = mode.vertfield(HX, dir, z);
				return CCmI/invommu0*a*fhx*e; 
				break;
			case c3Ey:
				fhy = mode.vertfield(HY, dir, z);
				fez = mode.vertfield(EZ, dir, z);
				return (CCmI/invommu0*a*fhy-dxa*fez)*e; 
				break;
			case c3Ez:
				fhz = mode.vertfield(HZ, dir, z);
				fey = mode.vertfield(EY, dir, z);
				return (CCmI/invommu0*a*fhz+dxa*fey)*e; 
				break;
			case c3Hx:
				fex = mode.vertfield(EX, dir, z);
				return CCI/invomep0*mode.wg.eps(z)*a*fex*e; 
				break;
			case c3Hy:
				fey = mode.vertfield(EY, dir, z);
				fhz = mode.vertfield(HZ, dir, z);
				return (CCI/invomep0*mode.wg.eps(z)*a*fey-dxa*fhz)*e; 
				break;
			case c3Hz:
				fez = mode.vertfield(EZ, dir, z);
				fhy = mode.vertfield(HY, dir, z);
				return (CCI/invomep0*mode.wg.eps(z)*a*fez+dxa*fhy)*e; 
				break;
			default: return CC0;
		}
		return CC0;
		break;

	case FIM: 
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<fimode.cRi) return CC0;
		if(r>fimode.cRe) return CC0;
		switch(cp)
		{
			case f3Ex: return  fimode.field(EX, x, z); break;
			case f3Ey: return  fimode.field(EZ, x, z); break;
			case f3Ez: return  fimode.field(EY, x, z); break;
			case f3Hx: return -fimode.field(HX, x, z); break;
			case f3Hy: return -fimode.field(HZ, x, z); break;
			case f3Hz: return -fimode.field(HY, x, z); break;
			case c3Ex: return CCI/invommu0*fimode.field(HX, x, z)-CCI*(ky-fimode.beta)*fimode.field(EY, x, z); break;
			case c3Ey: return CCI/invommu0*fimode.field(HZ, x, z);                                             break;
			case c3Ez: return CCI/invommu0*fimode.field(HY, x, z)+CCI*(ky-fimode.beta)*fimode.field(EX, x, z); break;
			case c3Hx: return CCI/invomep0*fimode.fib.eps(x, z)*fimode.field(EX, x, z)+CCI*(ky-fimode.beta)*fimode.field(HY, x, z); break;
			case c3Hy: return CCI/invomep0*fimode.fib.eps(x, z)*fimode.field(EZ, x, z);                                             break;
			case c3Hz: return CCI/invomep0*fimode.fib.eps(x, z)*fimode.field(EY, x, z)-CCI*(ky-fimode.beta)*fimode.field(HX, x, z); break;
			default: return CC0;
		}
		return CC0;
		break;
	default: 
		cmterror("HcmtElement::el_field_vec: element type not implemented");
		return CC0;
		break;

	}
	return CC0;
}

/* basis fields, elementwise, for use in the HCMT solver, all components */
Cvector HcmtElement::el_field(int k, double x, double z)
{
	Cvector f(12);
	f.init(CC0);
	double a, dza, dxa, dra, dsa;
	Complex e, fey, fhx, fhz, fhy, fex, fez;
	Complex fhr, fhp, drfey, drfhr, drfhp, cdsa;
	Complex fer, fep, drfhy, drfer, drfep;
	double dzr, dzp, dxr, dxp;
	double r, knr, rq;
	Complex feyr, dfeyr, ddfeyr, fhyr, dfhyr, ddfhyr;
	double aB, dsaB;
	double s, xs, zs; 
	double cosp, sinp, p;
	Complex phi, dxphi, dzphi;
	Complex fEy0, fHx0, fHz0, dzEy0, dxEy0, dzHx0mdxHz0;
	Complex fHy0, fEx0, fEz0, dzHy0, dxHy0, dzEx0mdxEz0;
	Complex ft, mfEx, mfEz, mfHx, mfHz, rfEx, rfEz, rfHx, rfHz;
	double eps, R;
	Complex fa, fb, fc, fd;
	int i, j;
	double ra, rb, rc, rd;

	if(pol != TE) goto el_field_ac_TM;
	
	switch(type)
	{
	case MHF: case MHB:
		e = expi(-d_beta*z);
		a = elf(k, z);
		dza = d_elf(k, z);
		fey = mode.cfield(EY, dir, x);
		fhx = mode.cfield(HX, dir, x);
		fhz = mode.cfield(HZ, dir, x);
		f(fcEy) = a*fey*e;
		f(fcHx) = a*fhx*e;
		f(fcHz) = a*fhz*e;
		f(dzEy) = (CCmI*d_beta*a+dza)*fey*e;
		f(dxEy) = CCmI*a*fhz*e/invommu0;
		f(dzHxmdxHz) = (CCmI*d_beta*a+dza)*fhx*e
			      -CCI*invommu0*a*(d_beta*d_beta
			      -k0*k0*mode.wg.eps(x))*fey*e;
		return f;
		break;

	case MVU: case MVD:
		e = expi(-d_beta*x);
		a = elf(k, x);
		dxa = d_elf(k, x);
		fey = mode.vertfield(EY, dir, z);
		fhx = mode.vertfield(HX, dir, z);
		fhz = mode.vertfield(HZ, dir, z);
		f(fcEy) = a*fey*e;
		f(fcHx) = a*fhx*e;
		f(fcHz) = a*fhz*e;
		f(dzEy) = CCI*a*fhx*e/invommu0;
		f(dxEy) = (CCmI*d_beta*a+dxa)*fey*e;
		f(dzHxmdxHz) = CCmI*invommu0*a*(d_beta*d_beta
		               -k0*k0*mode.wg.eps(z))*fey*e
			       -(CCmI*d_beta*a+dxa)*fhz*e;
		return f;
		break;

	case SRC: 
		dx = x-x0; dz = z-z0;
		rq = dx*dx+dz*dz; r = sqrt(rq);
		knr = kn*r; if(knr <= 0.01) return f;
		e = expi(-ordd*Complex(dx, dz).arg());
		  feyr = hankel_H2(orda, knr);
		 dfeyr = d_hankel_H2(orda, knr)*kn;
		ddfeyr = (ordd*ordd/rq-kqnq)*feyr-dfeyr/r;
		a = elf(k, r);
		dra = d_elf(k, r);
		f(fcEy) = a*feyr*e;
		f(fcHx) = a*CCmI*invommu0*(dfeyr*dz/r-CCI*ordd*feyr*dx/rq)*e;
		f(fcHz) = a*CCI*invommu0*(dfeyr*dx/r+CCI*ordd*feyr*dz/rq)*e;
		f(dzEy) = ((dra*feyr+a*dfeyr)*dz/r-CCI*ordd*a*feyr*dx/rq)*e;
		f(dxEy) = ((dra*feyr+a*dfeyr)*dx/r+CCI*ordd*a*feyr*dz/rq)*e;
		f(dzHxmdxHz) = CCmI*invommu0*(dra*dfeyr
			     +a*(ddfeyr+dfeyr/r-ordd*ordd/rq*feyr))*e;
		return f;
		break;

	case GBM: 
		s = sal*(x-x0)+cal*(z-z0);
		aB = elf(k, s)*baB;
		dsaB = d_elf(k, s)*baB;
		xs = cal*(x-x0)-sal*(z-z0);
		zs = s;
		phi = gbm_phi(kn, width, xs, zs);
		dxphi = gbm_dxphi(kn, width, xs, zs);
		dzphi = gbm_dzphi(kn, width, xs, zs);
		fEy0 = phi;
		dzEy0 = dxphi*(-sal)+dzphi*cal;
		dxEy0 = dxphi*cal+dzphi*sal;
		fHx0 = CCmI*invommu0*dzEy0;
		fHz0 = CCI*invommu0*dxEy0;
		f(fcEy) = aB*fEy0;
		f(fcHx) = aB*fHx0;
		f(fcHz) = aB*fHz0;
		f(dzEy) = dsaB*cal*fEy0+aB*dzEy0;
		f(dxEy) = dsaB*sal*fEy0+aB*dxEy0;
		f(dzHxmdxHz) = dsaB*(cal*fHx0-sal*fHz0)
		               +aB*CCI/invomep0*nb*nb*fEy0;
		return f;
		break;

	case RES: 
		xs = x-x0;
		zs = z-z0;
		mfEx =  mode.field(EY, xs);
		mfEz =  mode.field(EY, zs);
		mfHx =  mode.field(HZ, xs);
		mfHz =  mode.field(HZ, zs);
		rfEx = rmode.field(EY, xs);
		rfEz = rmode.field(EY, zs);
		rfHx = rmode.field(HZ, xs);
		rfHz = rmode.field(HZ, zs);
		f(fcEy) = mfEx*rfEz-rfEx*mfEz;
		f(fcHx) = CCmI*(mfEx*rfHz-rfEx*mfHz);
		f(fcHz) = CCI*(mfHx*rfEz-rfHx*mfEz);
		f(dzEy) = (mfEx*rfHz-rfEx*mfHz)/invommu0;
		f(dxEy) = (mfHx*rfEz-rfHx*mfEz)/invommu0;
		ft = CC0;
		eps = rmode.wg.eps(zs); 
		ft += mfEx*(rmode.betaq-k0*k0*eps)*rfEz;
		eps = mode.wg.eps(zs); 
		ft -= rfEx*(mode.betaq-k0*k0*eps)*mfEz;
		eps = mode.wg.eps(xs); 
		ft += (mode.betaq-k0*k0*eps)*mfEx*rfEz;
		eps = rmode.wg.eps(xs); 
		ft -= (rmode.betaq-k0*k0*eps)*rfEx*mfEz;
		f(dzHxmdxHz) = CCmI*invommu0*ft;
		return f;
		break;

	case CCC: case CCA:
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<BM_R_CROP_IN ) return f;
		if(r>BM_R_CROP_OUT) return f;
		if(r<1.0e-8) return f;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;
		R = bmode.R;

		a = celf(k, p);
		if(fabs(a) < 1.0e-10) return f; 
		e = expi(-d_kappa*R*p);
		fey = bmode.field(EY, r);
		if(dir == FORW) fhr =  bmode.field(HR, r);
		else            fhr = -bmode.field(HR, r);
		fhp = bmode.field(HP, r);
		fhx = fhr*cosp-fhp*sinp;
		fhz = fhr*sinp+fhp*cosp;
		dzr = dz/r;
		dxr = dx/r;
		if(    (0.25*PI<p && p<0.75*PI) 
		    || (1.25*PI<p && p<1.75*PI))
			dzp = cosp*sinp/dz;	
		else    dzp = cosp*cosp/dx;
		if(    (             p<0.25*PI) 
		    || (0.75*PI<p && p<1.25*PI) 
		    || (1.75*PI<p             ))
			dxp = -cosp*sinp/dx;	
		else    dxp = -sinp*sinp/dz;
		dsa = d_celf(k, p);
		cdsa = CCmI*d_kappa*R*a+dsa;
		drfey = CCmI*bmode.field(HP, r)/invommu0;
		drfhr = (CCI*d_gamma*R*fhp-fhr)/r;
		drfhp = fhp/(-r)
		      - CCI*invommu0*(k0*k0*bmode.wg.eps(r-bmode.R)
		      -d_gamma*d_gamma*R*R/r/r)
		      *fey;
		f(fcEy) = a*fey*e;
		f(fcHx) = a*fhx*e;
		f(fcHz) = a*fhz*e;
		f(dzEy) = (cdsa*dzp*fey+a*drfey*dzr)*e;
		f(dxEy) = (cdsa*dxp*fey+a*drfey*dxr)*e;
		f(dzHxmdxHz) = ((dzp*cdsa-dxp*a)*fhx-(dzp*a+dxp*cdsa)*fhz
			        +a*drfhr*(dzr*cosp-dxr*sinp)
			        -a*drfhp*(dzr*sinp+dxr*cosp))*e;	
		return f;
		break;

	case WGM:
		s = 1.0; if(wgmode.dir == BACK) s = -1.0;
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<WM_R_CROP_IN ) return f;
		if(r>WM_R_CROP_OUT) return f;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;
		e = expi(-((double) wgmode.m)*p*s);
		  fey = wgmode.profile(r);
		drfey = wgmode.d_profile(r);

		f(fcEy) = fey*e;
		f(dzEy) = (drfey*sinp - CCI*((double) wgmode.m)*s/r*cosp*fey)*e;
		f(fcHx) = CCmI*f(dzEy)*wgmode.invommu0();
		f(dxEy) = (drfey*cosp + CCI*((double) wgmode.m)*s/r*sinp*fey)*e;
		f(fcHz) = CCI*f(dxEy)*wgmode.invommu0();
		f(dzHxmdxHz) = CCI*wgmode.c_omega()*EP0*wgmode.cav.eps(x,z)*fey*e;

/*
if(f(fcEy).abs() > 1.0e-10)
{
	fprintf(stderr, "[%d] (%g, %g)\n", k, x, z);
	double ddz = 1.0e-6;
	double ddx = 1.0e-6;
	dzEy0 = ( el_field(fcEy, k, x, z+ddz/2.0)
	         -el_field(fcEy, k, x, z-ddz/2.0))/ddz;
	dxEy0 = ( el_field(fcEy, k, x+ddx/2.0, z)
	         -el_field(fcEy, k, x-ddx/2.0, z))/ddx;
	dzHx0mdxHz0 = 
	        ( el_field(fcHx, k, x, z+ddz/2.0)
	         -el_field(fcHx, k, x, z-ddz/2.0))/ddz 
	       -( el_field(fcHz, k, x+ddx/2.0, z)
	         -el_field(fcHz, k, x-ddx/2.0, z))/ddx;
	fprintf(stderr, "dzEy:      |%g| +- %g\n", 
	                 f(dzEy).abs(), (f(dzEy)-dzEy0).abs());
	fprintf(stderr, "dxEy:      |%g| +- %g\n", 
	                 f(dxEy).abs(), (f(dxEy)-dxEy0).abs());
	fprintf(stderr, "dzHxmdxHz: |%g| +- %g\n", 
	                 f(dzHxmdxHz).abs(), (f(dzHxmdxHz)-dzHx0mdxHz0).abs());
}
*/
		return f;
		break;

	case MAT: 
		i = (int)(floor((x-x0)/dx));
		if(i == numx) i=numx-1;
		j = (int)(floor((z-z0)/dz));
		if(j == numz) j=numz-1;
		ra = (x0+((double)(i+1))*dx-x)/dx*(z0+((double)(j+1))*dz-z)/dz;
		rb = (x0+((double)(i+1))*dx-x)/dx*(z-((double)(j  ))*dz-z0)/dz;
		rc = (x-((double)(i  ))*dx-x0)/dx*(z0+((double)(j+1))*dz-z)/dz;
		rd = (x-((double)(i  ))*dx-x0)/dx*(z-((double)(j  ))*dz-z0)/dz;
		fa = FDfEy(i,   j  );
		fb = FDfEy(i,   j+1);
		fc = FDfEy(i+1, j  );
		fd = FDfEy(i+1, j+1);
		f(fcEy) = fa*ra+fb*rb+fc*rc+fd*rd;
		fa = FDfHx(i,   j  );
		fb = FDfHx(i,   j+1);
		fc = FDfHx(i+1, j  );
		fd = FDfHx(i+1, j+1);
		f(fcHx) = fa*ra+fb*rb+fc*rc+fd*rd;
		fa = FDfHz(i,   j  );
		fb = FDfHz(i,   j+1);
		fc = FDfHz(i+1, j  );
		fd = FDfHz(i+1, j+1);
		f(fcHz) = fa*ra+fb*rb+fc*rc+fd*rd;
		fa = FDfdzEy(i,   j  );
		fb = FDfdzEy(i,   j+1);
		fc = FDfdzEy(i+1, j  );
		fd = FDfdzEy(i+1, j+1);
		f(dzEy) = fa*ra+fb*rb+fc*rc+fd*rd;
		fa = FDfdxEy(i,   j  );
		fb = FDfdxEy(i,   j+1);
		fc = FDfdxEy(i+1, j  );
		fd = FDfdxEy(i+1, j+1);
		f(dxEy) = fa*ra+fb*rb+fc*rc+fd*rd;
		fa = FDfdzHxmdxHz(i,   j  );
		fb = FDfdzHxmdxHz(i,   j+1);
		fc = FDfdzHxmdxHz(i+1, j  );
		fd = FDfdzHxmdxHz(i+1, j+1);
		f(dzHxmdxHz) = fa*ra+fb*rb+fc*rc+fd*rd;
		return f;
		break;
	case FIM: 
		cmterror("HcmtElement::el_field_ac, TE: FIM not implemented");
		return f;
		break;
	}
	return f;

	el_field_ac_TM: ;;

	switch(type)
	{
	case MHF: case MHB:
		e = expi(-d_beta*z);
		a = elf(k, z);
		dza = d_elf(k, z);
		fhy = mode.cfield(HY, dir, x);
		fex = mode.cfield(EX, dir, x);
		fez = mode.cfield(EZ, dir, x);
		eps = mode.wg.eps(x);
		f(fcHy) = a*fhy*e;
		f(fcEx) = a*fex*e;
		f(fcEz) = a*fez*e;
		f(dzHy) = (CCmI*d_beta*a+dza)*fhy*e;
		f(dxHy) = CCI*a*fez*e*eps/invomep0;
		f(dzExmdxEz) = (CCmI*d_beta*a+dza)*fex*e
			      +CCI*invomep0/eps*a*(d_beta*d_beta
			      -k0*k0*eps)*fhy*e;
		return f;
		break;

	case MVU: case MVD:
		e = expi(-d_beta*x);
		a = elf(k, x);
		dxa = d_elf(k, x);
		fhy = mode.vertfield(HY, dir, z);
		fex = mode.vertfield(EX, dir, z);
		fez = mode.vertfield(EZ, dir, z);
		eps = mode.wg.eps(z);
		f(fcHy) = a*fhy*e;
		f(fcEx) = a*fex*e;
		f(fcEz) = a*fez*e;
		f(dzHy) = CCmI*a*fex*e*eps/invomep0;
		f(dxHy) = (CCmI*d_beta*a+dxa)*fhy*e;
		f(dzExmdxEz) = CCI*invomep0/eps*a*(d_beta*d_beta
		               -k0*k0*eps)*fhy*e
			       -(CCmI*d_beta*a+dxa)*fez*e;
		return f;
		break;

	case SRC: 
	cmterror("HcmtElement::el_field_ac: SRC, TM not yet implemented");
		return f;
		break;

	case GBM: 
		s = sal*(x-x0)+cal*(z-z0);
		aB = elf(k, s)*baB;
		dsaB = d_elf(k, s)*baB;
		xs = cal*(x-x0)-sal*(z-z0);
		zs = s;
		phi = gbm_phi(kn, width, xs, zs);
		dxphi = gbm_dxphi(kn, width, xs, zs);
		dzphi = gbm_dzphi(kn, width, xs, zs);
		fHy0 = phi;
		dzHy0 = dxphi*(-sal)+dzphi*cal;
		dxHy0 = dxphi*cal+dzphi*sal;
		fEx0 = CCI*invomep0/(nb*nb)*dzHy0;
		fEz0 = CCmI*invomep0/(nb*nb)*dxHy0;
		f(fcHy) = aB*fHy0;
		f(fcEx) = aB*fEx0;
		f(fcEz) = aB*fEz0;
		f(dzHy) = dsaB*cal*fHy0+aB*dzHy0;
		f(dxHy) = dsaB*sal*fHy0+aB*dxHy0;
		f(dzExmdxEz) = dsaB*(cal*fEx0-sal*fEz0)
		               +aB*CCmI/invommu0*fHy0;
		return f;
		break;

	case RES: 
	cmterror("HcmtElement::el_field_ac: RES, TM not yet implemented");
		return f;
		break;

	case CCC: case CCA:
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<BM_R_CROP_IN ) return f;
		if(r>BM_R_CROP_OUT) return f;
		if(r<1.0e-8) return f;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;
		R = bmode.R;

		a = celf(k, p);
		if(fabs(a) < 1.0e-10) return f; 
		e = expi(-d_kappa*R*p);
		fhy = bmode.field(HY, r);
		if(dir == FORW) fer =  bmode.field(ER, r);
		else            fer = -bmode.field(ER, r);
		fep = bmode.field(EP, r);
		fex = fer*cosp-fep*sinp;
		fez = fer*sinp+fep*cosp;
		dzr = dz/r;
		dxr = dx/r;
		if(    (0.25*PI<p && p<0.75*PI) 
		    || (1.25*PI<p && p<1.75*PI))
			dzp = cosp*sinp/dz;	
		else    dzp = cosp*cosp/dx;
		if(    (             p<0.25*PI) 
		    || (0.75*PI<p && p<1.25*PI) 
		    || (1.75*PI<p             ))
			dxp = -cosp*sinp/dx;	
		else    dxp = -sinp*sinp/dz;
		dsa = d_celf(k, p);
		cdsa = CCmI*d_kappa*R*a+dsa;
		drfhy = CCI*bmode.field(EP, r)*bmode.wg.eps(r-bmode.R)/invomep0;
		drfer = (CCI*d_gamma*R*fep-fer)/r;
		drfep = fep/(-r)
		      + CCI*invomep0*(k0*k0
		      -d_gamma*d_gamma*R*R/r/r/bmode.wg.eps(r-bmode.R))
		      *fhy;
		f(fcHy) = a*fhy*e;
		f(fcEx) = a*fex*e;
		f(fcEz) = a*fez*e;
		f(dzHy) = (cdsa*dzp*fhy+a*drfhy*dzr)*e;
		f(dxHy) = (cdsa*dxp*fhy+a*drfhy*dxr)*e;
		f(dzExmdxEz) = ((dzp*cdsa-dxp*a)*fex-(dzp*a+dxp*cdsa)*fez
			        +a*drfer*(dzr*cosp-dxr*sinp)
			        -a*drfep*(dzr*sinp+dxr*cosp))*e;	
		return f;
		break;

	case WGM: 
		s = 1.0; if(wgmode.dir == BACK) s = -1.0;
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<WM_R_CROP_IN ) return f;
		if(r>WM_R_CROP_OUT) return f;
		cosp = dx/r;
		sinp = dz/r;
		p = atan2(dz, dx); if(p<0) p+=2.0*PI;
		e = expi(-((double) wgmode.m)*p*s);
		  fhy = wgmode.profile(r);
		drfhy = wgmode.d_profile(r);
		eps = wgmode.cav.eps(x,z);

		f(fcHy) = fhy*e;
		f(dzHy) = (drfhy*sinp - CCI*((double) wgmode.m)*s/r*cosp*fhy)*e;
		f(fcEx) = CCI*f(dzHy)*wgmode.invomep0()/eps;
		f(dxHy) = (drfhy*cosp + CCI*((double) wgmode.m)*s/r*sinp*fhy)*e;
		f(fcEz) = CCmI*f(dxHy)*wgmode.invomep0()/eps;
		f(dzExmdxEz) = CCmI*wgmode.c_omega()*MU0*fhy*e;
		return f;
		break;

	case MAT:
		cmterror("HcmtElement::el_field_ac: MAT, TM not yet implemented");
		return f;
		break;
	case FIM: 
		cmterror("HcmtElement::el_field_ac, TM: FIM not implemented");
		return f;
		break;
	}
	return f;
}

/* basis fields, elementwise, for use in the vectorial HCMT solver, all components */
Cvector HcmtElement::el_field_vec(int k, double x, double z)
{
	Cvector f(12);
	f.init(CC0);
	double a, dxa, dza;
	Complex e, fhx, fhy, fhz, fex, fey, fez;
	Complex fhx0, fhy0, fhz0, fex0, fey0, fez0;
	double r, dx, dz;

	if(vM != 1) { cmterror("HcmtElement::el_field_vec_ac: vM?"); return f; }

	switch(type)
	{
	case MHF: case MHB:
		a = elf(k, z);
		if(fabs(a) < 1.0e-10) return f; 
		// dza = d_elf(k, z);
		dza = d_elf(k, z)+d_alpha*a;
		// e = expi(-d_kz*(z-v_z0));
		e = expi(-d_beta*(z-v_z0));
		fex = mode.cfield(EX, dir, x);
		fey = mode.cfield(EY, dir, x);
		fez = mode.cfield(EZ, dir, x);
		fhx = mode.cfield(HX, dir, x);
		fhy = mode.cfield(HY, dir, x);
		fhz = mode.cfield(HZ, dir, x);
		f(f3Ex) = a*fex*e;
		f(f3Ey) = a*fey*e;
		f(f3Ez) = a*fez*e;
		f(f3Hx) = a*fhx*e;
		f(f3Hy) = a*fhy*e;
		f(f3Hz) = a*fhz*e;
		f(c3Ex) = (CCmI/invommu0*a*fhx-dza*fey)*e;
		f(c3Ey) = (CCmI/invommu0*a*fhy+dza*fex)*e;
		f(c3Ez) =  CCmI/invommu0*a*fhz*e;
		f(c3Hx) = (CCI/invomep0*mode.wg.eps(x)*a*fex-dza*fhy)*e;
		f(c3Hy) = (CCI/invomep0*mode.wg.eps(x)*a*fey+dza*fhx)*e;
		f(c3Hz) =  CCI/invomep0*mode.wg.eps(x)*a*fez*e;
		return f;
		break;

	case MVU: case MVD:
		a = elf(k, x);
		if(fabs(a) < 1.0e-10) return f; 
		// dxa = d_elf(k, x);
		dxa = d_elf(k, x)+d_alpha*a;
		// e = expi(-d_kz*(x-v_z0));
		e = expi(-d_beta*(x-v_z0));
		fex = mode.vertfield(EX, dir, z);
		fey = mode.vertfield(EY, dir, z);
		fez = mode.vertfield(EZ, dir, z);
		fhx = mode.vertfield(HX, dir, z);
		fhy = mode.vertfield(HY, dir, z);
		fhz = mode.vertfield(HZ, dir, z);
		f(f3Ex) = a*fex*e;
		f(f3Ey) = a*fey*e;
		f(f3Ez) = a*fez*e;
		f(f3Hx) = a*fhx*e;
		f(f3Hy) = a*fhy*e;
		f(f3Hz) = a*fhz*e;
		f(c3Ex) = CCmI/invommu0*a*fhx*e;
		f(c3Ey) = (CCmI/invommu0*a*fhy-dxa*fez)*e;
		f(c3Ez) = (CCmI/invommu0*a*fhz+dxa*fey)*e;
		f(c3Hx) = CCI/invomep0*mode.wg.eps(z)*a*fex*e;
		f(c3Hy) = (CCI/invomep0*mode.wg.eps(z)*a*fey-dxa*fhz)*e;
		f(c3Hz) = (CCI/invomep0*mode.wg.eps(z)*a*fez+dxa*fhy)*e;
		return f;
		break;

	case FIM: 
		dx = x-x0;
		dz = z-z0;
		r = sqrt(dx*dx+dz*dz);
		if(r<fimode.cRi) return f;
		if(r>fimode.cRe) return f;
		fimode.emfield(x, z, fex0, fey0, fez0, fhx0, fhy0, fhz0);
		fex =  fex0;
		fey =  fez0;
		fez =  fey0;
		fhx = -fhx0;
		fhy = -fhz0;
		fhz = -fhy0;
		f(f3Ex) = fex;
		f(f3Ey) = fey;
		f(f3Ez) = fez;
		f(f3Hx) = fhx;
		f(f3Hy) = fhy;
		f(f3Hz) = fhz;
		f(c3Ex) = CCmI/invommu0*fhx-CCI*(ky-fimode.beta)*fez;
		f(c3Ey) = CCmI/invommu0*fhy;
		f(c3Ez) = CCmI/invommu0*fhz+CCI*(ky-fimode.beta)*fex;
		f(c3Hx) = CCI/invomep0*fimode.fib.eps(x, z)*fex-CCI*(ky-fimode.beta)*fhz;
		f(c3Hy) = CCI/invomep0*fimode.fib.eps(x, z)*fey;
		f(c3Hz) = CCI/invomep0*fimode.fib.eps(x, z)*fez+CCI*(ky-fimode.beta)*fhx;
		return f;
		break;
	default: 
		cmterror("HcmtElement::el_field_vec: element type not implemented");
		return f;
		break;

	}
	return f;
}

/* prepare the element for the HCMT solver */
void HcmtElement::prepare(Interval cwx, Interval cwz)
{
	given = Ivector(Ns+1);
	given.init(NO);
	u.init(CC0);
	switch(type)
	{	
	case MHF: case MVU:
		u(0) = ain;
		given(0) = YES;
		break;
	case MHB: case MVD:
		u(Ns) = ain;
		given(Ns) = YES;
		break;
	case RES:
		break;
	case CCC: case CCA:
		break;
	case WGM: 
		break;
	case MAT: 
		break;
	case FIM: 
		break;
	default: 
		cmterror("HcmtElement::prepare: type");
		break;
	}
	nseg = Ivector(Ns+1);
	seg = new Isegment*[Ns+1];
	for(int k=0; k<=Ns; ++k)
	{
		seg[k] = new Isegment[MaxNumIsegs];
		nseg(k) = fillintsegs(k, seg[k], cwx, cwz);
	}
	return;
}

/* determine the segments for piecewise evaluation of the FEM integrals */
int HcmtElement::fillintsegs(int k, Isegment *lseg, Interval cwx, Interval cwz)
{
	int lnd = 0; 
	double cwx0 = cwx.x0;
	double cwx1 = cwx.x1;
	double cwz0 = cwz.x0;
	double cwz1 = cwz.x1;

	Isegment tseg;
	double p0, p1;
	double ebeg, ecen, eend;

//	double s, xea, xeb, zea, zeb, x, z;
	double sa, sb, x, z;
	int ano;

	switch(type)
	{
	case MHF: case MHB:
		p0 = mext.x0;
		p1 = mext.x1;
		if(p0<cwx0) p0 = cwx0;
		if(p1>cwx1) p1 = cwx1;
		ebeg = sbeg(k);
		ecen = scen(k);
		eend = send(k);
		if(k == 0)
		{
			tseg.za   = cwz0;
			tseg.zb   = ecen;
			tseg.lbt  = CST; tseg.lbcx = p0;
			tseg.ubt  = CST; tseg.ubcx = p1;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k > 0)
		{
			tseg.za   = ebeg;
			tseg.zb   = ecen;
			tseg.lbt  = CST; tseg.lbcx = p0;
			tseg.ubt  = CST; tseg.ubcx = p1;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k < Ns)
		{
			tseg.za   = ecen;
			tseg.zb   = eend;
			tseg.lbt  = CST; tseg.lbcx = p0;
			tseg.ubt  = CST; tseg.ubcx = p1;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k == Ns)
		{
			tseg.za   = ecen;
			tseg.zb   = cwz1;
			tseg.lbt  = CST; tseg.lbcx = p0;
			tseg.ubt  = CST; tseg.ubcx = p1;
			lseg[lnd] = tseg; ++lnd;
		}
		break;

	case MVU: case MVD:
		p0 = mext.x0;
		p1 = mext.x1;
		if(p0<cwz0) p0 = cwz0;
		if(p1>cwz1) p1 = cwz1;
		ebeg = sbeg(k);
		ecen = scen(k);
		eend = send(k);
		if(k == 0)
		{
			tseg.za   = p0;
			tseg.zb   = p1;
			tseg.lbt  = CST; tseg.lbcx = cwx0;
			tseg.ubt  = CST; tseg.ubcx = ecen;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k > 0)
		{
			tseg.za   = p0;
			tseg.zb   = p1;
			tseg.lbt  = CST; tseg.lbcx = ebeg;
			tseg.ubt  = CST; tseg.ubcx = ecen;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k < Ns)
		{
			tseg.za   = p0;
			tseg.zb   = p1;
			tseg.lbt  = CST; tseg.lbcx = ecen;
			tseg.ubt  = CST; tseg.ubcx = eend;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k == Ns)
		{
			tseg.za   = p0;
			tseg.zb   = p1;
			tseg.lbt  = CST; tseg.lbcx = ecen;
			tseg.ubt  = CST; tseg.ubcx = cwx1;
			lseg[lnd] = tseg; ++lnd;
		}
		break;

	case SRC:
		tseg.za   = z0-scen(k);
		tseg.zb   = z0-sbeg(k);
		tseg.lbt  = CRV;
		tseg.lbcx = x0;
		tseg.lbcz = z0;
		tseg.lbr  = scen(k);
		tseg.lbsg = -1.0;
		tseg.ubt  = CRV;
		tseg.ubcx = x0;
		tseg.ubcz = z0;
		tseg.ubr  = scen(k);
		tseg.ubsg =  1.0;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0-sbeg(k);
		tseg.zb   = z0+sbeg(k);
		tseg.lbt  = CRV;
		tseg.lbcx = x0;
		tseg.lbcz = z0;
		tseg.lbr  = scen(k);
		tseg.lbsg = -1.0;
		tseg.ubt  = CRV;
		tseg.ubcx = x0;
		tseg.ubcz = z0;
		tseg.ubr  = sbeg(k);
		tseg.ubsg = -1.0;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0-sbeg(k);
		tseg.zb   = z0+sbeg(k);
		tseg.lbt  = CRV;
		tseg.lbcx = x0;
		tseg.lbcz = z0;
		tseg.lbr  = sbeg(k);
		tseg.lbsg =  1.0;
		tseg.ubt  = CRV;
		tseg.ubcx = x0;
		tseg.ubcz = z0;
		tseg.ubr  = scen(k);
		tseg.ubsg =  1.0;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0+sbeg(k);
		tseg.zb   = z0+scen(k);
		tseg.lbt  = CRV;
		tseg.lbcx = x0;
		tseg.lbcz = z0;
		tseg.lbr  = scen(k);
		tseg.lbsg = -1.0;
		tseg.ubt  = CRV;
		tseg.ubcx = x0;
		tseg.ubcz = z0;
		tseg.ubr  = scen(k);
		tseg.ubsg =  1.0;
		lseg[lnd] = tseg; ++lnd;

		if(k < Ns)
		{
			tseg.za   = z0-send(k);
			tseg.zb   = z0-scen(k);
			tseg.lbt  = CRV;
			tseg.lbcx = x0;
			tseg.lbcz = z0;
			tseg.lbr  = send(k);
			tseg.lbsg = -1.0;
			tseg.ubt  = CRV;
			tseg.ubcx = x0;
			tseg.ubcz = z0;
			tseg.ubr  = send(k);
			tseg.ubsg =  1.0;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0-scen(k);
			tseg.zb   = z0+scen(k);
			tseg.lbt  = CRV;
			tseg.lbcx = x0;
			tseg.lbcz = z0;
			tseg.lbr  = send(k);
			tseg.lbsg = -1.0;
			tseg.ubt  = CRV;
			tseg.ubcx = x0;
			tseg.ubcz = z0;
			tseg.ubr  = scen(k);
			tseg.ubsg = -1.0;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0-scen(k);
			tseg.zb   = z0+scen(k);
			tseg.lbt  = CRV;
			tseg.lbcx = x0;
			tseg.lbcz = z0;
			tseg.lbr  = scen(k);
			tseg.lbsg =  1.0;
			tseg.ubt  = CRV;
			tseg.ubcx = x0;
			tseg.ubcz = z0;
			tseg.ubr  = send(k);
			tseg.ubsg =  1.0;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0+scen(k);
			tseg.zb   = z0+send(k);
			tseg.lbt  = CRV;
			tseg.lbcx = x0;
			tseg.lbcz = z0;
			tseg.lbr  = send(k);
			tseg.lbsg = -1.0;
			tseg.ubt  = CRV;
			tseg.ubcx = x0;
			tseg.ubcz = z0;
			tseg.ubr  = send(k);
			tseg.ubsg =  1.0;
			lseg[lnd] = tseg; ++lnd;
		}
		if(k == Ns)
		{
			tseg.za   = cwz0;
			tseg.zb   = z0-scen(k);
			tseg.lbt  = CST; tseg.lbcx = cwx0;
			tseg.ubt  = CST; tseg.ubcx = cwx1;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0-scen(k);
			tseg.zb   = z0+scen(k);
			tseg.lbt  = CST; 
			tseg.lbcx = cwx0;
			tseg.ubt  = CRV;
			tseg.ubcx = x0;
			tseg.ubcz = z0;
			tseg.ubr  = scen(k);
			tseg.ubsg = -1.0;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0-scen(k);
			tseg.zb   = z0+scen(k);
			tseg.lbt  = CRV;
			tseg.lbcx = x0;
			tseg.lbcz = z0;
			tseg.lbr  = scen(k);
			tseg.lbsg =  1.0;
			tseg.ubt  = CST; 
			tseg.ubcx = cwx1;
			lseg[lnd] = tseg; ++lnd;

			tseg.za   = z0+scen(k);
			tseg.zb   = cwz1;
			tseg.lbt  = CST; tseg.lbcx = cwx0;
			tseg.ubt  = CST; tseg.ubcx = cwx1;
			lseg[lnd] = tseg; ++lnd;
		}
		break;

	case GBM:
		ebeg = sbeg(k);
		ecen = scen(k);
		eend = send(k);
		tseg.lbt = STL; tseg.ubt = STL;
		tseg.x0 = x0; tseg.z0 = z0; tseg.alpha = alpha;
		if(k == 0) tseg.sa = s0-sqrt((cwx1-cwx0)*(cwx1-cwx0)+(cwz1-cwz0)*(cwz1-cwz0));
		else       tseg.sa = ebeg;
		tseg.sb = ecen;
		tseg.set_za_zb_al(cwx, cwz);
		lseg[lnd] = tseg; ++lnd;
		tseg.sa = ecen;
		if(k == Ns) tseg.sb = sN+sqrt((cwx1-cwx0)*(cwx1-cwx0)+(cwz1-cwz0)*(cwz1-cwz0));
		else        tseg.sb = eend;
		tseg.set_za_zb_al(cwx, cwz);
		lseg[lnd] = tseg; ++lnd;
		break;

	case RES: 
		tseg.za   = z0+mext.x0;
		tseg.zb   = z0-width;
		tseg.lbt  = CST; tseg.lbcx = x0+mext.x0;
		tseg.ubt  = CST; tseg.ubcx = x0+mext.x1;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0-width;
		tseg.zb   = z0+width;
		tseg.lbt  = CST; tseg.lbcx = x0+mext.x0;
		tseg.ubt  = CST; tseg.ubcx = x0-width;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0-width;
		tseg.zb   = z0+width;
		tseg.lbt  = CST; tseg.lbcx = x0-width;
		tseg.ubt  = CST; tseg.ubcx = x0+width;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0-width;
		tseg.zb   = z0+width;
		tseg.lbt  = CST; tseg.lbcx = x0+width;
		tseg.ubt  = CST; tseg.ubcx = x0+mext.x1;
		lseg[lnd] = tseg; ++lnd;

		tseg.za   = z0+width;
		tseg.zb   = z0+mext.x1;
		tseg.lbt  = CST; tseg.lbcx = x0+mext.x0;
		tseg.ubt  = CST; tseg.ubcx = x0+mext.x1;
		lseg[lnd] = tseg; ++lnd;

		// for(int i=0; i<=lnd-1; ++i) (lseg[i]).write(stderr);
		break;

	case CCC: case CCA:

		tseg.cwx0 = cwx.x0;
		tseg.cwx1 = cwx.x1;
		tseg.cwz0 = cwz.x0;
		tseg.cwz1 = cwz.x1;
		tseg.x0 = x0;
		tseg.z0 = z0;
		sa = sbeg(k); 
		sb = scen(k); 
		ano = tseg.set_LPC_slopes(sa, sb);
		lseg[lnd] = tseg; ++lnd;
		if(ano != 0)
		{
			ano = tseg.set_LPC_slopes(sa, sb);
			lseg[lnd] = tseg; ++lnd;
			if(ano != 0) 
			cmterror("Isegment: fillintsegs, CCC: confused (1)");
		}
		sa = scen(k); 
		sb = send(k); 
		ano = tseg.set_LPC_slopes(sa, sb);
		lseg[lnd] = tseg; ++lnd;
		if(ano != 0)
		{
			ano = tseg.set_LPC_slopes(sa, sb);
			lseg[lnd] = tseg; ++lnd;
			if(ano != 0) 
			cmterror("Isegment: fillintsegs, CCC: confused (2)");
		}

// fprintf(stderr, "\n\n El: %d, s in [%g, %g]\n", k, sbeg(k)/PI, send(k)/PI);
// for(int i=0; i<=lnd-1; ++i) (lseg[i]).write(stderr);

/*
		s = sbeg(k); 
		if(fabs(s-0.5*PI) < 1.0e-6)      zea = cwz1;
		else if(fabs(s-1.5*PI) < 1.0e-6) zea = cwz0;
		else if(0.5*PI<s && s<1.5*PI)    zea = z0+(cwx0-x0)*tan(s);
		else                             zea = z0+(cwx1-x0)*tan(s);
		if(fabs(s) < 1.0e-6 || fabs(s-2.0*PI) < 1.0e-6)
		                                 xea = cwx1;
		else if(fabs(s-PI) < 1.0e-6)     xea = cwx0;
		else if(0.0<s && s<PI)           xea = x0+(cwz1-z0)/tan(s);
		else                             xea = x0+(cwz0-z0)/tan(s);
		s = send(k); 
		if(fabs(s-0.5*PI) < 1.0e-6)      zeb = cwz1;
		else if(fabs(s-1.5*PI) < 1.0e-6) zeb = cwz0;
		else if(0.5*PI<s && s<1.5*PI)    zeb = z0+(cwx0-x0)*tan(s);
		else                             zeb = z0+(cwx1-x0)*tan(s);
		if(fabs(s) < 1.0e-6 || fabs(s-2.0*PI) < 1.0e-6)
		                                 xeb = cwx1;
		else if(fabs(s-PI) < 1.0e-6)     xeb = cwx0;
		else if(0.0<s && s<PI)           xeb = x0+(cwz1-z0)/tan(s);
		else                             xeb = x0+(cwz0-z0)/tan(s);

		z = z0; if(zea < z) z = zea; if(zeb < z) z = zeb; 
		if(z<cwz0) z = cwz0;
		tseg.za = z;
		z = z0; if(zea > z) z = zea; if(zeb > z) z = zeb;
		if(z>cwz1) z = cwz1;
		tseg.zb = z;

		x = x0; if(xea < x) x = xea; if(xeb < x) x = xeb; 
		if(x<cwx0) x = cwx0;
		tseg.lbt  = CST; tseg.lbcx = x;
		x = x0; if(xea > x) x = xea; if(xeb > x) x = xeb; 
		if(x>cwx1) x = cwx1;
		tseg.ubt  = CST; tseg.ubcx = x;
		lseg[lnd] = tseg; ++lnd;
// fprintf(stderr, "\n\n El: %d, s in [%g, %g]\n", k, sbeg(k)/PI, send(k)/PI);
// for(int i=0; i<=lnd-1; ++i) (lseg[i]).write(stderr);
*/

		break;

	case WGM: 
		z = z0-WM_R_CROP_OUT;
		tseg.za = z; if(z < cwz0) tseg.za = cwz0;
		z = z0+WM_R_CROP_OUT;
		tseg.zb = z; if(z > cwz1) tseg.zb = cwz1;
		tseg.lbt = CST; 
		x = x0-WM_R_CROP_OUT;
		tseg.lbcx = x; if(x < cwx0) tseg.lbcx = cwx0;
		tseg.ubt = CST; 
		x = x0+WM_R_CROP_OUT;
		tseg.ubcx = x; if(x > cwx1) tseg.ubcx = cwx1;
		lseg[lnd] = tseg; ++lnd;
		break;

	case MAT:
		tseg.za   = cwz0;
		tseg.zb   = cwz1;
		tseg.lbt  = CST; tseg.lbcx = cwx0;
		tseg.ubt  = CST; tseg.ubcx = cwx1;
		lseg[lnd] = tseg; ++lnd;
		break;

	case FIM: 
		z = z0-fimode.cRe;
		tseg.za = z; if(z < cwz0) tseg.za = cwz0;
		z = z0+fimode.cRe;
		tseg.zb = z; if(z > cwz1) tseg.zb = cwz1;
		tseg.lbt = CST; 
		x = x0-fimode.cRe;
		tseg.lbcx = x; if(x < cwx0) tseg.lbcx = cwx0;
		tseg.ubt = CST; 
		x = x0+fimode.cRe;
		tseg.ubcx = x; if(x > cwx1) tseg.ubcx = cwx1;
		lseg[lnd] = tseg; ++lnd;
		break;

	default: 
		cmterror("HcmtElement::fillintsegs: type not implemented");
		break;
	}
	if(lnd >= MaxNumIsegs-1) 
		cmterror("HcmtElement::fillintsegs: MaxNumIsegs exceeded");
	return lnd;
}

/* polish the element, upon completed solver procedure */
void HcmtElement::polish()
{
	switch(type)
	{	
	case MHB: case MVD:
		aout = u(0);
		break;
	case FIM: 
		aout = u(0);
		break;
	default: 
		aout = u(Ns);
		break;
	}
	if(seg != NULL)
	{

		for(int k=0; k<=Ns; ++k) 
		{
			if(seg[k] != NULL) 
			{
				delete[] seg[k];
				seg[k] = NULL;
			}
		}	
		delete[] seg;
		seg = NULL;
		nseg.strip();
	}
	return;
}



/* write element information to stderr */
void HcmtElement::info() const
{
	char c0, c1, c2; hcmtbftchr(type, c0, c1, c2);
	switch(type)
	{	
	case MHF: case MHB: case MVU: case MVD:
		if(mode.vM == 0)
			fprintf(stderr, 
				"%c%c%c T%c%c%c lambda=%g neff=%g |a_in|^2=%g",  
				c0, c1, c2, polCHR(pol), 
				dig10(mode.ord), dig1(mode.ord), 
				2.0*PI/k0, mode.neff, ain.sqabs());
		else 
			fprintf(stderr, 
				"%c%c%c VEC T%c%c%c lambda=%g neff=%g ky=%g kz=%g+i%g |a_in|^2=%g",  
				c0, c1, c2, polCHR(pol), 
				dig10(mode.ord), dig1(mode.ord), 
				2.0*PI/k0, mode.neff, mode.ky.re, d_kz.re, d_kz.im, ain.sqabs());
		break;
	case SRC:
		fprintf(stderr, 
		        "%c%c%c T%c ord=%d x0=%g z0=%g lambda=%g |a_in|^2=%g",
	                c0, c1, c2, polCHR(pol), ord, x0, z0,
		        2.0*PI/k0, ain.sqabs());
		break;
	case GBM:
		fprintf(stderr, 
		"%c%c%c T%c lambda=%g w=%g n=%g al=%g x0=%g z0=%g |a_in|^2=%g",
	          c0, c1, c2, polCHR(pol), 2.0*PI/k0, width, nb, alpha/PI*180.0,
		  x0, z0, ain.sqabs());
		break;
	case RES:
		fprintf(stderr, 
		"%c%c%c T%c lambda=%g w=%g x0=%g z0=%g",
	          c0, c1, c2, polCHR(pol), 2.0*PI/k0, width*2.0, x0, z0);
		break;
	case CCC: case CCA:
		fprintf(stderr, 
		        "%c%c%c T%c lambda=%g beta/k0=%g alpha/k0=%g kappa*R=%g",  
	                c0, c1, c2, polCHR(pol), 
		        2.0*PI/k0, d_beta/k0, d_alpha/k0, d_kappa*bmode.R);
		break;
	case WGM: 
		fprintf(stderr, 
		        "%c%c%c %s lambda_r=%g alpha=%g x0=%g z0=%g",  
	                c0, c1, c2, wgmode.ids, val_lambda(wgmode.omega), wgmode.alpha, x0, z0);
		break;
	case MAT:
		fprintf(stderr, 
		        "%c%c%c T%c lambda=%g x=[%g|%d|%g], z=[%g|%d|%g]",  
	                c0, c1, c2, polCHR(pol), 2.0*PI/k0, x0, numx, xn, z0, numz, zn);
		break;
	case FIM: 
		fprintf(stderr, 
		        "%c%c%c %s neff=%g beta=%g x0=%g z0=%g",  
	                c0, c1, c2, fimode.ids, fimode.neff, fimode.beta, x0, z0);
		break;
	}
	switch(type)
	{	
	case RES: case WGM: case MAT: case FIM: 
		fprintf(stderr, "\n    Amp.\n");
		break;
	default: 
		fprintf(stderr, 
			"\n    FEM-disc.: [%g, %g] #%d\n", s0, sN, Ns+1);
		break;
	}
	return;
}

/* the optical field associated with the element, including FEM coefficients */
Complex HcmtElement::field(Fcomp cp, double x, double z)
{
	Complex f = CC0;
	Complex v;
	if(vM) 
	{
		C3cp c = fcomp2c3cp(cp);
		for(int k=0; k<=Ns; ++k)
		{
			v = u(k);	
			if(v.re != 0.0 || v.im != 0.0) f += v*el_field_vec(c, k, x, z);
		}
		return f;
	}
	FKcp c = fcomp2fkcp(cp);
	for(int k=0; k<=Ns; ++k)
	{
		v = u(k);	
		if(v.re != 0.0 || v.im != 0.0) f += v*el_field(c, k, x, z);
	}
	return f;
}


/* local field "strength",  id: one of mE, mH, qE, qH, mW */
double HcmtField::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 = eps(x, z);
		return  0.25*(  MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			+EP0*eps_xz*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
	default:
		cmterror("HcmtField::lfs, FSid?");
                break;
	}
	return 0.0;
}

/* ------------------------------------------------------------------------
   container for the full hybrid CMT model, including the solver procedure 
   ------------------------------------------------------------------------ */

/* initialize, generic */
HcmtField::HcmtField()
{
	vM = 0;
	ky = 0.0;
	sttype = REC;
	pol = TE;
	lambda = 1.0;
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
	nel = 0;
	evec = NULL;
	cwx = Interval(-1.0, 1.0);
	cwz = Interval(-1.0, 1.0);
	nispwl = 10;
	om.strip();
	smp.strip();
}
/* initialize, sttype = REC */
HcmtField::HcmtField(SegWgStruct s, Polarization p, double wl,
                     Interval Cwx, Interval Cwz, int Nispwl)
{
	vM = 0;
	ky = 0.0;
	sttype = REC;
	str = s;
	pol = p;
	lambda = wl;
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
	nel = 0;
	evec = NULL;
	cwx = Cwx;
	cwz = Cwz;
	nispwl = Nispwl;
	om.strip();
	smp.strip();
}
/* initialize, sttype = OVL */
HcmtField::HcmtField(OvlStruct s, Polarization p, double wl,
                     Interval Cwx, Interval Cwz, int Nispwl)
{
	vM = 0;
	ky = 0.0;
	sttype = OVL;
	ovl = s;
	pol = p;
	lambda = wl;
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
	nel = 0;
	evec = NULL;
	cwx = Cwx;
	cwz = Cwz;
	nispwl = Nispwl;
	Dvector hx, hz;
	hx.strip();
	hz.strip();
	for(int i=0; i<=ovl.no-1; ++i)
	{
		switch(ovl(i).type)
		{	
		case HLAY:
			hx.append(ovl(i).lbd);
			hx.append(ovl(i).ubd);
			break;
		case VLAY:
			hz.append(ovl(i).lbd);
			hz.append(ovl(i).ubd);
			break;
		case RING: 
			hx.append(ovl(i).xo+ovl(i).lbd);
			hx.append(ovl(i).xo+ovl(i).ubd);
			hx.append(ovl(i).xo-ovl(i).lbd);
			hx.append(ovl(i).xo-ovl(i).ubd);
			hz.append(ovl(i).zo+ovl(i).lbd);
			hz.append(ovl(i).zo+ovl(i).ubd);
			hz.append(ovl(i).zo-ovl(i).lbd);
			hz.append(ovl(i).zo-ovl(i).ubd);
			break;
		case DISK:
			hx.append(ovl(i).xo+ovl(i).ubd);
			hx.append(ovl(i).xo-ovl(i).ubd);
			hz.append(ovl(i).zo+ovl(i).ubd);
			hz.append(ovl(i).zo-ovl(i).ubd);
			break;
		case RECT:
			hx.append(ovl(i).lbd);
			hx.append(ovl(i).ubd);
			hz.append(ovl(i).xo);
			hz.append(ovl(i).zo);
			break;
		}
	}
	hx.sort();
	int done;
	done = 0;
	int rl = 0;
	while(hx.nel >= 2 && done != 1)
        {
                done = 1;
                for(int l=0; l<=hx.nel-2; ++l)
                {
                        if(done == 1 && fabs(hx(l)-hx(l+1)) < COMPTOL_HX)
                        {
                                done = 0;
                                rl = l;
                        }
                }
                if(done != 1)
                {
                        for(int l=rl; l<=hx.nel-2; ++l) hx(l) = hx(l+1);
                        --hx.nel;
                }
        }
	if(hx.nel < 2)
	{	
		if(hx.nel < 1) { hx.append(-1.0); hx.append(1.0); }
		else hx.append(hx(0)+1.0);
	}
	hz.sort();
	done = 0;
	rl = 0;
	while(hz.nel >= 2 && done != 1)
        {
                done = 1;
                for(int l=0; l<=hz.nel-2; ++l)
                {
                        if(done == 1 && fabs(hz(l)-hz(l+1)) < COMPTOL_HX)
                        {
                                done = 0;
                                rl = l;
                        }
                }
                if(done != 1)
                {
                        for(int l=rl; l<=hz.nel-2; ++l) hz(l) = hz(l+1);
                        --hz.nel;
                }
        }
	if(hz.nel < 2)
	{	
		if(hz.nel < 1) { hz.append(-1.0); hz.append(1.0); }
		else hz.append(hz(0)+1.0);
	}
	Circuit c(hx.nel-1, hz.nel-1);
	c.hx = hx;	
	c.hz = hz;	
	c.n.init(ovl.bgn);
	for(int j=0; j<=hx.nel; ++j)
	{
		double xa, xb;
		if(j<=0)      xa = -AWAY; else xa = hx(j-1);
		if(j>=hx.nel) xb =  AWAY; else xb = hx(j);
		for(int k=0; k<=hz.nel; ++k)
		{
			double za, zb;
			if(k<=0)      za = -AWAY; else za = hz(k-1);
			if(k>=hz.nel) zb =  AWAY; else zb = hz(k);
			c.n(j, k) = ovl.n(0.5*(xa+xb), 0.5*(za+zb));
		}
	}
	c.lambda = lambda;
	str = SegWgStruct(c);
	om.strip();
	smp.strip();
}
/* initialize, sttype = MIC */
HcmtField::HcmtField(double w, double g, double R, double d, 
                     double x0, double z0, double nb, double ng,
                     Polarization p, double wl,
                     Interval Cwx, Interval Cwz, int Nispwl)
{
	vM = 0;
	ky = 0.0;
	sttype = MIC;
	st_w = w;
	st_g = g;
	st_R = R;
	st_d = d;
	st_x0 = x0;
	st_z0 = z0;
	st_nb = nb;
	st_ng = ng;
	st_eb = nb*nb;
	st_eg = ng*ng;
	
	Waveguide awg(3);
	awg.hx(0) = st_x0-st_R-st_g-st_w;
	awg.hx(1) = st_x0-st_R-st_g;
	awg.hx(2) = st_x0+st_R+st_g;
	awg.hx(3) = st_x0+st_R+st_g+st_w;
	awg.n(0) = st_nb;
	awg.n(1) = st_ng;
	awg.n(2) = st_nb;
	awg.n(3) = st_ng;
	awg.n(4) = st_nb;
	awg.lambda = wl;
	str = SegWgStruct(1);
	str.hz(0) = st_z0-st_R;
	str.hz(1) = st_z0+st_R;
	str(0) = awg;
	str(1) = awg;
	str(2) = awg;

	pol = p;
	lambda = wl;
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
	nel = 0;
	evec = NULL;
	cwx = Cwx;
	cwz = Cwz;
	nispwl = Nispwl;
	om.strip();
	smp.strip();
}

/* map sttype -> local permittivity */
double HcmtField::eps(double x, double z) const
{
	switch(sttype)
	{
		case REC:
			return str.eps(x, z);
			break;
		case OVL:
			return ovl.eps(x, z);
			break;
		case MIC:
			if(x < st_x0-st_R-st_g-st_w) return st_eb;
			if(x > st_x0+st_R+st_g+st_w) return st_eb;
			if(x < st_x0-st_R-st_g) return st_eg;
			if(x > st_x0+st_R+st_g) return st_eg;
			if(x < st_x0-st_R) return st_eb;
			if(x > st_x0+st_R) return st_eb;
			if(z < st_z0-st_R) return st_eb;
			if(z > st_z0+st_R) return st_eb;
			double dx = x-st_x0;
			double dz = z-st_z0;
			double r = sqrt(dx*dx+dz*dz);
			if(r > st_R) return st_eb;
			if(r > st_R-st_d) return st_eg;
			return st_eb;
			break;
	}
	cmterror("HcmtField::eps: sttype");
	return 1.0;
}

/* destroy */
HcmtField::~HcmtField()
{
	if(evec != NULL) delete[] evec;
        evec = NULL;
        nel = 0;
}

/* copy constructor */
HcmtField::HcmtField(const HcmtField& hc)
{
	vM = hc.vM;
	ky = hc.ky;
	sttype = hc.sttype;
	nel = hc.nel;
	str = hc.str;
	ovl = hc.ovl;
	pol = hc.pol;
	lambda = hc.lambda;
	k0 = hc.k0;
	invommu0 = hc.invommu0;
	invomep0 = hc.invomep0;
        evec = new HcmtElement[nel];
        HcmtElement* ap = hc.evec;
        HcmtElement* ep =    evec;
        for(int i=0; i<=nel-1; ++i) *ep++ = *ap++;
	cwx = hc.cwx;
	cwz = hc.cwz;
	nispwl = hc.nispwl;
	st_w = hc.st_w;
	st_g = hc.st_g;
	st_R = hc.st_R;
	st_d = hc.st_d;
	st_x0 = hc.st_x0;
	st_z0 = hc.st_z0;
	st_nb = hc.st_nb;
	st_ng = hc.st_ng;
	st_eb = hc.st_eb;
	st_eg = hc.st_eg;
}

/* assignment */
HcmtField& HcmtField::operator=(const HcmtField& hc)
{
        if(this != &hc)
        {
                if(evec != NULL) delete[] evec;
		vM = hc.vM;
		ky = hc.ky;
		sttype = hc.sttype;
                nel = hc.nel;
		str = hc.str;
		ovl = hc.ovl;
		pol = hc.pol;
		lambda = hc.lambda;
		k0 = hc.k0;
		invommu0 = hc.invommu0;
		invomep0 = hc.invomep0;
                evec = new HcmtElement[nel];
                HcmtElement *ap = hc.evec;
                HcmtElement *ep =    evec;
                for(int i=0; i<=nel-1; ++i) *ep++ = *ap++;
		cwx = hc.cwx;
		cwz = hc.cwz;
		nispwl = hc.nispwl;
		st_w = hc.st_w;
		st_g = hc.st_g;
		st_R = hc.st_R;
		st_d = hc.st_d;
		st_x0 = hc.st_x0;
		st_z0 = hc.st_z0;
		st_nb = hc.st_nb;
		st_ng = hc.st_ng;
		st_eb = hc.st_eb;
		st_eg = hc.st_eg;
        }
        return *this;
}

/* delete all field entries */
void HcmtField::clear()
{
        if(evec != NULL) delete[] evec;
        evec = NULL;
        nel = 0;
}

/* add a constituting field;
   return value: the access index of fe in (*this) */
int HcmtField::addcf(HcmtElement fe)
{
	if(fabs(fe.k0-k0)>1.0e-10 && fe.type != RES 
	                          && fe.type != WGM 
	                          && fe.type != MAT) 
		cmterror("HcmtField::addcf: wavelength mismatch");
	if(fe.vM != vM) cmterror("HcmtField::addcf: vectorzation mismatch");
	if(vM == 0)
	{
		if(fe.pol != pol) cmterror("HcmtField::addcf: polarization mismatch");
	}
	HcmtElement *avec;
        avec = new HcmtElement[nel+1];
        HcmtElement* ap = avec;
        HcmtElement* ep = evec;
        for(int i=0; i<=nel-1; ++i) *ap++ = *ep++;
        *ap = fe;
        if(evec != NULL) delete[] evec;
        evec = avec;
        ++nel;
	(*ap).ky = ky;
        return nel-1;
}

/* add two reciprocal constituting fields, guided modes
   t: MHF for modes that travel forward / backward along z
      MVU for modes that travel upward / downward along x
   m: the mode object
   FEM discretization of the amplitude as a function of
   the x- or z-coordinate in [s_0, s_n], stepsize d_s, regular
   mif, mib (output): element indices of the of the forward/upwards
                      and backward/downwards traveling waves */
void HcmtField::addchannel(HcmtBft t, Mode m, 
                           double s_0, double s_n, double d_s,
                           int& mif, int& mib)
{
	HcmtElement f, b;
	switch(t)
	{	
	case MHF:
		f = HcmtElement(MHF, m, s_0, s_n, d_s);
		b = HcmtElement(MHB, m, s_0, s_n, d_s);
		break;
	case MVU:
		f = HcmtElement(MVU, m, s_0, s_n, d_s);
		b = HcmtElement(MVD, m, s_0, s_n, d_s);
		break;
	default: 
		cmterror("HcmtField, addchannel, type");
		break;
	}
	mif = addcf(f); 
	mib = addcf(b);
	(evec[mif]).recEli = mib;
	(evec[mib]).recEli = mif;
	return;
}

/* Variational approach: integration over products of elements */

/* integral kernel */
Cvector HcmtField::fkernel(int l, int m, double x, double z) const
{
	if(vM != 0) return fkernel_vec(l, m, x, z); 
	Cvector s(2); s.init(CC0);
	Cvector lf, mf;
	Complex ley, lhx, lhz, mey, mhx, mhz;
	Complex lhy, lex, lez, mhy, mex, mez;
	double e = eps(x, z);
	Complex inndep0(0.0, e/invomep0);
	Complex idmu0(0.0, 1.0/invommu0);

	if(cckernel < 0) goto cc_supermode_kernel;
	if(cckernel == 1) goto cc_kernel;

	if(pol != TE) goto fkernel_TM;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	ley = lf(fcEy);
	lhx = lf(fcHx);
	lhz = lf(fcHz);
	mey = mf(fcEy);
	mhx = mf(fcHx);
	mhz = mf(fcHz);
	s(0) += ley * mf(dzHxmdxHz);
	s(0) -= lhx * mf(dzEy);
	s(0) += lhz * mf(dxEy);
	s(0) -= ley * mey*inndep0;
	s(0) += lhx * mhx*idmu0;
	s(0) += lhz * mhz*idmu0;

	s(1) += mey * lf(dzHxmdxHz);
	s(1) -= mhx * lf(dzEy);
	s(1) += mhz * lf(dxEy);
	s(1) -= mey * ley*inndep0;
	s(1) += mhx * lhx*idmu0;
	s(1) += mhz * lhz*idmu0;

	return s;	

	fkernel_TM: ;;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	lhy = lf(fcHy);
	lex = lf(fcEx);
	lez = lf(fcEz);
	mhy = mf(fcHy);
	mex = mf(fcEx);
	mez = mf(fcEz);
	s(0) += lhy * mf(dzExmdxEz);
	s(0) -= lex * mf(dzHy);
	s(0) += lez * mf(dxHy);
	s(0) += lhy * mhy*idmu0;
	s(0) -= lex * mex*inndep0;
	s(0) -= lez * mez*inndep0;

	s(1) += mhy * lf(dzExmdxEz);
	s(1) -= mex * lf(dzHy);
	s(1) += mez * lf(dxHy);
	s(1) += mhy * lhy*idmu0;
	s(1) -= mex * lex*inndep0;
	s(1) -= mez * lez*inndep0;

	return s;	

	cc_kernel: ;;

	if(pol != TE) goto cc_fkernel_TM;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	ley = lf(fcEy);
	lhx = lf(fcHx);
	lhz = lf(fcHz);
	mey = mf(fcEy);
	mhx = mf(fcHx);
	mhz = mf(fcHz);
	s(0) += ley.conj() * mf(dzHxmdxHz);
	s(0) += lhx.conj() * mf(dzEy);
	s(0) -= lhz.conj() * mf(dxEy);
	s(0) -= ley.conj() * mey*inndep0;
	s(0) -= lhx.conj() * mhx*idmu0;
	s(0) -= lhz.conj() * mhz*idmu0;
	s(1) += mey.conj() * lf(dzHxmdxHz);
	s(1) += mhx.conj() * lf(dzEy);
	s(1) -= mhz.conj() * lf(dxEy);
	s(1) -= mey.conj() * ley*inndep0;
	s(1) -= mhx.conj() * lhx*idmu0;
	s(1) -= mhz.conj() * lhz*idmu0;

	return s;	

	cc_fkernel_TM: ;;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	lhy = lf(fcHy);
	lex = lf(fcEx);
	lez = lf(fcEz);
	mhy = mf(fcHy);
	mex = mf(fcEx);
	mez = mf(fcEz);
	s(0) -= lhy.conj() * mf(dzExmdxEz);
	s(0) -= lex.conj() * mf(dzHy);
	s(0) += lez.conj() * mf(dxHy);
	s(0) -= lhy.conj() * mhy*idmu0;
	s(0) -= lex.conj() * mex*inndep0;
	s(0) -= lez.conj() * mez*inndep0;

	s(1) -= mhy.conj() * lf(dzExmdxEz);
	s(1) -= mex.conj() * lf(dzHy);
	s(1) += mez.conj() * lf(dxHy);
	s(1) -= mhy.conj() * lhy*idmu0;
	s(1) -= mex.conj() * lex*inndep0;
	s(1) -= mez.conj() * lez*inndep0;

	return s;	

	cc_supermode_kernel: ;;

	if(pol != TE) goto cc_supermode_kernel_TM;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	ley = lf(fcEy);
	lhx = lf(fcHx);
	lhz = lf(fcHz);
	mey = mf(fcEy);
	mhx = mf(fcHx);
	mhz = mf(fcHz);
	if(cckernel == -1)
	{
		s(0) += ley.conj() * mf(dzHxmdxHz);
		s(0) += lhx.conj() * mf(dzEy);
		s(0) -= lhz.conj() * mf(dxEy);
		s(1) += mey.conj() * lf(dzHxmdxHz);
		s(1) += mhx.conj() * lf(dzEy);
		s(1) -= mhz.conj() * lf(dxEy);
		return s;	
	}
	if(cckernel == -2)
	{
		inndep0 = Complex(0.0, e*EP0);
		idmu0 = Complex(0.0, MU0);
		s(0) += ley.conj() * mey*inndep0;
		s(0) += lhx.conj() * mhx*idmu0;
		s(0) += lhz.conj() * mhz*idmu0;
		s(1) += mey.conj() * ley*inndep0;
		s(1) += mhx.conj() * lhx*idmu0;
		s(1) += mhz.conj() * lhz*idmu0;
		return s;	
	}
	return s;	

	cc_supermode_kernel_TM: ;;

	lf = (evec[U_bfi(l)]).el_field(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field(U_eli(m), x, z);
	lhy = lf(fcHy);
	lex = lf(fcEx);
	lez = lf(fcEz);
	mhy = mf(fcHy);
	mex = mf(fcEx);
	mez = mf(fcEz);
	if(cckernel == -1)
	{
		s(0) -= lhy.conj() * mf(dzExmdxEz);
		s(0) -= lex.conj() * mf(dzHy);
		s(0) += lez.conj() * mf(dxHy);
		s(1) -= mhy.conj() * lf(dzExmdxEz);
		s(1) -= mex.conj() * lf(dzHy);
		s(1) += mez.conj() * lf(dxHy);
		return s;	
	}
	if(cckernel == -2)
	{
		inndep0 = Complex(0.0, e*EP0);
		idmu0 = Complex(0.0, MU0);
		s(0) += lhy.conj() * mhy*idmu0;
		s(0) += lex.conj() * mex*inndep0;
		s(0) += lez.conj() * mez*inndep0;
		s(1) += mhy.conj() * lhy*idmu0;
		s(1) += mex.conj() * lex*inndep0;
		s(1) += mez.conj() * lez*inndep0;
		return s;	
	}
	return s;	
}

/* integral kernel, vectorial variant */
Cvector HcmtField::fkernel_vec(int l, int m, double x, double z) const
{
	Cvector s(2); s.init(CC0);
	if(vM != 1) { cmterror("HcmtField::fkernel_vec: vM?"); return s; }
	if(cckernel == 0) { cmterror("HcmtField::fkernel_vec: cckernel=0"); return s; }
	Cvector lf, mf;
	Complex lex, ley, lez, lhx, lhy, lhz;
	Complex mex, mey, mez, mhx, mhy, mhz;
	double e = eps(x, z);
	Complex idep0(0.0, e/invomep0);
	Complex idmu0(0.0, 1.0/invommu0);

	lf = (evec[U_bfi(l)]).el_field_vec(U_eli(l), x, z);
	if(m == l) mf = lf;
	else       mf = (evec[U_bfi(m)]).el_field_vec(U_eli(m), x, z);
	lex = lf(f3Ex);
	ley = lf(f3Ey);
	lez = lf(f3Ez);
	lhx = lf(f3Hx);
	lhy = lf(f3Hy);
	lhz = lf(f3Hz);
	mex = mf(f3Ex);
	mey = mf(f3Ey);
	mez = mf(f3Ez);
	mhx = mf(f3Hx);
	mhy = mf(f3Hy);
	mhz = mf(f3Hz);

	if(cckernel < 0) goto cc_supermode_kernel;

	s(0) = (lex.conj()*mf(c3Hx)+ley.conj()*mf(c3Hy)+lez.conj()*mf(c3Hz))
	     - (lhx.conj()*mf(c3Ex)+lhy.conj()*mf(c3Ey)+lhz.conj()*mf(c3Ez))
	     - (lex.conj()*mex     +ley.conj()*mey     +lez.conj()*mez)*idep0
	     - (lhx.conj()*mhx     +lhy.conj()*mhy     +lhz.conj()*mhz)*idmu0;
	s(1) = (mex.conj()*lf(c3Hx)+mey.conj()*lf(c3Hy)+mez.conj()*lf(c3Hz))
	     - (mhx.conj()*lf(c3Ex)+mhy.conj()*lf(c3Ey)+mhz.conj()*lf(c3Ez))
	     - (mex.conj()*lex     +mey.conj()*ley     +mez.conj()*lez)*idep0
	     - (mhx.conj()*lhx     +mhy.conj()*lhy     +mhz.conj()*lhz)*idmu0;
	return s;	

	cc_supermode_kernel: ;;

	if(cckernel == -1)
	{
		s(0) = (lex.conj()*mf(c3Hx)+ley.conj()*mf(c3Hy)+lez.conj()*mf(c3Hz))
		     - (lhx.conj()*mf(c3Ex)+lhy.conj()*mf(c3Ey)+lhz.conj()*mf(c3Ez));
		s(1) = (mex.conj()*lf(c3Hx)+mey.conj()*lf(c3Hy)+mez.conj()*lf(c3Hz))
		     - (mhx.conj()*lf(c3Ex)+mhy.conj()*lf(c3Ey)+mhz.conj()*lf(c3Ez));
		return s;	
	}
	if(cckernel == -2)
	{
		idep0 = Complex(0.0, e*EP0);
		idmu0 = Complex(0.0, MU0);
		s(0) = (lex.conj()*mex     +ley.conj()*mey     +lez.conj()*mez)*idep0
		     + (lhx.conj()*mhx     +lhy.conj()*mhy     +lhz.conj()*mhz)*idmu0;
		s(1) = (mex.conj()*lex     +mey.conj()*ley     +mez.conj()*lez)*idep0
		     + (mhx.conj()*lhx     +mhy.conj()*lhy     +mhz.conj()*lhz)*idmu0;
		return s;	
	}
	return s;	
}

/* integrate over x, one continuous piece */
Cvector HcmtField::fpintegx(int l, int m, double xa, double xb, double z) const
{
	double xr, xm, dx, x;
	static double xc[] = {0.0, 0.1488743389, 0.4333953941,
		              0.6794095682, 0.8650633666, 0.9739065285};
	static double wc[] = {0.0, 0.2955242247, 0.2692667193,
		              0.2190863625, 0.1494513491, 0.0666713443};
	Cvector sum(2); sum.init(CC0); 
	double a, b, step;
	int numx; 
	numx = (int) ceil((xb-xa)/lambda*nispwl);
	step = (xb-xa)/((double) numx);
	Cvector s, fval; 
	for(int i=0; i<=numx-1; ++i)
	{
		a  = xa+((double) i)*step;
		b  = a+step;
		xm = 0.5*(b+a);
		xr = 0.5*(b-a);
		s = Cvector(2); s.init(CC0);
		for(int j=1; j<=5; j++) 
		{
			dx = xr*xc[j];
			   x  = xm+dx;
			fval  = fkernel(l, m, x, z);
			   x  = xm-dx;
			fval.addeq(fkernel(l, m, x, z));
			s.addeq(mult(fval, wc[j]));
		}
		sum.addeq(mult(s, xr));
	}
	return sum;
}

/* integrate over x, piecewise over parts of the structure */
Cvector HcmtField::fintegx(int l, int m, double xa, double xb, double z) const
{
	Cvector sum(2); sum.init(CC0); 
	if(xb-xa < HDIST) return sum;
	Waveguide wg = str(str.segidx(z));
	double x0, x1;
	x0 = xa; x1 = xa;
	int li = wg.layeridx(x0+HDIST);
	Interval lay = wg.layer(li);
	while(x1 < xb-HDIST)
	{
		x1 = lay.x1;
		if(x1 > xb) x1 = xb;
		sum.addeq(fpintegx(l, m, x0, x1, z));
		x0 = x1;
		++li;
		lay = wg.layer(li);
	}
	return sum;
}

/* integrate over z, one continuous piece */
Cvector HcmtField::fpinteg(int l, int m, 
                           const Isegment& lseg, const Isegment& mseg, 
                           double za, double zb) const
{
	Cvector sum(2); sum.init(CC0); 
	double xa, xb, txa, txb;
	double zr, zm, dz, z;
	static double zc[] = {0.0, 0.1488743389, 0.4333953941,
		              0.6794095682, 0.8650633666, 0.9739065285};
	static double wc[] = {0.0, 0.2955242247, 0.2692667193,
		              0.2190863625, 0.1494513491, 0.0666713443};
	double a, b, step;
	int numz; 
	numz = (int) ceil((zb-za)/lambda*nispwl);
	step = (zb-za)/((double) numz);
	Cvector s, fval; 
	for(int i=0; i<=numz-1; ++i)
	{
		a  = za+((double) i)*step;
		b  = a+step;
		zm = 0.5*(b+a);
		zr = 0.5*(b-a);
		s = Cvector(2); s.init(CC0);
		for(int j=1; j<=5; j++) 
		{
			dz = zr*zc[j];

			z  = zm+dz;
			 xa = lseg.lbxpos(z);
			 xb = lseg.ubxpos(z);
			txa = mseg.lbxpos(z);
			txb = mseg.ubxpos(z);
			if(txa > xa) xa = txa;
			if(txb < xb) xb = txb;
		 	fval = fintegx(l, m, xa, xb, z);

			z  = zm-dz;
			 xa = lseg.lbxpos(z);
			 xb = lseg.ubxpos(z);
			txa = mseg.lbxpos(z);
			txb = mseg.ubxpos(z);
			if(txa > xa) xa = txa;
			if(txb < xb) xb = txb;
			fval.addeq(fintegx(l, m, xa, xb, z));
			s.addeq(mult(fval,wc[j]));
		}
		sum.addeq(mult(s,zr));
	}
	return sum;
}

/* HCMT system, functional kernel, 
   integrate over z, piecewise over parts of the structure */
Cvector HcmtField::finteg(int l, int m, 
                          const Isegment& lseg, const Isegment& mseg) const
{
	Cvector sum(2); sum.init(CC0); 
	double za = lseg.za;
	double zb = lseg.zb;
	if(mseg.za > za) za = mseg.za;
	if(mseg.zb < zb) zb = mseg.zb;
	if(zb-za < HDIST) return sum;
	double z0, z1;
	z0 = za; z1 = za;
	int si = str.segidx(z0+HDIST);
	Interval sg = str.segment(si);
	while(z1 < zb-HDIST)
	{
		z1 = sg.x1;
		if(z1 > zb) z1 = zb;
		sum.addeq(fpinteg(l, m, lseg, mseg, z0, z1));
		z0 = z1;
		++si;
		sg = str.segment(si);
	}
	return sum;
}

/* HCMT system, quadratic boundary terms */
Cvector HcmtField::btbeval(int l, int m) const
{
	Cvector b(2); b.init(CC0); 
	int lbfi = U_bfi(l);
	int mbfi = U_bfi(m);
	if(mbfi != (evec[lbfi]).recEli) return b;
	int leli = U_eli(l);
	int meli = U_eli(m);
	if(leli != meli) return b; 
	int lns = (evec[lbfi]).Ns;
	if(leli != 0 && leli != lns) return b; 
	Complex ov;
	switch((evec[lbfi]).type)
	{
		case MHF: case MHB: case MVU: case MVD:
			ov = -2.0*overlap_ncns((evec[lbfi]).mode, FORW,
			                      (evec[lbfi]).mode, FORW);
			b(0) = ov;
			b(1) = ov;
			return b;
			break;
		default: 
			cmterror("HcmtField::btbeval: type");
			break;
	}
	return b;
}

/* HCMT system, linear boundary terms, inhomogenities */
Complex HcmtField::btreval(int l) const
{
	int bfi = U_bfi(l);
	int ns = (evec[bfi]).Ns;
	int eli = U_eli(l);
	if(eli != 0 && eli != ns) return CC0; 
	Complex r = CC0; 
	switch((evec[bfi]).type)
	{
		case MHF: 
			if(eli == ns) 
			{
				int re = (evec[bfi]).recEli;
				Complex a = (evec[re]).ain;
				Complex ov = 
					overlap_ncns((evec[bfi]).mode, FORW,
					             (evec[bfi]).mode, FORW);
				return -8.0*a*ov;
			}
			return CC0;
			break;
		case MHB:
			if(eli == 0) 
			{
				int re = (evec[bfi]).recEli;
				Complex a = (evec[re]).ain;
				Complex ov = 
				        overlap_ncns((evec[bfi]).mode, FORW,
				                     (evec[bfi]).mode, FORW);
				return -8.0*a*ov;
			}
			return CC0;
			break;
		case MVU: 
			if(eli == ns) 
			{
				int re = (evec[bfi]).recEli;
				Complex a = (evec[re]).ain;
				Complex ov = 
				        overlap_ncns((evec[bfi]).mode, FORW,
				                     (evec[bfi]).mode, FORW);
				return -8.0*a*ov;
			}
			return CC0;
			break;
		case MVD:
			if(eli == 0) 
			{
				int re = (evec[bfi]).recEli;
				Complex a = (evec[re]).ain;
				Complex ov = 
				        overlap_ncns((evec[bfi]).mode, FORW,
				                     (evec[bfi]).mode, FORW);
				return -8.0*a*ov;
			}
			return CC0;
			break;
		default: 
			cmterror("HcmtField::btreval: type");
			break;
	}
	return r;
}

/* the HCMT solver, complex conjugate fields, Galerkin procedure */
void HcmtField::init()
{
	fprintf(stderr, "\n");
	fprintf(stderr, "------------- HCMT ---------------------------------------------'06, '20 ---\n");
	if(vM == 0)
	{
		fprintf(stderr, 
			"T%c  x in (%g, %g)  z in (%g, %g)  nispwl: %d\n",
			polCHR(pol), cwx.x0, cwx.x1, cwz.x0, cwz.x1, nispwl);
	}
	else 
	{
	fprintf(stderr, 
	        "VEC  ky = %g  x in (%g, %g)  z in (%g, %g)  nispwl: %d\n",
	        ky, cwx.x0, cwx.x1, cwz.x0, cwz.x1, nispwl);
	}
	Circuit rc = str.circuit();
switch(sttype)
{
	case REC:
	fprintf(stderr, "------------- circuit ------------------------------------------------------\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");
}
	break;
	case OVL:
	fprintf(stderr, "------------- overlay structure --------------------------------------------\n");
	ovl.write(stderr);
	break;
	case MIC:
	fprintf(stderr, "------------- circular microcavity -----------------------------------------\n");
fprintf(stderr, "w: %6.10g  ", st_w);
fprintf(stderr, "g: %6.10g  ", st_g);
fprintf(stderr, "R: %6.10g  ", st_R);
fprintf(stderr, "d: %6.10g  ", st_d);
fprintf(stderr, "\n");
fprintf(stderr, "x0: %6.10g  ", st_x0);
fprintf(stderr, "z0: %6.10g  ", st_z0);
fprintf(stderr, "nb: %6.10g  ", st_nb);
fprintf(stderr, "ng: %6.10g  ", st_ng);
fprintf(stderr, "\n");
	break;
	default: 
	cmterror("HcmtField::solve: invalid sttype");
	break;
}
fprintf(stderr, "lambda: %.10g  k0: %g", rc.lambda, val_k0(rc.lambda));
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
	fprintf(stderr, "------------- elements -----------------------------------------------------\n");
	for(int j=0; j<=nel-1; ++j) (evec[j]).info();
	fprintf(stderr, "----------------------------------------------------------------------------\n");
	
	NumU=0; for(int j=0; j<=nel-1; ++j) NumU += (evec[j]).u.nel;
	U_bfi = Ivector(NumU);
	U_eli = Ivector(NumU);
	U_val = Cvector(NumU);

	u = 0; 
	g = NumU-1;	
	for(int j=0; j<=nel-1; ++j)
	{
		(evec[j]).prepare(cwx, cwz);
		for(int k=0; k<=(evec[j]).Ns; ++k)
		{
			if((evec[j]).given(k)==NO)
			{
				U_bfi(u) = j;
				U_eli(u) = k;
				U_val(u) = (evec[j]).u(k);
				++u;
			}
			else
			{
				U_bfi(g) = j;
				U_eli(g) = k;
				U_val(g) = (evec[j]).u(k);
				--g;
			}
		}
		
	}
	if(u!=g+1) cmterror("HcmtField::solve: unknowns, numbering");

	g = NumU-u;
	// unknowns: u, indices 0 -> u-1 
	// given:    g, indices u -> NumU-1
	return;
}

/* the HCMT solver, complex conjugate fields, Galerkin procedure */
void HcmtField::setupsystem(Cmatrix& M)
{
	fprintf(stderr, "HybCMTsolver_Gal(%d, %du, %dg):\n", NumU, u, g);

	M = Cmatrix(NumU, NumU);
	M.init(CC0);

	cckernel = 1;

	int lnd, mnd, lbfi, mbfi, leli, meli;
	for(int l=0; l<=NumU-1; ++l)
	{
		lbfi = U_bfi(l);
		leli = U_eli(l);
		lnd = (evec[lbfi]).nseg(leli);
		for(int m=l; m<=NumU-1; ++m)
		{
			mbfi = U_bfi(m);
			meli = U_eli(m);
			mnd = (evec[mbfi]).nseg(meli);
			Cvector mc(2); mc.init(CC0);
			for(int lr=0; lr<=lnd-1; ++lr)
			{
				for(int mr=0; mr<=mnd-1; ++mr)
				{
			   		mc.addeq(finteg(l, m, 
				                (evec[lbfi]).seg[leli][lr],
				                (evec[mbfi]).seg[meli][mr]));
				}
			}
			M(l, m) = mc(0);
			M(m, l) = mc(1);
		}
		fprintf(stderr, ".");
	}
	fprintf(stderr, " .\n");
	return;
}

void HcmtField::solve(Cmatrix M)
{
	fprintf(stderr, "*");
	Cvector uv;
	Cvector gv(g);
	for(int i=0; i<=g-1; ++i) gv(i) = U_val(u+i);

	Cmatrix Fu, Fg, Fut, FutFu;
	Fu = M.submatrix(u+g, u, 0, 0);
	Fg = M.submatrix(u+g, g, 0, u);

	Fut = Fu;
	Fut.adjoint();
	FutFu = mult(Fut, Fu);
	Cvector rhs;
	rhs = mult(Fut, mult(Fg, gv));
	fprintf(stderr, "|");
	uv = FutFu.CDsolve(rhs);
	uv.multeq(-1.0);
	fprintf(stderr, ".");

	for(int i=0; i<=u-1; ++i) U_val(i) = uv(i);
	for(int r=0; r<=NumU-1; ++r) (evec[U_bfi(r)]).u(U_eli(r)) = U_val(r);
	for(int j=0; j<=nel-1; ++j) (evec[j]).polish();
	fprintf(stderr, " Ok.\n");
	return;
}


void HcmtField::solve()
{
	Cmatrix M; 
	init();
	setupsystem(M);
	solve(M);
	return;
}


/* the HCMT solver, non complex conjugated fields,
                    variational scheme incl. boundary terms;
                !!! bidirectional versions of all fields required,
                <-> addchannel, field references */
void HcmtField::varsolve()
{
	fprintf(stderr, "\n");
	fprintf(stderr, "------------- variational HCMT ------------------------------------- '07 ---\n");
	fprintf(stderr, 
	        "T%c  x in (%g, %g)  z in (%g, %g)  nispwl: %d\n",
	        polCHR(pol), cwx.x0, cwx.x1, cwz.x0, cwz.x1, nispwl);
	Circuit rc = str.circuit();
switch(sttype)
{
	case REC:
	fprintf(stderr, "------------- circuit ------------------------------------------------------\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");
}
	break;
	case MIC:
	fprintf(stderr, "------------- circular microcavity -----------------------------------------\n");
fprintf(stderr, "w: %6.10g  ", st_w);
fprintf(stderr, "g: %6.10g  ", st_g);
fprintf(stderr, "R: %6.10g  ", st_R);
fprintf(stderr, "d: %6.10g  ", st_d);
fprintf(stderr, "\n");
fprintf(stderr, "x0: %6.10g  ", st_x0);
fprintf(stderr, "z0: %6.10g  ", st_z0);
fprintf(stderr, "nb: %6.10g  ", st_nb);
fprintf(stderr, "ng: %6.10g  ", st_ng);
fprintf(stderr, "\n");
	cmterror("HcmtField::varsolve: MIC not implemented");
	break;
	default: 
	cmterror("HcmtField::varsolve: invalid sttype");
	break;
}
fprintf(stderr, "lambda: %.10g  k0: %g", rc.lambda, val_k0(rc.lambda));
if(rc.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
	fprintf(stderr, "------------- elements -----------------------------------------------------\n");
	for(int j=0; j<=nel-1; ++j) (evec[j]).info();
	fprintf(stderr, "----------------------------------------------------------------------------\n");

	if(vM != 0) cmterror("HcmtField::varsolve: vectorial variant not implemented");
	
	NumU=0; for(int j=0; j<=nel-1; ++j) NumU += (evec[j]).u.nel;
	U_bfi = Ivector(NumU);
	U_eli = Ivector(NumU);
	U_val = Cvector(NumU);

	int u = 0; 
	int g = NumU-1;	
	for(int j=0; j<=nel-1; ++j)
	{
		(evec[j]).prepare(cwx, cwz);
		for(int k=0; k<=(evec[j]).Ns; ++k)
		{
			if((evec[j]).given(k)==NO)
			{
				U_bfi(u) = j;
				U_eli(u) = k;
				U_val(u) = (evec[j]).u(k);
				++u;
			}
			else
			{
				U_bfi(g) = j;
				U_eli(g) = k;
				U_val(g) = (evec[j]).u(k);
				--g;
			}
		}
		
	}
	if(u!=g+1) cmterror("HcmtField::solve: unknowns, numbering");

	g = NumU-u;
	// unknowns: u, indices 0 -> u-1 
	// given:    g, indices u -> NumU-1

	fprintf(stderr, "HybCMTsolver_var(%d, %du, %dg):\n", NumU, u, g);

	Cmatrix F(NumU, NumU);
	F.init(CC0);
	Cvector R(NumU);
	R.init(CC0);

	cckernel = 0;

	int lnd, mnd, lbfi, mbfi, leli, meli;
	for(int l=0; l<=NumU-1; ++l)
	{
		lbfi = U_bfi(l);
		leli = U_eli(l);
		lnd = (evec[lbfi]).nseg(leli);
		for(int m=l; m<=NumU-1; ++m)
		{
			mbfi = U_bfi(m);
			meli = U_eli(m);
			mnd = (evec[mbfi]).nseg(meli);
			Cvector mc(2); mc.init(CC0);
			for(int lr=0; lr<=lnd-1; ++lr)
			{
				for(int mr=0; mr<=mnd-1; ++mr)
				{
			   		mc.addeq(finteg(l, m, 
				                (evec[lbfi]).seg[leli][lr],
				                (evec[mbfi]).seg[meli][mr]));
				}
			}
			F(l, m) = mc(0);
			F(m, l) = mc(1);
			
			mc = btbeval(l, m);
			F(l, m) += mc(0);
			F(m, l) += mc(1);

		}
		fprintf(stderr, ".");
		R(l) = btreval(l);
	}
	fprintf(stderr, "|");

	Cvector uv(u);
        Cvector gv(g);
        for(int i=0; i<=g-1; ++i) gv(i) = U_val(u+i);

	Cmatrix Ft, FpFt;
	Ft = F; Ft.transpose();
	FpFt = add(F, Ft);

	Cmatrix Mu, Mg, Mut, MutMu;
	Mu = FpFt.submatrix(u+g, u, 0, 0);
	Mg = FpFt.submatrix(u+g, g, 0, u);

	Mut = Mu;
	Mut.adjoint();
	MutMu = mult(Mut, Mu);
	fprintf(stderr, ".");

	Cvector rhs;
	rhs = mult(Mut, sub(R, mult(Mg, gv)));
	uv = MutMu.CDsolve(rhs);

	fprintf(stderr, ".");

	for(int i=0; i<=u-1; ++i) U_val(i) = uv(i);
	for(int r=0; r<=NumU-1; ++r) (evec[U_bfi(r)]).u(U_eli(r)) = U_val(r);
	for(int j=0; j<=nel-1; ++j) (evec[j]).polish();
	fprintf(stderr, " Ok.\n");
}


/* HCMT supermode analysis, complex conjugate fields, Galerkin procedure;
   returns a vector of supermode eigenfrequencies */
Cvector HcmtField::supermodes()
{
	Cmatrix A; 
	Cmatrix B; 
	init();
	fprintf(stderr, "HCMT_supermodes(%d, %du, %dg):\n", NumU, u, g);
	for(int i=u; i<=NumU-1; ++i) U_val(i) = 0;

	A = Cmatrix(u, u);
	A.init(CC0);
	B = Cmatrix(u, u);
	B.init(CC0);

	int lnd, mnd, lbfi, mbfi, leli, meli;
	for(int l=0; l<=u-1; ++l)
	{
		lbfi = U_bfi(l);
		leli = U_eli(l);
		lnd = (evec[lbfi]).nseg(leli);
		for(int m=l; m<=u-1; ++m)
		{
			mbfi = U_bfi(m);
			meli = U_eli(m);
			mnd = (evec[mbfi]).nseg(meli);
			Cvector ac(2); ac.init(CC0);
			Cvector bc(2); bc.init(CC0);
			for(int lr=0; lr<=lnd-1; ++lr)
			{
				for(int mr=0; mr<=mnd-1; ++mr)
				{
					cckernel = -1;
			   		ac.addeq(finteg(l, m, 
				                (evec[lbfi]).seg[leli][lr],
				                (evec[mbfi]).seg[meli][mr]));
					cckernel = -2;
			   		bc.addeq(finteg(l, m, 
				                (evec[lbfi]).seg[leli][lr],
				                (evec[mbfi]).seg[meli][mr]));
				}
			}
			A(l, m) = ac(0);
			A(m, l) = ac(1);
			B(l, m) = bc(0);
			B(m, l) = bc(1);
		}
		fprintf(stderr, ".");
	}
	fprintf(stderr, "*");

	B.inverse();
	smp = mult(B, A);
	int info;
	om = smp.geneigenv(info);
	
	if(info == 0) fprintf(stderr, ".\n");
	else cmterror("geneigenv: Failed");
	
	for(int m=0; m<=u-1; ++m)
	{

		fprintf(stderr, "[%d] lam=%g, om=%g, al=%g, Q=%g, D_lam=%g\n",
                m, val_lambda(om(m).re), om(m).re, om(m).im, theQ(om(m)), FWHM_lambda(om(m)));
		if(u <= 10)
		{
			for(int j=0; j<=u-1; ++j) 
			{
				Complex a = (smp.col(m))(j);
				HcmtBft t = (evec[U_bfi(j)]).type;
				if(t == WGM) 
					fprintf(stderr, "   %s: |a| = %g  a = %g+i%g\n", (evec[U_bfi(j)]).wgmode.ids, a.abs(), a.re, a.im);
			}
		}
	}
	fprintf(stderr, "Ok.\n");
	return om;
}

/* filter the set of supermodes: returns numbers of smodes 
   with resonant wavelengths between minwl and maxwl, and
   with attenuation below limal */
Ivector HcmtField::sm_filter(double minwl, double maxwl, double limal)
{
	Ivector rsmi;
	rsmi.strip();
	for(int j=0; j<=om.nel-1; ++j)
	{
		double l = val_lambda(om(j).re);
		if(l >= minwl && l <= maxwl && om(j).im < limal)
			rsmi.append(j);
	}
	return rsmi;
}

/* setup internal representation of supermode no. m, 
   return value: the respective: eigenfrequency */
Complex HcmtField::supermode(int m)
{
	if(m < 0 || m > om.nel-1) cmterror("supermode: m");
	if(m >= smp.c || m >= smp.r) cmterror("supermode: m");
	lambda = val_lambda(om(m).re);
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
	for(int j=0; j<=str.nz+1; ++j) str(j).lambda = lambda;
	ovl.lambda = lambda;
	for(int i=0; i<=u-1; ++i) U_val(i) = (smp.col(m))(i);
	for(int r=0; r<=NumU-1; ++r) (evec[U_bfi(r)]).u(U_eli(r)) = U_val(r);
	for(int j=0; j<=nel-1; ++j) (evec[j]).polish();
	return om(m);
}

/* the resulting field */
Complex HcmtField::field(Fcomp cp, double x, double z) const
{
	Complex f;
	for(int m=0; m<=nel-1; ++m) f += (evec[m]).field(cp, x, z);
	return f;
}

/* values on a rectangular mesh */
Cmatrix HcmtField::field(Fcomp cp, Interval wx, int npx, 
                                   Interval wz, int npz) const
{
	Cmatrix f(npx, npz);
	double dx = (wx.x1-wx.x0)/((double)(npx-1));
	double dz = (wz.x1-wz.x0)/((double)(npz-1));
	double x, z;
	z = wz.x0;
	for(int j=0; j<=npz-1; ++j)
	{
		x = wx.x0;
		for(int i=0; i<=npx-1; ++i) 
		{ 
			f(i, j) = field(cp, x, z); 
			x += dx;
		}
		z += dz;
	}
	return f;
}

/* multiply the entire field by a phase factor exp(i phi),
   the internal representaion is modified accordingly */
void HcmtField::adjustphase(double phi)
{
	Complex a;
	a = expi(phi);
	for(int m=0; m<=nel-1; ++m) (evec[m]).u.multeq(a);
	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 HcmtField::adjustphase(double x, double z)
{
	Fcomp cp = principalcomp(pol);
	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 the window given by xbeg, xend, zbeg, zend */ 
void HcmtField::adjustphase(double xbeg, double xend, double zbeg, double zend,
                            int numx, int numz)
{
	Fcomp cp = principalcomp(pol);
	double x, z;
	double xs, zs;
	xs = (xend-xbeg)/((double) numx);
	zs = (zend-zbeg)/((double) numz);
	x = xbeg;
	double xm = 0.0, zm = 0.0, maxf = 0.0;
	double f;
	for(int xi=0; xi <=numx; ++xi)
	{
		z = zbeg;
		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;
}

/* overlap at position z between the HCMT field and
   mode m, traveling in direction d along the z-axis 
   coarse approximation for fields with pronounced discontinuities;
   the routines assume continuous fields with continuous derivatives */
Complex HcmtField::ovlkernel_x(const Mode& m, Propdir d, 
                               double x, double z) const
{
	Complex k = CC0;
	k += m.cfield(EX, d, x).conj()*field(HY, x, z);
	k += m.cfield(HY, d, x).conj()*field(EX, x, z);
	k -= m.cfield(EY, d, x).conj()*field(HX, x, z);
	k -= m.cfield(HX, d, x).conj()*field(EY, x, z);
	return k;
}
Complex HcmtField::overlap_x(const Mode& m, Propdir d, double z) const
{
	double xr, xm, dx, x;
	static double xc[] = {0.0, 0.1488743389, 0.4333953941,
		              0.6794095682, 0.8650633666, 0.9739065285};
	static double wc[] = {0.0, 0.2955242247, 0.2692667193,
		              0.2190863625, 0.1494513491, 0.0666713443};
	Complex sum = CC0; 
	double a, b, step;
	int numx; 
	Interval ext;
	ext = m.extent(MCROP);
	double xa = ext.x0;
	double xb = ext.x1;
	numx = (int) (ceil((xb-xa)/lambda*nispwl)*10.0);
	step = (xb-xa)/((double) numx);
	Complex s, fval; 
	for(int i=0; i<=numx-1; ++i)
	{
		a  = xa+((double) i)*step;
		b  = a+step;
		xm = 0.5*(b+a);
		xr = 0.5*(b-a);
		s  = CC0;
		for(int j=1; j<=5; j++) 
		{
			dx = xr*xc[j];
			   x  = xm+dx;
			fval  = ovlkernel_x(m, d, x, z);
			   x  = xm-dx;
			fval += ovlkernel_x(m, d, x, z);
			s += wc[j]*fval;
		}
		sum += s*xr;
	}
	return sum/4.0;
}

/* overlap at position x between the HCMT field and
   mode m, traveling in direction d along the x-axis; 
   coarse approximation for fields with pronounced discontinuities,
   the routines assume continuous fields with continuous derivatives */
Complex HcmtField::ovlkernel_z(const Mode& m, Propdir d, 
                               double x, double z) const
{
	Complex k = CC0;
	k += m.vertfield(EY, d, z).conj()*field(HZ, x, z);
	k += m.vertfield(HZ, d, z).conj()*field(EY, x, z);
	k -= m.vertfield(EZ, d, z).conj()*field(HY, x, z);
	k -= m.vertfield(HY, d, z).conj()*field(EZ, x, z);
	return k;
}
Complex HcmtField::overlap_z(const Mode& m, Propdir d, double x) const
{
	double zr, zm, dz, z;
	static double zc[] = {0.0, 0.1488743389, 0.4333953941,
		              0.6794095682, 0.8650633666, 0.9739065285};
	static double wc[] = {0.0, 0.2955242247, 0.2692667193,
		              0.2190863625, 0.1494513491, 0.0666713443};
	Complex sum = CC0; 
	double a, b, step;
	int numz; 
	Interval ext;
	ext = m.extent(MCROP);
	double za = ext.x0;
	double zb = ext.x1;
	numz = (int) (ceil((zb-za)/lambda*nispwl)*10.0);
	step = (zb-za)/((double) numz);
	Complex s, mf, bf, fval; 
	for(int i=0; i<=numz-1; ++i)
	{
		a  = za+((double) i)*step;
		b  = a+step;
		zm = 0.5*(b+a);
		zr = 0.5*(b-a);
		s  = CC0;
		for(int j=1; j<=5; j++) 
		{
			dz = zr*zc[j];
			   z  = zm+dz;
			fval  = ovlkernel_z(m, d, x, z);
			   z  = zm-dz;
			fval += ovlkernel_z(m, d, x, z);
			s += wc[j]*fval;
		}
		sum += s*zr;
	}
	return sum/4.0;
}

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

/* 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 HcmtField::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");
	switch(sttype)
	{
		case REC: 
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
		case OVL: 
			mlout_gengeoxz(dat, ovl, xbeg, xend, zbeg, zend);
		break;
		case MIC:
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
	}
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix cfld = field(cp, Interval(xbeg, xend), npx, 
	                         Interval(zbeg, zend), npz);
	Dmatrix fld = realize(cfld, foa); 
	mlout_fld(dat, npx, npz, cp, fld);
	name[8] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* image 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 HcmtField::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) 
		cmterror("HcmtField::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");
	switch(sttype)
	{
		case REC: 
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
		case OVL: 
			mlout_gengeoxz(dat, ovl, xbeg, xend, zbeg, zend);
		break;
		case MIC:
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
	}
	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;
}

/* 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 HcmtField::movie(double xbeg, double xend, double zbeg, double zend,
                      int npx, int npz, int ntfr, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "an____M.m";
	double tmax, dt;
	Fcomp cp = principalcomp(pol);
        name[2] = fldchr(cp);
        name[3] = cpchr(cp);      
	name[4] = ext0;
	name[5] = ext1;
	tmax = val_period_T(lambda);
	dt = tmax/ntfr;
	fprintf(stderr, 
	   "movie(%d x [%g, %g] x [%g, %g]), T=%g fs, dt=%g fs >> %s\n", 
	   ntfr, xbeg, xend, zbeg, zend, tmax, dt, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	switch(sttype)
	{
		case REC: 
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
		case OVL: 
			mlout_gengeoxz(dat, ovl, xbeg, xend, zbeg, zend);
		break;
		case MIC:
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
	}
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix fld = field(cp, Interval(xbeg, xend), npx, 
	                        Interval(zbeg, zend), npz); 
	double famp = fld.max();
	mlout_fld(dat, npx, npz, cp, fld.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npz, cp, fld.im());
	mlout_fldtoim(dat, cp);
	name[7] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_genmovie(cp, name, desc, dat, ntfr, dt, val_omega(lambda), famp); 
	fclose(dat);
	return;
}
		

/* export full field data 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 HcmtField::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, pol, lambda, vM);
	switch(sttype)
	{
		case REC: 
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
		case OVL: 
			mlout_gengeoxz(dat, ovl, xbeg, xend, zbeg, zend);
		break;
		case MIC:
			mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
		break;
	}
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix f;
	Fcomp cp;
	if(vM)
	{
		cp = EX;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EZ;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HX;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HY;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HZ;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);

	}
	else 
	{
		if(pol == TE)
		{
			cp = EX;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
			cp = EY;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);
			cp = EZ;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
			cp = HX;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);
			cp = HY;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
			cp = HZ;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);

		}
		else 	
		{
			cp = EX;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);
			cp = EY;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
			cp = EZ;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);
			cp = HX;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
			cp = HY;
			f = field(cp, Interval(xbeg, xend), npx, 
				      Interval(zbeg, zend), npz); 
			mlout_fld(dat, npx, npz, cp, f.re());
			mlout_fldtore(dat, cp);
			mlout_fld(dat, npx, npz, cp, f.im());
			mlout_fldtoim(dat, cp);
			cp = HZ;
			mlout_0fld(dat, cp);
			mlout_fldtore(dat, cp);
			mlout_0fld(dat, cp);
			mlout_fldtoim(dat, cp);
		}
	}

	Dmatrix n(npx, npz);
	double dx = (xend-xbeg)/((double)(npx-1));
	double dz = (zend-zbeg)/((double)(npz-1));
	double x, z;
	z = zbeg;
	for(int j=0; j<=npz-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) { n(i, j) = sqrt(eps(x, z)); x += dx; }
		z += dz;
	}
	mlout_n(dat, npx, npz, n);
	mlout_fldviewer(dat, name);
	fclose(dat);
	return;
}

/* --- vHCMT variant -------------------------------------------------- */

/* initialize, sttype = OVL */
HcmtField::HcmtField(OvlStruct s, double gky,  double wl,
                     Interval Cwx, Interval Cwz, int Nispwl)
{
	(*this) = HcmtField(s, TE, wl, Cwx, Cwz, Nispwl);
	vM = 1;
	ky = gky;
}

/* initialization, guided mode supported by a circular step-index fiber 
   t: FIM
   fim: the corresponding fiber mode object */
HcmtElement::HcmtElement(HcmtBft t, FIMode fim)
{
	switch(t)
	{
	case FIM: 
		type = t;
		break;
	default: 
		cmterror("HcmtElement constructor, fiber mode");
	}
	fimode = fim;
	vM = 1;
	pol = TE;	
	k0 = val_k0(fim.wg.lambda);
	invommu0 = val_invommu0(fim.wg.lambda);
	invomep0 = val_invomep0(fim.wg.lambda);
	x0 = fim.x0;
	z0 = fim.y0;
	u.strip();
	seg = NULL;
	Ns = 0;
	u = Cvector(Ns+1);
	u.init(CC0);
	pramp.strip();
	prpos.strip();
}
/* ... in an HcmtBundle: observe amplitude evolution 
   at ny+1 equidistant positions ya <= y <= yb */
HcmtElement::HcmtElement(HcmtBft t, FIMode fim, double ya, double yb, int ny)
{
	(*this) = HcmtElement(t, fim);
	double s = (yb-ya)/((double) ny);
	pramp = Cvector(ny+1); 
	prpos = Dvector(ny+1);
	pramp.init(CC0);
	for(int i=0; i<=ny; ++i) prpos(i) = ya + ((double) i)*s;
}

/* waveguide elements, normalized slab modes, 
   potentially non-real kz: output power */
double HcmtElement::pout()
{
	switch(type)
	{	
	case MHF: case MHB: case MVU: case MVD:
		if(vM)
		{
			if(fabs(d_kz.im) <= 1.0e-100) return aout.sqabs();
			return 0.0;
		}
		return aout.sqabs();
		break;
	default: 
		cmterror("HcmtElement::pout: non-waveguide element");
		return 0.0;
		break;
	}
	return 0.0;
}

/* change wavenumber parameter */
void HcmtElement::setky(double gky)
{
	vM = 1;
	ky = gky;
	switch(type)
	{
	case MHF: case MHB: case MVU: case MVD:
		mode.vectorify(Complex(gky, 0.0));
		if(dir == FORW)
		{
			d_kz = mode.kz;
			d_beta  =  mode.kz.re;
			d_alpha = -mode.kz.im;
		}
		else 
		{
			d_kz = -mode.kz;
			d_beta  = -mode.kz.re;
			d_alpha =  mode.kz.im;
		}
		break;
	case FIM:
		break;
	default:
		cmterror("HcmtElement::setky: vM, element type");
		break;
	}
	return;
}

/* change global wavenumber parameter */
void HcmtField::setky(double gky)
{
	vM = 1;
	ky = gky;
	for(int j=0; j<=nel-1; ++j) (evec[j]).setky(gky);
	return;
}

/* ------------------------------------------------------------------------   
   - Gaussian bundles of vectorial HCMT solutions                         - 
   ------------------------------------------------------------------------ */

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

/* ------------------------------------------------------------------- 
 * bundles of vHCMT solutions
 */

/* full setup:
      pre-filled field object f;
   excitation:  
      a guided mode, HcmtElement eli, 
      coming in with unit amplitude, 
      at primary angle of incidence t (^o), 
      beam of full width w (1/e, field),
      represented by na angles between t-dt and t+dt,
      dt determined such that spectral contributions with relative 
      amplitude below dtlev are neglected; 
      positioned at focus coordinates fcy, fcz 
   numerical integration needs to resolve a resonance 
      in angular range [tr-dtr, tr+dtr] with nr equidisttant steps */
vHcmtBundle::vHcmtBundle(HcmtField f, 
                         int eli, 
                         double t, double w, double fcy, double fcz, int na, double dtlev,
			 double tr, double dtr, int nr)
{
	if(f.sttype != OVL) cbundleerror("vHcmtBundle::vHcmtBundle(): sttype not supported");
	fld = f;
	str = fld.ovl;
	k0 = val_k0(str.lambda);

	if(fld(eli).type != MHF) cbundleerror("vHcmtBundle::vHcmtBundle(): inputelement type not supported");
	inel = eli;
	Nin = fld(inel).mode.neff;

	theta0 = t;
	bhw = w/2.0;
	y0 = fcy;
	z0 = fcz;

	Nky = na;
	Wky = 2.0/bhw;
	Dky = Wky*sqrt(-log(dtlev));


	cwx = fld.cwx;
	cwz = fld.cwz;
	nispwl = fld.nispwl;

	ky0 = k0*Nin*sin(theta0/180.0*PI);

	kyR0 = k0*Nin*sin((tr-dtr)/180.0*PI);
	kyR1 = k0*Nin*sin((tr+dtr)/180.0*PI);
	NkyR = nr;

fprintf(stderr, "\n------------- vHCMT-bundle ----------------------------------------- '20 ---\n");
fprintf(stderr, "Nin = %.10g, lambda = %.10g mum, k0 = %.10g /mum\n", Nin, str.lambda, k0);
fprintf(stderr, "y0 = %g mum  ", y0);
fprintf(stderr, "z0 = %g mum\n", z0);
fprintf(stderr, "width_y = %g mum  ", 2.0*bhw);
fprintf(stderr, "width_cr = %g mum\n", 2.0*bhw*cos(theta0*PI/180.0));
fprintf(stderr, "theta0 = %g^o  ", theta0);
fprintf(stderr, "ky0 = %g /mum\n", ky0);
fprintf(stderr, "ky /mum: [%g :: %g : %g : %g :: %g]_%d\n", 
	ky0-Dky, ky0-Wky, ky0, ky0+Wky, ky0+Dky, Nky);
if(fabs((ky0-Dky)/k0/Nin)<1 && fabs((ky0+Dky)/k0/Nin)<1) 
{
	fprintf(stderr, "theta /^o: [%g :: %g : %g : %g :: %g]\n", 
	asin((ky0-Dky)/k0/Nin)/PI*180.0, 
	asin((ky0-Wky)/k0/Nin)/PI*180.0, 
	asin((ky0    )/k0/Nin)/PI*180.0, 
	asin((ky0+Wky)/k0/Nin)/PI*180.0, 
	asin((ky0+Dky)/k0/Nin)/PI*180.0);
}
else
{
	fprintf(stderr, "ky_max = %g /mum\n", k0*Nin);
	cbundleerror("vHcmtBundle::vHcmtBundle(): unphysical spectral range");
}
if(NkyR > 1)
{
	fprintf(stderr, "Res ky /mum: [%g, %g]_%d\n", kyR0, kyR1, NkyR);
	fprintf(stderr, "  theta /^o: [%g, %g]\n", tr-dtr, tr+dtr);
}
fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "x in (%g, %g)  z in (%g, %g)  nispwl: %d\n", cwx.x0, cwx.x1, cwz.x0, cwz.x1, nispwl);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	str.write(stderr);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	VA.clear();
	FA.clear();
	PA.clear();
	Pin = 0.0;
	eliR.strip();
	dkyR = 0.0;
}
/* ... no resonance expected */
vHcmtBundle::vHcmtBundle(HcmtField f, 
                         int eli, 
                         double t, double w, double fcy, double fcz, int na, double dtlev) 
{
	(*this) = vHcmtBundle(f, eli, t, w, fcy, fcz, na, dtlev, 0.0, 0.0, 0); 	
}
/* ...  numerical integration needs to resolve resonances in 
   the output amplitudes of elements with indices in reli, down to an angular stepsize dtr */
vHcmtBundle::vHcmtBundle(HcmtField f, 
                         int eli, 
                         double t, double w, double fcy, double fcz, int na, double dtlev,
	                 Ivector reli, double dtr) 
{
	(*this) = vHcmtBundle(f, eli, t, w, fcy, fcz, na, dtlev, 0.0, 0.0, 0); 	
fprintf(stderr, "\n");
	eliR = reli;
	dkyR = fabs(k0*Nin*sin((t+dtr/2.0)/180.0*PI) - k0*Nin*sin((t-dtr/2.0)/180.0*PI));
fprintf(stderr, "ResResolve: d_theta = %g^o, dky = %g /mum\n", dtr, dkyR);
fprintf(stderr, "   Elements: ");
for(int i=0; i<=eliR.nel-1; ++i) fprintf(stderr, " %d", eliR(i));
fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "\n");
}


/* field for a single lateral spectral component */
Complex vHcmtBundle::spfld(Fcomp cp, double x, double y, double z)
{
	return fld.field(cp, x, z)*expmi(ky*(y-y0));
}

/* stepwise vHCMT: fill views with field data */
void vHcmtBundle::solve()
{
	for(int n=0; n<=FA.num-1; ++n) FA(n).eps = str.eps(FA(n).x, FA(n).z);
	if(eliR.nel > 0)
	{
		solveR();
		return;
	}
	double s = 2.0*Dky/Nky; 
	if(NkyR <= 1)
	{
		solve(ky0-Dky, ky0+Dky, Nky);
		return;
	}
	if(kyR1 < ky0-Dky || kyR0 > ky0+Dky)
	{
		fprintf(stderr, "\n(?) Resonance outside of beam wavenumber range.\n\n");
		solve(ky0-Dky, ky0+Dky, Nky);
		solve(kyR0, kyR1, NkyR);
		return;
	}
	if(kyR0 <= ky0-Dky && ky0+Dky <= kyR1)
	{
		fprintf(stderr, "\n(?) Resonance wider than beam wavenumber range.\n");
		solve(kyR0, kyR1, NkyR);
		return;
	}
	if(kyR0 <= ky0-Dky && ky0-Dky <= kyR1 && kyR1 <= ky0+Dky)
	{
		fprintf(stderr, "\n(?) Resonance at left border of beam wavenumber range.\n");
		solve(kyR0, kyR1, NkyR);
		solve(kyR1, ky0+Dky, (int) ceil(((ky0+Dky)-kyR1)/s));
		return;
	}
	if(ky0-Dky <= kyR0 && kyR0 <= ky0+Dky && ky0+Dky <= kyR1)
	{
		fprintf(stderr, "\n(?) Resonance at right border of beam wavenumber range.\n");
		solve(ky0-Dky, kyR0, (int) ceil((kyR0-(ky0-Dky))/s));
		solve(kyR0, kyR1, NkyR);
		return;
	}
	if(ky0-Dky <= kyR0 && kyR1 <= ky0+Dky)
	{
		fprintf(stderr, "\n(!) Resonance inside beam wavenumber range.\n");
		solve(ky0-Dky, kyR0, (int) ceil((kyR0-(ky0-Dky))/s));
		solve(kyR0, kyR1, NkyR);
		solve(kyR1, ky0+Dky, (int) ceil(((ky0+Dky)-kyR1)/s));
		return;
	}
	cbundleerror("vHcmtBundle: solve(resonance): confused");
	return;
}

/* ... proceed from k0 to k1 in (roughly) nk steps */
void vHcmtBundle::solve(double ka, double kb, int nk)
{
	if(nk <= 6) nk = 7;
	fprintf(stderr, "\n +++   ky /mum = [%g, %g]_%d\n", ka, kb, nk); 
	fprintf(stderr, " +++ theta /^o = [%g, %g]\n", asin(ka/k0/Nin)/PI*180.0, asin(kb/k0/Nin)/PI*180.0);
	double kys = (kb-ka)/((double) nk);
	double A = (1.0/sqrt(2.0*PI*sqrt(PI/2.0)))/sqrt(Wky);
	for(int sci=0; sci<=nk; ++sci)
	{
		ky = ka+((double) sci)*kys;
		double theta = asin(ky/k0/Nin)/PI*180.0;
		double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
		// Complex amp = expmi(kz*(fld(inel).v_z0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		Complex amp = A*expmi(kz*(fld(inel).v_z0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		fprintf(stderr, "\n   ++[%d] theta = %g ^o, amp = %g\n", sci, theta, amp.abs()); 

		fld.setky(ky);
		fld(inel).ain = amp;
		fld.solve();

		double p = 0.0;
		for(int j=0; j<=fld.nel-1; ++j)
		{
			HcmtBft t = fld(j).type;
			if(t == MHF || t == MHB || t == MVU || t == MVD) p += fld(j).pout();
		}
		fprintf(stderr, "P_out/P_in = %g\n", p/amp.sqabs());

		double                      icoeff =       1.0*kys;
		if(sci == 0 || sci == nk  ) icoeff =   3.0/8.0*kys;
		if(sci == 1 || sci == nk-1) icoeff =   7.0/6.0*kys; 
		if(sci == 2 || sci == nk-2) icoeff = 23.0/24.0*kys; 

		for(int n=0; n<=VA.num-1; ++n)
		{ 
			double dv = (VA(n).va_end-VA(n).va_beg)
			            /((double) (VA(n).va_N-1));
			double dh = (VA(n).ha_end-VA(n).ha_beg)
			            /((double) (VA(n).ha_N-1));
			for(int i=0; i<=VA(n).va_N-1; ++i)
			{
				double v = VA(n).va_beg+((double) i)*dv;
				for(int j=0; j<=VA(n).ha_N-1; ++j)
				{
					double h = VA(n).ha_beg+((double) j)*dh;
					double x=0.0;
					double y=0.0;
					double z=0.0;
					VA(n).trf(v, h, x, y, z);
					VA(n).ex(i, j) += spfld(EX, x, y, z)*icoeff;
					VA(n).ey(i, j) += spfld(EY, x, y, z)*icoeff;
					VA(n).ez(i, j) += spfld(EZ, x, y, z)*icoeff;
					VA(n).hx(i, j) += spfld(HX, x, y, z)*icoeff;
					VA(n).hy(i, j) += spfld(HY, x, y, z)*icoeff;
					VA(n).hz(i, j) += spfld(HZ, x, y, z)*icoeff;
				}
			}
		}
		for(int n=0; n<=FA.num-1; ++n)
		{  
			double x = FA(n).x;
			double y = FA(n).y;
			double z = FA(n).z;
			FA(n).ex += spfld(EX, x, y, z)*icoeff;
			FA(n).ey += spfld(EY, x, y, z)*icoeff;
			FA(n).ez += spfld(EZ, x, y, z)*icoeff;
			FA(n).hx += spfld(HX, x, y, z)*icoeff;
			FA(n).hy += spfld(HY, x, y, z)*icoeff;
			FA(n).hz += spfld(HZ, x, y, z)*icoeff;
		}
		for(int n=0; n<=PA.num-1; ++n)
		{  
			Complex a = fld(PA(n).eli).aout;
			PA(n).Pout += a.sqabs()*icoeff;
		}
		Pin += amp.sqabs()*icoeff;
for(int i=0; i<=fld.nel-1; ++i)
{
	if(fld(i).type == FIM) 
	{
		Complex u = fld(i).u(0);
		fprintf(stderr, "%s: amp = %g +i %g, |amp|^2 = %g, icont = %g\n", fld(i).fimode.ids, u.re, u.im, u.sqabs(), u.sqabs()*icoeff);
// apptoxyf(fld(i).fimode.sids, theta, u.sqabs());
	}
}
	}
	return;
}

/* ... resolve amplitude resonances */
void vHcmtBundle::solveR()
{
	double ka = ky0-Dky;
	double kb = ky0+Dky;
	double kys = (kb-ka)/((double) Nky);
	fprintf(stderr, "\n ++E   ky /mum = [%g, %g]_%d\n", ka, kb, Nky); 
	fprintf(stderr, " ++E theta /^o = [%g, %g]\n", asin(ka/k0/Nin)/PI*180.0, asin(kb/k0/Nin)/PI*180.0);

	Dvector knp, anp;
	knp.strip();
	anp.strip();

	Dvector th, am;
	th.strip();
	am.strip();
	int done=0;
	while(done==0)
	{
		done = next(Interval(ka, kb), kys, dkyR, th, am);
		if(done==0)
		{
			ky = th(th.nel-1);
			double theta = asin(ky/k0/Nin)/PI*180.0;
			fprintf(stderr, "\n   ++[E] theta = %g ^o\n", theta); 

			fld.setky(ky);
			fld(inel).ain = CC1;
			fld.solve();
	
			double amax = 0.0;
			for(int j=0; j<=eliR.nel-1; ++j)
			{ 
				Complex u = fld(eliR(j)).aout;
				if(u.sqabs() > amax) amax = u.sqabs();
			}
			fprintf(stderr, "\n   ++[E]R theta = %g ^o, aout_max = %g\n", theta, amax); 
			am.append(amax);
		
			knp.append(ky);
			anp.append(amax);
		}
	}
	sort(knp, anp);
// toxyf("expky0", knp, anp);
	double amax = anp.max();
	double da = amax/Nky;
	Dvector knpn;
	knpn.strip();
// Dvector anpn;
// anpn.strip();
	for(int i=0; i<=knp.nel-2; ++i)
	{
		if(fabs(anp(i+1)-anp(i)) > da)
		{
			int nap = floor(fabs(anp(i+1)-anp(i))/da);
			double s = (knp(i+1)-knp(i))/((double)(nap+1));
			for(int j=1; j<=nap; ++j)
			{
				knpn.append(knp(i)+((double)j)*s);
// anpn.append(-1.0);
			}
		}
	}
	Dvector knpf(knp.nel+knpn.nel);
// Dvector anpf(anp.nel+anpn.nel);
	if(knpn.nel > 0)
	{
		knpf.setsubvector(knp, 0);
		knpf.setsubvector(knpn, knp.nel);
// anpf.setsubvector(anp, 0);
// anpf.setsubvector(anpn, anp.nel);
	}
	else 
	{
		knpf = knp;
// anpf = anp;
	}
// sort(knpf, anpf);
// toxyf("expky", knpf, anpf);

	for(int i=0; i<=fld.nel-1; ++i)
	{
		if(fld(i).type == FIM) fld(i).pramp.init(CC0);
	}
	
	knpf.sort();
	int nk = knpf.nel;
	fprintf(stderr, "\n +++   ky /mum = [%g, %g]_%d\n", knpf(0), knpf(nk-1), nk); 
	fprintf(stderr, " +++ theta /^o = [%g, %g]\n", asin(knpf(0)/k0/Nin)/PI*180.0, asin(knpf(nk-1)/k0/Nin)/PI*180.0);
	double A = (1.0/sqrt(2.0*PI*sqrt(PI/2.0)))/sqrt(Wky);
	for(int ik=0; ik<=nk-1; ++ik)
	{
		ky = knpf(ik);
		double theta = asin(ky/k0/Nin)/PI*180.0;
		double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
		// Complex amp = expmi(kz*(fld(inel).v_z0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		Complex amp = A*expmi(kz*(fld(inel).v_z0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		fprintf(stderr, "\n   ++[%d] theta = %g ^o, amp = %g\n", ik, theta, amp.abs()); 

		fld.setky(ky);
		fld(inel).ain = amp;
		fld.solve();

		double p = 0.0;
		for(int j=0; j<=fld.nel-1; ++j)
		{
			HcmtBft t = fld(j).type;
			if(t == MHF || t == MHB || t == MVU || t == MVD) p += fld(j).pout();
		}
		fprintf(stderr, "P_out/P_in = %g\n", p/amp.sqabs());

		double icoeff = 0.0;
		if(ik == 0)              icoeff = 0.5*(knpf(1)-knpf(0));
		if(ik >=1 && ik <= nk-2) icoeff = 0.5*(knpf(ik+1)-knpf(ik-1));
		if(ik == nk-1)           icoeff = 0.5*(knpf(nk-1)-knpf(nk-2));

		for(int n=0; n<=VA.num-1; ++n)
		{ 
			double dv = (VA(n).va_end-VA(n).va_beg)
			            /((double) (VA(n).va_N-1));
			double dh = (VA(n).ha_end-VA(n).ha_beg)
			            /((double) (VA(n).ha_N-1));
			for(int i=0; i<=VA(n).va_N-1; ++i)
			{
				double v = VA(n).va_beg+((double) i)*dv;
				for(int j=0; j<=VA(n).ha_N-1; ++j)
				{
					double h = VA(n).ha_beg+((double) j)*dh;
					double x=0.0;
					double y=0.0;
					double z=0.0;
					VA(n).trf(v, h, x, y, z);
					VA(n).ex(i, j) += spfld(EX, x, y, z)*icoeff;
					VA(n).ey(i, j) += spfld(EY, x, y, z)*icoeff;
					VA(n).ez(i, j) += spfld(EZ, x, y, z)*icoeff;
					VA(n).hx(i, j) += spfld(HX, x, y, z)*icoeff;
					VA(n).hy(i, j) += spfld(HY, x, y, z)*icoeff;
					VA(n).hz(i, j) += spfld(HZ, x, y, z)*icoeff;
				}
			}
		}
		for(int n=0; n<=FA.num-1; ++n)
		{  
			double x = FA(n).x;
			double y = FA(n).y;
			double z = FA(n).z;
			FA(n).ex += spfld(EX, x, y, z)*icoeff;
			FA(n).ey += spfld(EY, x, y, z)*icoeff;
			FA(n).ez += spfld(EZ, x, y, z)*icoeff;
			FA(n).hx += spfld(HX, x, y, z)*icoeff;
			FA(n).hy += spfld(HY, x, y, z)*icoeff;
			FA(n).hz += spfld(HZ, x, y, z)*icoeff;
		}
		for(int n=0; n<=PA.num-1; ++n)
		{  
			Complex a = fld(PA(n).eli).aout;
			PA(n).Pout += a.sqabs()*icoeff;
		}
		Pin += amp.sqabs()*icoeff;

		for(int i=0; i<=fld.nel-1; ++i)
		{
			if(fld(i).type == FIM) 
			{
				for(int j=0; j<=fld(i).prpos.nel-1; ++j)
				{
					double yp = fld(i).prpos(j);
					Complex a = fld(i).aout*expmi(ky*(yp-y0));
					fld(i).pramp(j) += a*icoeff;
				}
			}
		}

for(int i=0; i<=fld.nel-1; ++i)
{
	if(fld(i).type == FIM) 
	{
		Complex u = fld(i).u(0);
		fprintf(stderr, "%s: amp = %g +i %g, |amp|^2 = %g, icont = %g\n", fld(i).fimode.ids, u.re, u.im, u.sqabs(), u.sqabs()*icoeff);
// apptoxyf(fld(i).fimode.sids, theta, u.sqabs());
	}
}

	}
	return;
}

/* ------------------------------------------------------------------- 
 * containers for registration of guided modal output power 
 */

/* initialize */
cPort::cPort()
{
	eli = -1;
	id1 = '?';
	id2 = '?';
	Pout = 0.0;
}
cPort::cPort(int e, char c1, char c2)
{
	eli = e;
	id1  = c1;
	id2  = c2;
	Pout = 0.0;
}
   
/* ------------------------------------------------------------------- 
 * an array of containers for registration of guided modal output power 
 */
   
/* initialize */
cPortArray::cPortArray()
          : num(0)
{
	pvec = NULL;
}

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

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

/* delete all cPort entries */
void cPortArray::clear()
{
	if(pvec != NULL) delete[] pvec;
	pvec = NULL;
	num = 0; 
}

/* subscripting */
cPort& cPortArray::operator() (int i)  
{ 
	if(i<0 || i>=num) cbundleerror("cPortArray: () out of range");
	return pvec[i]; 
}
cPort cPortArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) cbundleerror("cPortArray: () out of range");
	return pvec[i]; 
}

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

/* ------------------------------------------------------------------- 
 * an array of containers for discretized fields 
 */
   
/* initialize */
cViewArray::cViewArray()
          : num(0)
{
	vvec = NULL;
}

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

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

/* delete all cView entries */
void cViewArray::clear()
{
	if(vvec != NULL) delete[] vvec;
	vvec = NULL;
	num = 0; 
}

/* subscripting */
cView& cViewArray::operator() (int i)  
{ 
	if(i<0 || i>=num) cbundleerror("cViewArray: () out of range");
	return vvec[i]; 
}
cView cViewArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) cbundleerror("cViewArray: () out of range");
	return vvec[i]; 
}

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

/* ------------------------------------------------------------------- 
 * containers for discretized fields, plot routines 
 */

/* initialize */
cView::cView()
{
	type = cXZ;
	ha_beg = -1.0;;
	ha_end = 1.0;
	ha_N   = 10;
	va_beg = -1.0;
	va_end = 1.0;
	va_N   = 10;
	pos = 0.0;
	ex = Cmatrix(va_N, ha_N); ex.init(CC0);
	ey = Cmatrix(va_N, ha_N); ey.init(CC0);
	ez = Cmatrix(va_N, ha_N); ez.init(CC0);
	hx = Cmatrix(va_N, ha_N); hx.init(CC0);
	hy = Cmatrix(va_N, ha_N); hy.init(CC0);
	hz = Cmatrix(va_N, ha_N); hz.init(CC0);
	va_chr = 'x';
	ha_chr = 'z';
	st = OvlStruct();
	pid0 = '-';
	pid1 = '-';
}

/* initialize: supply discretization information, setup matrices */
cView::cView(cVtype t, OvlStruct s, 
           double v_beg, double v_end, int v_N, 
           double h_beg, double h_end, int h_N, double p, 
           char ext0, char ext1)
{
	type = t;
	ha_beg = h_beg;
	ha_end = h_end;
	ha_N   = h_N;
	va_beg = v_beg;
	va_end = v_end;
	va_N   = v_N;
	pos = p;
	ex = Cmatrix(va_N, ha_N); ex.init(CC0);
	ey = Cmatrix(va_N, ha_N); ey.init(CC0);
	ez = Cmatrix(va_N, ha_N); ez.init(CC0);
	hx = Cmatrix(va_N, ha_N); hx.init(CC0);
	hy = Cmatrix(va_N, ha_N); hy.init(CC0);
	hz = Cmatrix(va_N, ha_N); hz.init(CC0);
	va_chr = 'v';
	ha_chr = 'h';
	switch(type)
	{
	case cXZ: va_chr = 'x'; ha_chr = 'z'; break;
	case cYZ: va_chr = 'y'; ha_chr = 'z'; break;
	case cXY: va_chr = 'x'; ha_chr = 'y'; break;
	}
	Waveguide wg;
	double d, r, dx, dz, ri, di; 
	Overlay oo, on;
	switch(type)
	{
	case cXZ: 
		st = s;
		break;
	case cYZ: 
		st = OvlStruct(s.bgeps,s.lambda);
		for(int j=0; j<=s.no-1; ++j)
		{
			oo = s(j);
			switch(oo.type)
			{
			case HLAY:
				if(pos >= oo.lbd && pos <= oo.ubd)
				{
					on = Overlay(VLAY, oo.n, 2.0*AWAY, 0.0);
					st.place(on);	
				}
				break;
			case VLAY:
				on = Overlay(VLAY, oo.n, oo.ubd-oo.lbd, 0.5*(oo.ubd+oo.lbd));
				st.place(on);	
				break;
			case RING:
				r = oo.ubd;
				dx = pos-oo.xo;
				if(fabs(dx) < r)
				{
					d = 2.0*sqrt(r*r-dx*dx);
					ri = oo.lbd;
					if(fabs(dx) > ri)
					{
						on = Overlay(VLAY, oo.n, d, oo.zo);	
						st.place(on);	
					}
					else 
					{
						di = 2.0*sqrt(ri*ri-dx*dx);
						on = Overlay(VLAY, oo.n, 0.5*(d-di), oo.zo-0.5*(d+di));	
						st.place(on);	
						on = Overlay(VLAY, oo.n, 0.5*(d-di), oo.zo+0.5*(d+di));	
						st.place(on);	
					}
				}
				break;
			case DISK:
				r = oo.ubd;
				dx = pos-oo.xo;
				if(fabs(dx) < r)
				{
					d = 2.0*sqrt(r*r-dx*dx);
					on = Overlay(VLAY, oo.n, d, oo.zo);	
					st.place(on);	
				}
				break;
			case RECT:
				if(pos >= oo.lbd && pos <= oo.ubd)
				{
					on = Overlay(VLAY, oo.n, oo.zo-oo.xo, 0.5*(oo.zo+oo.xo));
					st.place(on);	
				}
				break;
			default: 
				cbundleerror("cView(): YZ, struct");
				break;
			}
		}
		break;
	case cXY: 
		st = OvlStruct(s.bgeps,s.lambda);
		for(int j=0; j<=s.no-1; ++j)
		{
			oo = s(j);
			switch(oo.type)
			{
			case VLAY:
				if(pos >= oo.lbd && pos <= oo.ubd)
				{
					on = Overlay(HLAY, oo.n, 2.0*AWAY, 0.0);
					st.place(on);	
				}
				break;
			case HLAY:
				on = Overlay(HLAY, oo.n, oo.ubd-oo.lbd, 0.5*(oo.ubd+oo.lbd));
				st.place(on);	
				break;
			case RING:
				r = oo.ubd;
				dz = pos-oo.zo;
				if(fabs(dz) < r)
				{
					d = 2.0*sqrt(r*r-dz*dz);
					ri = oo.lbd;
					if(fabs(dz) > ri)
					{
						on = Overlay(HLAY, oo.n, d, oo.xo);	
						st.place(on);	
					}
					else 
					{
						di = 2.0*sqrt(ri*ri-dz*dz);
						on = Overlay(HLAY, oo.n, 0.5*(d-di), oo.xo-0.5*(d+di));	
						st.place(on);	
						on = Overlay(HLAY, oo.n, 0.5*(d-di), oo.xo+0.5*(d+di));	
						st.place(on);	
					}
				}
				break;
			case DISK:
				r = oo.ubd;
				dz = pos-oo.zo;
				if(fabs(dz) < r)
				{
					d = 2.0*sqrt(r*r-dz*dz);
					on = Overlay(HLAY, oo.n, d, oo.xo);	
					st.place(on);	
				}
				break;
			case RECT:
				if(pos >= oo.xo && pos <= oo.zo)
				{
					on = Overlay(HLAY, oo.n, oo.ubd-oo.lbd, 0.5*(oo.ubd+oo.lbd));
					st.place(on);	
				}
				break;
			default: 
				cbundleerror("cView(): XY, struct");
				break;
			}
		}
		break;
	}
	pid0 = ext0;
	pid1 = ext1;
}

/* transform v, h, pos to x, y, z coordinates */
void cView::trf(double v, double h, double& x, double& y, double& z) const
{
	switch(type)
	{
	case cXZ: 
		x = v; y = pos; z = h; return;
		break;
	case cYZ: 
		x = pos; y = v; z = h; return;
		break;
	case cXY: 
		x = v; y = h; z = pos; return;
		break;
	}
	return;
}

/* visualization: output to MATLAB .m files 
 * ext0, ext1: filename id characters for the m.file */

/* image plot corresponding to field component cp
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   foa: ORG, MOD, SQR, REP, IMP (ORG = REP) */
void cView::plot(Fcomp cp, Afo foa) const 
{
	FILE *dat;
	char name[15] = "__pl_____I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[4] = afochr(foa);
        name[5] = fldchr(cp);
        name[6] = cpchr(cp);      
	name[7] = pid0;
	name[8] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld;
	switch(cp)
	{
		case EX: fld = realize(ex, foa); break;
		case EY: fld = realize(ey, foa); break;
		case EZ: fld = realize(ez, foa); break;
		case HX: fld = realize(hx, foa); break;
		case HY: fld = realize(hy, foa); break;
		case HZ: fld = realize(hz, foa); break;
		default: cbundleerror("cView::plot: cp"); break;
	}
	mlout_fld(dat, va_N, ha_N, cp, fld);
	name[10] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(va_chr, ha_chr, cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* image plot of the field "strength", id: one of mE, mH, qE, qH, mW */
void cView::plot(FSid id) const
{
	if(va_N <= 1 || ha_N <= 1) 
		cbundleerror("cView::plot, FSid: va_N, ha_N"); 
	FILE *dat;
	char name[15] = "__pl____I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[4] = mqchr(id);
	name[5] = idchr(id);
	name[6] = pid0;
	name[7] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld(va_N, ha_N);
	double eps;
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) 
		{
			switch(id)
			{
			case mE:
				fld(i, j) = sqrt( ex(i, j).sqabs()
					         +ey(i, j).sqabs()
					         +ez(i, j).sqabs());
				break;
			case mH:
				fld(i, j) = sqrt( hx(i, j).sqabs()
					         +hy(i, j).sqabs()
					         +hz(i, j).sqabs());
				break;
			case qE:
				fld(i, j) = ( ex(i, j).sqabs()
					     +ey(i, j).sqabs()
					     +ez(i, j).sqabs());
				break;
			case qH:
				fld(i, j) = ( hx(i, j).sqabs()
					     +hy(i, j).sqabs()
					     +hz(i, j).sqabs());
				break;
			case mW:
				eps = st.eps(vp, hp);
				fld(i, j) = 0.25*(    MU0*(hx(i, j).sqabs()
							  +hy(i, j).sqabs()
							  +hz(i, j).sqabs())
						 +EP0*eps*(ex(i, j).sqabs()
							  +ey(i, j).sqabs()
							  +ez(i, j).sqabs()));
				break;
			default:
				cbundleerror("plot, FSid?");
				break;
			}
			vp += dv;
		}
		hp += dh;
	}
	mlout_fld(dat, va_N, ha_N, SZ, fld);
	name[9] = 0; 
	char desc[10];
	iddesc(id, desc);
	mlout_genimage(va_chr, ha_chr, SZ, MOD, name, desc, dat); 
	mlout_lowlevcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* electromagnetic energy density image */
void cView::plot() const
{
	if(va_N <= 1 || ha_N <= 1) 
		cbundleerror("cView::plot, energy density: va_N, ha_N"); 
	FILE *dat;
	char name[15] = "__plW__I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[5] = pid0;
	name[6] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld(va_N, ha_N);
	double eps;
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) 
		{
			
			eps = st.eps(vp, hp);
			fld(i, j) = 0.25*(    MU0*(hx(i, j).sqabs()
			                          +hy(i, j).sqabs()
			                          +hz(i, j).sqabs())
                                         +EP0*eps*(ex(i, j).sqabs()
			                          +ey(i, j).sqabs()
			                          +ez(i, j).sqabs()));
			vp += dv;
		}
		hp += dh;
	}
	mlout_fld(dat, va_N, ha_N, SZ, fld);
	name[8] = 0; 
	char desc[3];
	desc[0] = 'W'; desc[1] = 0;
	mlout_genimage(va_chr, ha_chr, SZ, MOD, name, desc, dat); 
	mlout_lincolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* animation of the light propagation
   the frames show images of component cp
   at ntfr equidistant times over one time period */
void cView::movie(Fcomp cp, int ntfr) const
{
	FILE *dat;
	char name[15] = "__an____M.m";
	double tmax, dt;
        name[0] = va_chr;
        name[1] = ha_chr;
        name[4] = fldchr(cp);
        name[5] = cpchr(cp);      
	name[6] = pid0;
	name[7] = pid1;
	tmax = val_period_T(st.lambda);
	dt = tmax/ntfr;
	fprintf(stderr, 
	   "movie_%c%c(%d x [%g, %g] x [%g, %g] @ %g), T=%g fs, dt=%g fs >> %s\n", 
	   va_chr, ha_chr, ntfr, va_beg, va_end, ha_beg, ha_end, pos, tmax, dt, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Cmatrix fld;
	switch(cp)
	{
		case EX: fld = ex; break;
		case EY: fld = ey; break;
		case EZ: fld = ez; break;
		case HX: fld = hx; break;
		case HY: fld = hy; break;
		case HZ: fld = hz; break;
		default: cbundleerror("cView::movie: cp"); break;
	}
	double famp = fld.max();
	mlout_fld(dat, va_N, ha_N, cp, fld.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, va_N, ha_N, cp, fld.im());
	mlout_fldtoim(dat, cp);
	name[9] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_genmovie(va_chr, ha_chr, cp, name, desc, dat, ntfr, dt, st.lambda, famp); 
	fclose(dat);
	return;
}

/* export full field data ("everything") into a viewer MATLAB file */
void cView::viewer() const
{
	FILE *dat;
	char name[13] = "__fld__A.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[5] = pid0;
	name[6] = pid1;
	fprintf(stderr, 
	   "viewer_%c%c([%g (%d) %g] x [%g (%d) %g] @ %g) >> %s\n", 
	   va_chr, ha_chr, va_beg, va_N, va_end, ha_beg, ha_N, ha_end, pos, name);
	dat = fopen(name, "w+");
	name[8] = 0; 
	mlout_viewertop(dat, name, va_chr, ha_chr, TE, st.lambda, 1);
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);

	mlout_fld(dat, va_N, ha_N, EX, ex.re()); mlout_fldtore(dat, EX);
	mlout_fld(dat, va_N, ha_N, EX, ex.im()); mlout_fldtoim(dat, EX);
	mlout_fld(dat, va_N, ha_N, EY, ey.re()); mlout_fldtore(dat, EY);
	mlout_fld(dat, va_N, ha_N, EY, ey.im()); mlout_fldtoim(dat, EY);
	mlout_fld(dat, va_N, ha_N, EZ, ez.re()); mlout_fldtore(dat, EZ);
	mlout_fld(dat, va_N, ha_N, EZ, ez.im()); mlout_fldtoim(dat, EZ);

	mlout_fld(dat, va_N, ha_N, HX, hx.re()); mlout_fldtore(dat, HX);
	mlout_fld(dat, va_N, ha_N, HX, hx.im()); mlout_fldtoim(dat, HX);
	mlout_fld(dat, va_N, ha_N, HY, hy.re()); mlout_fldtore(dat, HY);
	mlout_fld(dat, va_N, ha_N, HY, hy.im()); mlout_fldtoim(dat, HY);
	mlout_fld(dat, va_N, ha_N, HZ, hz.re()); mlout_fldtore(dat, HZ);
	mlout_fld(dat, va_N, ha_N, HZ, hz.im()); mlout_fldtoim(dat, HZ);

	Dmatrix n(va_N, ha_N);
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) { n(i, j) = st.n(vp, hp); vp += dv; }
		hp += dh;
	}
	mlout_n(dat, va_N, ha_N, n);
	mlout_fldviewer(dat, name, va_chr, ha_chr);
	fclose(dat);
	return;
}


/* ------------------------------------------------------------------- 
 * container for field levels at specific points 
 */

/* initialize */
cProbe::cProbe()
{
	x = 0.0;
	y = 0.0;
	z = 0.0;
	id = '?';
	eps = 1.0;
	ex = CC0;
	ey = CC0;
	ez = CC0;
	hx = CC0;
	hy = CC0;
	hz = CC0;
}
cProbe::cProbe(double cx, double cy, double cz, char ic)
{
	x = cx;
	y = cy;
	z = cz;
	id = ic;
	eps = 1.0;
	ex = CC0;
	ey = CC0;
	ez = CC0;
	hx = CC0;
	hy = CC0;
	hz = CC0;
}

// field levels
// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
// foa: one of MOD, ORG, SQR, REP, IMP
// id: the field "strength", one of mE, mH, qE, qH, mW 
Complex cProbe::field(Fcomp cp) const
{
	Complex f;
	switch(cp)
	{
	case EX: return ex; break;
	case EY: return ey; break;
	case EZ: return ez; break;
	case HX: return hx; break;
	case HY: return hy; break;
	case HZ: return hz; break;
	case SX: 
		f = ey*hz.conj()-ez*hy.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	case SY: 
		f = ez*hx.conj()-ex*hz.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	case SZ: 
		f = ex*hy.conj()-ey*hx.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	default: 
		cbundleerror("cProbe::field, cp?");
		break;
	
	}
	return CC0;
}
double cProbe::field(Fcomp cp, Afo foa) const
{
	return realize(field(cp), foa);
}
double cProbe::field(FSid id) const
{
	switch(id)
	{
	case mE:
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
		break;
	case mH:
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
		break;
	case qE:
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
		break;
	case qH:
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
		break;
	case mW:
		return  0.25*(  MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			+EP0*eps*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default:
		cbundleerror("cProbe::field, FSid?");
                break;
	}
	return 0.0;
}
   
/* ------------------------------------------------------------------- 
 * an array of probes
 */
   
/* initialize */
cProbeArray::cProbeArray()
          : num(0)
{
	pvec = NULL;
}

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

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

/* delete all cProbe entries */
void cProbeArray::clear()
{
	if(pvec != NULL) delete[] pvec;
	pvec = NULL;
	num = 0; 
}

/* subscripting */
cProbe& cProbeArray::operator() (int i)  
{ 
	if(i<0 || i>=num) cbundleerror("cProbeArray: () out of range");
	return pvec[i]; 
}
cProbe cProbeArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) cbundleerror("cProbeArray: () out of range");
	return pvec[i]; 
}

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



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

/* 
 * conventional coupled mode theory, 
 * codirectional guided wave propagation along parallel waveguide cores
 */

/* minimum distance between boundaries to be considered separated  */
#define HDIST 1.0e-8

/* subtract the permittivity profile in waveguide wp from that in wu  
   the PERMITTIVITY difference is returned in the Dmatrix n of wd */
#define MAXNUMBD 100
Waveguide permdifference(Waveguide wu, Waveguide wp)
{
	Waveguide wd;

	double x;
	double p, u;
	int i, ip, iu, id;
	Dvector th(MAXNUMBD);
	Interval r;
	int l;

	iu = 0; ip = 0; id = 0;
	while((iu < wu.nx+1 || ip < wp.nx+1) && id < MAXNUMBD-1) 
	{
		u = wu.layer(iu).x1;
		p = wp.layer(ip).x1;
		if(fabs(u-p) < HDIST) { x = u; ++ip; ++iu; }
		else
		{
			if(u < p) { x = u; ++iu; }
			else      { x = p; ++ip; }
		}
		th(id) = x; 
		++id;
	}
	if(id>=MAXNUMBD-1) cmterror("MAXNUMBD exceeded");
	--id; 
	if(id<=-1) cmterror("empty structures (x)");
	wd.nx = id; 
	wd.hx = Dvector(id+1);
	for(i=0; i<=id; ++i) wd.hx(i) = th(i);

	wd.n = Dvector(wd.nx+2);
	for(l=0; l<=wd.nx+1; ++l)
	{
			r = wd.layer(l); 
			x = 0.5*(r.x0+r.x1);
			u = wu.eps(x);
			p = wp.eps(x);
			wd.n(l) = u-p; 	
	}

	wd.lambda = wu.lambda;
	wd.special = 1;
	return wd;
}

/* setup coupling matrices */
void CMTmats(ModeArray& umode,     // modes to be coupled
	     Waveguide cpwg,       // the entire coupler structure
	     Dmatrix& cmtS,        // output: power coupling matrix
	     Dmatrix& cmtBK)       // output: coupling coefficients 
{
	int n;
	n = umode.num;
	Waveguide wd;
	Dmatrix s(n, n);
	Dmatrix k(n, n);
	Dmatrix tmp(n, n);
	Dvector p(n);
	Dvector b(n);
	Dmatrix a(n, n);
	int i, j, l;
	
	double wl;
	Polarization pl;
	Mode modei;
	Mode modej;

	Interval r;

	if(n<=1) cmterror("Grrr: n");
	modei = umode(0);
	pl = modei.pol;
	wl = modei.wg.lambda;
	if(cpwg.bdmatch(modei.wg) != 0) cmterror("wg.hx, wg.hy");
	p(0) = modei.power();
	b(0) = modei.beta;
	for(i=1; i<=n-1; ++i)
	{
		modei = umode(i);
		if(modei.pol != pl) cmterror("pol");
		if(modei.wg.lambda != wl) cmterror("lambda");
		if(cpwg.bdmatch(modei.wg) != 0) cmterror("wg.hx");
		p(i) = modei.power(); 
		b(i) = modei.beta; 
	}

	r = Interval(-AWAY, AWAY);
	for(i=0; i<=n-1; ++i)
	{
		modei = umode(i);
		for(j=0; j<=n-1; ++j)
		{
			if(i==j) tmp(i,j) = p(i);
			else     tmp(i,j) = 0.5*
				(modei.integrate(EX, umode(j), HY, r)
				-modei.integrate(EY, umode(j), HX, r));	
		}
	}
	for(i=0; i<=n-1; ++i)
	{
		for(j=0; j<=n-1; ++j)
		{
			s(i,j) = 0.5*(tmp(i,j)+tmp(j,i))/sqrt(p(i)*p(j));
		}
	}
			
	for(i=0; i<=n-1; ++i)
	{
		modei = umode(i);
		wd = permdifference(cpwg, modei.wg);
		for(j=0; j<=n-1; ++j) 
		{
			modej = umode(j);
			tmp(i,j) = 0.0;
			for(l=0; l<=wd.nx+1; ++l)
			{
				if(fabs(wd.n(l)) > 1.0e-10)
				{
					r = wd.layer(l);
					tmp(i,j) += wd.n(l)*(
				modei.integrate(EX, modej, EX, r)
			      + modei.integrate(EY, modej, EY, r)
			      + modei.integrate(EZ, modej, EZ, r));
				}
			}
		}
	}
	for(i=0; i<=n-1; ++i)
	{
		for(j=0; j<=n-1; ++j) 
		{
			k(i,j) = 1.0/val_invomep0(wl)
			         /4.0/sqrt(p(i)*p(j))
				 *0.5*(tmp(i,j)+tmp(j,i));
			k(i,j) += s(i,j)
			          *0.5*(b(i) + b(j));
		}
	}

	cmtS  = s;
	cmtBK = k; 
	return;
}

/* Coupled mode theory:
   setup coupling matrices, calculate supermode propagation constants
   and normalized amplitude vectors */
void CMTsetup(ModeArray& umode,     // modes to be coupled
	      Waveguide cpwg,       // the entire coupler structure
	      Dmatrix& sigma,       // output: power coupling matrix
	      Dvector& smpc,        // output: supermode propagation constants
	      Dmatrix& smav)        // output: supermode amplitude vectors 
{
	int n;
	n = umode.num;
	Dmatrix s(n, n);
	Dmatrix k(n, n);
	Dmatrix a(n, n);
	Dvector b(n);

	fprintf(stderr, "* CMT setup (%d) - ", n);

	CMTmats(umode, cpwg, s, k);

	solvegeneralEVproblem(k, s, b, a);

	sigma = s;	
	smpc = b; 
	smav = a;

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

	return; 
}


/* coupled mode transfer matrix */
Cmatrix CMTtrafm(double len,     // length of the coupling segment
                                 // as output from CMTsetup:
                 Dmatrix& sigma, // power coupling matrix
                 Dvector& smpc,  // supermode propagation constants
                 Dmatrix& smav)  // supermode amplitude vectors
{
	int n;
	n = sigma.r;

	Dvector pow(n);
        for(int k=0; k<=n-1; ++k)
		pow(k) = scalp(smav.col(k), mult(sigma, smav.col(k)));

	Cmatrix t(n, n);
	Complex efac;
	Dmatrix a, at;
	t.init(CC0);
	for(int s=0; s<=n-1; ++s) 
	{
		a = smav.submatrix(n, 1, 0, s);
		at = a;
		at.transpose();
		efac = expi(-smpc(s)*len);
		t = add(t, mult(mult(a, at), efac/pow(s)));
	}
	t = mult(t, sigma);
	return t;
}