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

/* 
 * structure.cpp
 * Multilayer slab waveguides, waveguide sequences, 
 * rectangular optical circuits. 
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"matlvis.h"

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

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

/* a normalized parameter that specifies whether localized fields related 
   to pieces of a waveguide are used to establish approximate modal 
   basis sets */
double DECOUPLED = 1.0e-6;

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

/* interval from the waveguide cross section */
/* initialize */
Interval::Interval()
{
	x0 = 0.0;
	x1 = 1.0;
}
Interval::Interval(double xa, double xb)
{
	if(xa <= xb) { x0 = xa; x1 = xb; }
	else         { x0 = xb; x1 = xa; }
}

/* output to FILE dat */
void Interval::write(FILE *dat) const
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "[%g, %g]\n", x0, x1);
	}
	else     
	{
		fputdouble(x0, dat);
		fputdouble(x1, dat);
	}
	return;
}

/* input from FILE dat */
void Interval::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		structerror("Interval: read");
	}
	else 
	{
		x0 = fgetdouble(dat);
		x1 = fgetdouble(dat);
	}
	return;
}

/* compare the present object to i, return 1 if equal, otherwise return 0 */
int Interval::equal(const Interval& i) const
{
	if(this == &i) return 1;
	if(fabs(x0-i.x0) > COMPTOL_IV) return 0; 
	if(fabs(x1-i.x1) > COMPTOL_IV) return 0; 
	return 1;
}

/* restrict iv to the part inside the present object,
   return value: 0: no overlap, 1: something from iv remains */
int Interval::clip(Interval& iv) const
{
	if(x0 >= x1) return 0;
	if(iv.x0 >= iv.x1) return 0;
	if(iv.x1 <= x0) return 0;
	if(iv.x0 >= x1) return 0;
	if(x0 > iv.x0) iv.x0 = x0;
	if(x1 < iv.x1) iv.x1 = x1;
	return 1;
}


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

/* waveguide geometry + permittivity profile; vacuum wavelength */

/* initialize */
Waveguide::Waveguide()
{
	nx = 0;
	hx = Dvector();
	lambda = 1.0;
	n = Dvector();
	special = 0;
}
Waveguide::Waveguide(int vnx) 
{
	nx = vnx;
	hx = Dvector(nx+1);
	lambda = 1.0;
	n = Dvector(nx+2);
	special = 0;
}
Waveguide::Waveguide(int vnx, Dvector vhx, double l, Dvector vn)
{
	nx = vnx;
	hx = vhx;
	lambda = l;
	n = vn;
	special = 0;
}

/* free allocated memory */
void Waveguide::strip()
{
	nx = 0;
	hx.strip();
	n.strip();
}

/* get Interval index corresponding to position x */ 
int Waveguide::layeridx(double x) const
{
	int l=0;
	while(l<=nx && x>hx(l)) ++l;
	return l;
}

/* get layer boundaries corresponding to index l */  
Interval Waveguide::layer(int l) const
{
	Interval i;
	if(l <= 0) i.x0 = -AWAY;
	else i.x0 = hx(l-1);
	if(l >= nx+1) i.x1 = AWAY;
	else i.x1 = hx(l);
	return i;
}

/* get layer boundaries corresponding to position x */ 
Interval Waveguide::layer(double x) const
{
	return layer(layeridx(x));
}

/* permittivity on layer l */
double Waveguide::eps(int l) const
{
	double e;
	e = n(l);
	if(special) return e;
	else        return e*e;
}

/* permittivity at position x */
double Waveguide::eps(double x) const
{
	double e;
	int l;
	l = layeridx(x);
	e = n(l);
	if(special) return e;
	else        return e*e;
}

/* upper bound for effective mode permittivities, default value */
double Waveguide::defaultepseffmax() const
{
	int l;
	double eeffmax = -MAXREFIND*MAXREFIND;
	for(l=1; l<=nx; ++l) { if(eeffmax < eps(l)) eeffmax = eps(l); }
	return eeffmax;
}

/* lower bound for effective mode permittivities, default value */
double Waveguide::defaultepseffmin() const
{
	double eeffmin;
	eeffmin = eps(0);
	if(eeffmin<eps(nx+1)) eeffmin = eps(nx+1);
	return eeffmin;
}

/* translate: hx -> hx+dx */
void Waveguide::translate(double dx)
{
	int l;
	for(l=0; l<=nx; ++l) hx(l) += dx;
	return;
}

/* output to FILE dat */
void Waveguide::write(FILE *dat) const
{
	int i;
	if(dat == stderr || dat == stdout)
	{
		comment("Waveguide", dat);
		fprintf(dat, "Nx: %d, Hx: ", nx);
		for(i=0; i<=nx; ++i) fprintf(dat, "%g ", hx(i));
		fprintf(dat, "\n");
		fprintf(dat, "lambda: %g\n", lambda);
		fprintf(dat, "n:\n");
		for(i=0; i<=nx+1; ++i) fprintf(dat, "%g  ", n(i));
		fprintf(dat, "\n");
		fprintf(dat, "eps:\n");
		for(i=0; i<=nx+1; ++i) fprintf(dat, "%g  ", eps(i));
		fprintf(dat, "\n");
		fprintf(dat, "special: %d\n", special);
	}
	else
	{
		comment("Waveguide", dat);
		comment("nx", dat);
		fputint(nx, dat);
		comment("hx[0..nx]", dat);
		hx.write(dat);
		comment("lambda", dat);
		fputdouble(lambda, dat);
		comment("n[0..nx+1]", dat);
		n.write(dat);
		comment("special", dat);
		fputint(special, dat);
	}
	return;
}

/* input from FILE dat */
void Waveguide::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		structerror("Waveguide: read");
	}
	nx = fgetint(dat);
	hx.read(dat);
	lambda = fgetdouble(dat);
	n.read(dat);
	special = fgetint(dat);
	return;
}

/* test matching of boundary positions with another waveguide,
   all boundaries of this object should appear in w */
int Waveguide::bdmatch(Waveguide w) const
{
        int j, k, e;
        for(j=0; j<=w.nx; ++j)
        {
                e = 0;
                for(k=0; k<=nx; ++k)
                {
                        if(fabs(hx(k)-w.hx(j))<=HDIST) e = 1;
                }
                if(e != 1) return 1;
        }
	return 0;
}      

/* compare the present object to wg, return 1 if equal, otherwise return 0 */
int Waveguide::equal(const Waveguide& wg) const 
{
	if(this == &wg) return 1;
	if(nx != wg.nx) return 0;
	if(fabs(lambda-wg.lambda) > COMPTOL_LAMBDA) return 0; 
	for(int l=0; l<=nx+1; ++l)
		if(fabs(n(l)-wg.n(l)) > COMPTOL_N) return 0; 
	for(int l=0; l<=nx; ++l)
		if(fabs(hx(l)-wg.hx(l)) > COMPTOL_HX) return 0; 
	if(special != wg.special) return 0;
	return 1;
}

/* test for constant refractive index, return 1, if constant, 0 otherwise */
int Waveguide::constmed() const 
{
	double tn = n(0);
        for(int l=1; l<=nx+1; ++l)
        if(fabs(tn - n(l))> COMPTOL_N) return 0;
	return 1;
}

/* remove unnecessary boundaries from waveguide, 
   down to a minimum of minnx inner layers */
void Waveguide::optimize(int minnx)
{
	int done = 0;
	int rl = 0;
	while(nx >= (minnx+1) && done != 1)	
	{	
		done = 1;
		for(int l=0; l<=nx; ++l)
		{
			if(done == 1 && fabs(n(l)-n(l+1)) < COMPTOL_N)
			{
				done = 0;
				rl = l; 
			}
		}
		if(done != 1)
		{
			for(int l=rl; l<=nx; ++l) n(l) = n(l+1);
			for(int l=rl; l<=nx-1; ++l) hx(l) = hx(l+1);
			--nx;
		}
	}
	Dvector n_n = n.subvector(nx+2, 0);
	Dvector n_hx = hx.subvector(nx+1, 0);
	n = n_n;
	hx = n_hx;
	return;
}
/* ... down to a minimum of one inner layer */
void Waveguide::optimize()
{
	optimize(1);
	return;
}

/* check consistency of the waveguide object: error&exit, if not ok */
void Waveguide::consistency() const
{
	if(nx <= 0 || nx >= 1000) structerror("Waveguide: consistency: nx");
	if(hx.nel != nx+1) structerror("Waveguide: consistency, hx: #elements");
	for(int b=0; b<=nx-1; ++b)	
	{
		if(hx(b) > hx(b+1) - HDIST)
		structerror("Waveguide: consistency, hx: order/distance");
	}	
	if(lambda < 0.001 && lambda > 100.0)
		structerror("Waveguide: consistency, lambda: magnitude");
	if(n.nel != nx+2) structerror("Waveguide: consistency, n: #elements");
	if(special)
	{
		for(int l=0; l<=nx+1; ++l)	
		{
			if(   eps(l)<-MAXREFIND*MAXREFIND 
			   || eps(l)> MAXREFIND*MAXREFIND)
		structerror("Waveguide: consistency(special), eps: magnitude");
		}
	}
	else
	{
		for(int l=0; l<=nx+1; ++l)	
		{
			if(n(l)<MINREFIND || n(l)>MAXREFIND)
			structerror("Waveguide: consistency, n: magnitude");
		}
	}
	return;
}

/* check the symmetry of the waveguide, return value l
   l == 0: no symmetry
   l >= 1: symmetric structure with respect to layer l */
int Waveguide::checksymmetry() const
{
	if(nx % 2 == 0) return 0;
	double c = (hx(nx)+hx(0))/2.0;
	int cl = layeridx(c);
	for(int l=0; l<=nx; ++l)
	{
		if(fabs((c-hx(l)) - (hx(nx-l)-c))> COMPTOL_HX) return 0;
	}
	for(int l=0; l<=nx+1; ++l)
	{
		if(fabs(n(l)-n(nx-l+1)) > COMPTOL_N) return 0;
	}
	return cl;
}

/* split the present waveguide into two, at layer l
   the pieces are stored in wg1 and wg2 */
void Waveguide::split(int sl, Waveguide& wg1, Waveguide& wg2) const
{
	if(sl == 1)
	{
		wg1 = Waveguide(1);
		wg1.hx(0) = hx(0)-LIMBDSEP/2.0;
		wg1.hx(1) = hx(0);
		wg1.n(0) = n(0);
		wg1.n(1) = n(0);
		wg1.n(2) = n(1);
		wg1.lambda = lambda;
		wg1.special = special;
	}
	else
	{
		int nnx = sl-1;
		wg1 = Waveguide(nnx);
		for(int l=0; l<=nnx; ++l) wg1.hx(l) = hx(l);
		for(int l=0; l<=nnx+1; ++l) wg1.n(l) = n(l);
		wg1.lambda = lambda;
		wg1.special = special;
	}
	if(sl == nx)
	{
		wg2 = Waveguide(1);
		wg2.hx(0) = hx(nx);
		wg2.hx(1) = hx(nx)+LIMBDSEP/2.0;
		wg2.n(0) = n(nx);
		wg2.n(1) = n(nx+1);
		wg2.n(2) = n(nx+1);
		wg2.lambda = lambda;
		wg2.special = special;
	}
	else
	{
		int nnx = nx-sl;
		wg2 = Waveguide(nnx);
		for(int l=0; l<=nnx; ++l) wg2.hx(l) = hx(sl+l);
		for(int l=0; l<=nnx+1; ++l) wg2.n(l) = n(sl+l);
		wg2.lambda = lambda;
		wg2.special = special;
	}
	return;
}

/* expand the waveguide towards a structure 
   that is symmetric with respect to xs; xs > hx(nx) */
Waveguide Waveguide::expand(double xs)
{
	if(xs < hx(nx)+HDIST) structerror("Waveguide: expand: xs");
	int nnx = 2*nx+1;
	Waveguide wg(nnx);
	for(int b=0; b<=nx; ++b) 
	{
		wg.hx(b) = hx(b);
		wg.hx(nnx-b) = xs+(xs-hx(b));
	}
	for(int l=0; l<=nx; ++l)
	{
		 wg.n(l) = n(l);
		 wg.n(nnx+1-l) = n(l);
	}
	wg.n(nx+1) = n(nx+1);
	wg.lambda = lambda;
	wg.special = special;
	return wg;
}

/* mirror the waveguide at xs; xs > hx(nx) */
Waveguide Waveguide::mirror(double xs)
{
	if(xs < hx(nx)+HDIST) structerror("Waveguide: mirror: xs");
	Waveguide wg(nx);
	for(int b=0; b<=nx; ++b) 
	{
		wg.hx(nx-b) = xs+(xs-hx(b));
	}
	for(int l=0; l<=nx+1; ++l)
	{
		 wg.n(nx+1-l) = n(l);
	}
	wg.lambda = lambda;
	wg.special = special;
	return wg;
}

/* determine a limit for the effective permittivity, beyond which the inner 
   layer l optically decouples the neighbouring structures */
double Waveguide::decoupepseff(int l)
{
	if(l<=0 || l>=nx+1) structerror("Waveguide: decoupepseff: l");
	double d = hx(l)-hx(l-1);
	double lnD = log(DECOUPLED);
	lnD = lnD/2/PI;
	return eps(l)+lnD*lnD*lambda*lambda/d/d;
}

/* determine a limit for the effective permittivity, beyond which the lower 
   boundary at position p becomes irrelevant */
double Waveguide::lbdecepseff(double p)
{
	if(p >= hx(0)-HDIST) structerror("Waveguide: lbdecepseff: p");
	double d = hx(0)-p;
	d *= 2.0;
	double lnD = log(DECOUPLED);
	lnD = lnD/2/PI;
	return eps(0)+lnD*lnD*lambda*lambda/d/d;
}

/* determine a limit for the effective permittivity, beyond which the upper 
   boundary at position p becomes irrelevant */
double Waveguide::ubdecepseff(double p)
{
	if(p <= hx(nx)+HDIST) structerror("Waveguide: ubdecepseff: p");
	double d = p-hx(nx);
	d *= 2.0;
	double lnD = log(DECOUPLED);
	lnD = lnD/2/PI;
	return eps(nx+1)+lnD*lnD*lambda*lambda/d/d;
}

/* refractive index profile -> MATLAB .m file  
   disp: output region on the x axis
   ext0, ext1: filename id characters */
void Waveguide::plot(Interval disp, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "n__.m";
	name[1] = ext0;
	name[2] = ext1;
	fprintf(stderr, "n [%g, %g] >> %s\n", disp.x0, disp.x1, name);
	if(special) fprintf(stderr, "  (special: eps<->n !)\n");
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Waveguide profile");
	name[3] = 0;
	double minn, maxn;
	minn = n.min();
	maxn = n.max();
	if(special) 
	{
		minn = sgn_sqrt(minn);
		maxn = sgn_sqrt(maxn);
	}
	mlout_geo(dat, (*this), minn, maxn); 
	mlout_refindplot(name, dat, disp); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}
void Waveguide::plot(char ext0, char ext1) const
{
	Interval disp;
	double w = (hx(nx)-hx(0))/((double) nx);
	disp.x0 = hx(0)-w; 
	disp.x1 = hx(nx)+w; 
	plot(disp, ext0, ext1); 
	return;
}

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

/* segmented planar waveguide: geometry + permittivity profile;
                               vacuum wavelength ... */         

/* initialize */
SegWgStruct::SegWgStruct()
{
	nz = -1;
	wgvec = new Waveguide[nz+2];
	Waveguide wg(1);
	wg.n(0) =  1.0;
	wg.n(1) =  1.0;
	wg.n(2) =  1.0;
	wg.hx(0) = -1.0;
	wg.hx(1) =  1.0;
	wg.lambda = 1.0;
	wg.special = 0;
	for(int i=0; i<=nz+1; ++i) wgvec[i] = wg;
	hz = Dvector(nz+1);
	for(int i=0; i<=nz; ++i) hz(i) = ((double) i);
}
SegWgStruct::SegWgStruct(int vnz)
{
	if(vnz<-1) structerror("SegWgStruct: init < -1"); 
	nz = vnz;
	wgvec = new Waveguide[nz+2];
	Waveguide wg(1);
	wg.n(0) =  1.0;
	wg.n(1) =  1.0;
	wg.n(2) =  1.0;
	wg.hx(0) = -1.0;
	wg.hx(1) =  1.0;
	wg.lambda = 1.0;
	wg.special = 0;
	for(int i=0; i<=nz+1; ++i) wgvec[i] = wg;
	hz = Dvector(nz+1);
	for(int i=0; i<=nz; ++i) hz(i) = ((double) i);
}
/* convert a Circuit to a SegWgStruct 
   opt = 1: merge equal segments & layers (default) */
SegWgStruct::SegWgStruct(Circuit rc, int opt)
{
	nz = rc.nz;
	hz = rc.hz;
	wgvec = new Waveguide[nz+2];
	for(int s=0; s<=nz+1; ++s)
	{
		Waveguide wg(rc.nx);
		wg.hx = rc.hx;
		wg.n = rc.n.col(s);
		wg.lambda = rc.lambda;
		wg.special = rc.special;
		if(opt == 1) wg.optimize();
		wgvec[s] = wg;
	}
	if(opt == 1) optimize();
}
SegWgStruct::SegWgStruct(Circuit rc)
{
	wgvec = NULL;
	(*this) = SegWgStruct(rc, 1);
}

/* destroy */
SegWgStruct::~SegWgStruct()
{
	if(wgvec != NULL) delete[] wgvec;
	wgvec = NULL;
	nz = 0; 
	hz.strip();
}

/* copy constructor */
SegWgStruct::SegWgStruct(const SegWgStruct& swg)
{
	nz = swg.nz;
	hz = swg.hz;
	wgvec = new Waveguide[nz+2];
	for(int i=0; i<=nz+1; ++i) wgvec[i] = swg.wgvec[i]; 
}
	
/* assignment */
SegWgStruct& SegWgStruct::operator=(const SegWgStruct& swg)
{
	if(this != &swg)
	{
		if(wgvec != NULL) delete[] wgvec;
		wgvec = NULL;
		nz = swg.nz; 
		hz = swg.hz; 
		wgvec = new Waveguide[nz+2];
		for(int i=0; i<=nz+1; ++i) wgvec[i] = swg.wgvec[i]; 
	}
	return *this;
}

/* delete all segment entries */
void SegWgStruct::clear()
{
	nz = -1; 
	hz.strip();
	if(wgvec != NULL) delete[] wgvec; 
	wgvec = NULL;
}

/* set all segments to wg */
void SegWgStruct::init(Waveguide wg)
{
	for(int i=0; i<=nz+1; ++i) wgvec[i] = wg; 
}

/* subscripting 
Waveguide& SegWgStruct::operator() (int s)
{ 	
	if(s < 0) structerror("SegWgStruct: access ( < 0)"); 
	if(s > nz+1) structerror("SegWgStruct: access ( > nz+1)");
	return wgvec[s]; 
}
Waveguide SegWgStruct::operator() (int s) const 
{
	if(s < 0) structerror("SegWgStruct: access ( < 0)"); 
	if(s > nz+1) structerror("SegWgStruct: access ( > nz+1)");
	return wgvec[s]; 
}
*/

/* input from FILE dat */
void SegWgStruct::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		structerror("SegWgStruct: read");
	}
	else
	{
		if(wgvec != NULL) delete[] wgvec;
		wgvec = NULL;
		nz = fgetint(dat);
		if(nz<-1) structerror("SegWgStruct: read: nz < -1");
		hz.read(dat);
		wgvec = new Waveguide[nz+2];
		for(int i=0; i<=nz+1; ++i) (wgvec[i]).read(dat);
	}
	return;
}

/* output to FILE dat */
void SegWgStruct::write(FILE *dat) const
{
	comment("SegWgStruct", dat);
	comment("nz", dat); fputint(nz, dat);
	comment("hz", dat); hz.write(dat); 
	for(int i=0; i<=nz+1; ++i) (wgvec[i]).write(dat);
	return;
}

/* vacuum wavelength: should be unique in all segments */
double SegWgStruct::lambda()
{
	return wgvec[0].lambda;
}


/* check consistency of the SegWgStruct object: error&exit, if not ok */
void SegWgStruct::consistency() const
{
	double l;
	l = wgvec[0].lambda;
	for(int i=1; i<=nz+1; ++i)
	{
		if(wgvec[i].lambda !=l) 
			structerror("SegWgStruct: lambda not unique");
	}
	int sp; 
	sp = wgvec[0].special;
	for(int i=1; i<=nz+1; ++i)
	{
		if(wgvec[i].special !=sp) 
			structerror("SegWgStruct: wg.special not unique");
	}
	for(int i=0; i<=nz+1; ++i) wgvec[i].consistency(); 
	return;
}

/* get segment index corresponding to position z */ 
int SegWgStruct::segidx(double z) const
{
	int s=0;
	while(s<=nz && z>hz(s)) ++s;
	return s;
}

/* get segment boundaries corresponding to index s */  
Interval SegWgStruct::segment(int s) const
{
	Interval i;
//	if(s < 0) structerror("SegWgStruct: segment( < 0)"); 
//	if(s > nz+1) structerror("SegWgStruct: segment( > nz+1)");
	if(s <= 0) i.x0 = -AWAY;
	else i.x0 = hz(s-1);
	if(s >= nz+1) i.x1 = AWAY;
	else i.x1 = hz(s);
	return i;
}

/* get segment boundaries corresponding to position z */ 
Interval SegWgStruct::segment(double z) const
{
	return segment(segidx(z));
}

/* get refractive index at position x, z */
double SegWgStruct::n(double x, double z) const
{
	int s = segidx(z);
	Waveguide wg = wgvec[s];
	double nv = wg.n(wg.layeridx(x));
	if(wg.special) return sgn_sqrt(nv);
	return nv;
}

/* get permittivity at position x, z */
double SegWgStruct::eps(double x, double z) const
{
	int s = segidx(z);
	Waveguide wg = wgvec[s];
	double nv = wg.n(wg.layeridx(x));
	if(wg.special) return nv;
	return nv*nv;
}

/* translate: hz -> hz+dz */
void SegWgStruct::ztranslate(double dz)
{
	for(int s=0; s<=nz; ++s) hz(s) += dz;
	return;
}

/* translate: hx -> hx+dx for all segments */
void SegWgStruct::xtranslate(double dx)
{
	for(int s=0; s<=nz+1; ++s) (wgvec[s]).translate(dx); 
	return;
}

/* combine equal adjacent segments */
void SegWgStruct::optimize()
{
	int done = 0;
	int rs = 0;
	while(nz >= 2 && done != 1)	
	{	
		done = 1;
		for(int s=0; s<=nz; ++s)
		{
			if(done == 1 && (wgvec[s]).equal((wgvec[s+1])) != 0)
			{
				done = 0;
				rs = s; 
			}
		}
		if(done != 1)
		{
			for(int s=rs; s<=nz; ++s) wgvec[s] = wgvec[s+1];
			for(int s=rs; s<=nz-1; ++s) hz(s) = hz(s+1);
			--nz;
		}
	}
	Waveguide *n_wgvec;
	n_wgvec = new Waveguide[nz+2];	
	for(int s=0; s<=nz+1; ++s) n_wgvec[s] = wgvec[s];
	delete[] wgvec;
        wgvec = n_wgvec;
	Dvector n_hz = hz.subvector(nz+1, 0);
	hz = n_hz;
	return;
}

/* QUEP: construct the corresponding vertical segmented waveguide sequence */
SegWgStruct SegWgStruct::rotate() const
{
	consistency();
	int spl = wgvec[0].special;
	int b=0;
	for(int s=0; s<=nz+1; ++s) b += wgvec[s].nx+1; 
	Dvector tbd(b);
	Dvector bd;
	b=0;
	for(int s=0; s<=nz+1; ++s)
	{
		tbd.setsubvector(wgvec[s].hx, b);		
		b += wgvec[s].nx+1;
	}
	tbd.sort();
	bd = tbd;
	if(bd.nel<2) structerror("SegWgStruct: rotated: (1) #x_bounds < 2.");
	tbd(0) = bd(0); b=1;
	for(int i=1; i<=bd.nel-1; ++i)
	{
		if(fabs(tbd(b-1)-bd(i)) > HDIST)
		{
			tbd(b) = bd(i);
			++b;
		}
	}
	bd = tbd.subvector(b, 0);
	if(bd.nel<2) structerror("SegWgStruct: rotated: (2) #x_bounds < 2.");

	int nnz = bd.nel-1;	

	SegWgStruct nst(nnz);
	nst.hz = bd;
	Waveguide twg(nz);
	twg.lambda = wgvec[0].lambda;
	twg.hx = hz;
	twg.n.init(1.0);
	nst.init(twg);
	
	Interval itv;
	double x, z, r;
	for(int s=0; s<=nst.nz+1; ++s)
	{	
		itv = nst.segment(s);
		x = 0.5*(itv.x0+itv.x1);
		for(int l=0; l<=twg.nx+1; ++l)
		{
			itv = twg.layer(l);
	                z = 0.5*(itv.x0+itv.x1);
			int idx = wgvec[segidx(z)].layeridx(x);
			r = wgvec[segidx(z)].n(idx);
			nst(s).n(l) = r;
		}
	}
	for(int s=0; s<=nst.nz+1; ++s) nst(s).optimize();
	nst.optimize();
	for(int s=0; s<=nst.nz+1; ++s) nst(s).special = spl;
	return nst;
}

/* QUEP: enlarge the sequence by identical input and output segments */
SegWgStruct SegWgStruct::expand(Interval w) const
{
	int nnz = nz+2;
	SegWgStruct nst(nnz);
	if(w.x0 >= hz(0)) 
		structerror("SegWgStruct: expand: illegal window (1)");
	if(w.x1 <= hz(nz)) 
		structerror("SegWgStruct: expand: illegal window (2)");
	nst.hz(0) = w.x0;
	for(int b=0; b<=nz; ++b) nst.hz(b+1) = hz(b);
	nst.hz(nnz) = w.x1;
	nst(0) = wgvec[0];
	for(int b=0; b<=nz+1; ++b) nst(b+1) = wgvec[b];
	nst(nnz+1) = wgvec[nz+1];
	return nst;
}

/* QUEP: remove the additional outer segments */
SegWgStruct SegWgStruct::shrink() const
{
	if(nz <= 2) structerror("SegWgStruct: shrink: nz<=2");
	int nnz = nz-2;
	SegWgStruct nst(nnz);
	for(int b=0; b<=nnz; ++b) nst.hz(b) = hz(b+1);
	for(int b=0; b<=nnz+1; ++b) nst(b) = wgvec[b+1];
	return nst;
}

/* determine a default window around the "inner" structure,
   e.g. for purposes of plotting */
void SegWgStruct::defaultwindow(Interval& wx, Interval& wz) const
{
	wx = wz = Interval(-1.0, 1.0);
	if(nz < -1) return;
	Waveguide wg = wgvec[0];
	double x0 = wg.hx(0);
	double x1 = wg.hx(wg.nx);
	double wl = wg.lambda;
	double emin = wg.defaultepseffmin();
	double xt;
	for(int j=1; j<=nz+1; ++j)
	{
		wg = wgvec[j];
		xt = wg.hx(0);
		if(xt < x0) x0 = xt;
		xt = wg.hx(wg.nx);
		if(xt > x1) x1 = xt;
		if(wg.eps(0) < emin) emin = wg.eps(0);
		if(wg.eps(wg.nx+1) < emin) emin = wg.eps(wg.nx+1);
	}
	if(emin < 1.0) emin=1.0;
	double d = 1.5*wl/sqrt(emin);
	x0 = x0-d;
	x1 = x1+d;
	double z0 = -d;
	double z1 = d;
	if(nz >= 0) { z0 = hz(0)-d; z1 = hz(nz)+d; }
	wx = Interval(x0, x1);
	wz = Interval(z0, z1);
	return;
}

/* convert this to a Circuit object */
Circuit SegWgStruct::circuit() const
{
	SegWgStruct h;
	if(nz >= 1) h = (*this);
	else
	{
		h = SegWgStruct(1);
		if(nz == -1) 
		{
			Waveguide wg = wgvec[0];
			h(0) = wg;
			h(1) = wg;
			h(2) = wg;	
			double d;
			d = wg.hx(wg.nx) - wg.hx(0);
			h.hz(0) = -d/2.0;
			h.hz(1) =  d/2.0;
		}
		else 
		{
			if(nz == 0) 
			{
				Waveguide wg = wgvec[0];
				h(0) = wg;
				h(1) = wg;
				h(2) = wgvec[1];
				double d;
				d = wg.hx(wg.nx) - wg.hx(0);
				h.hz(0) = hz(0)-d;
				h.hz(1) = hz(0);
			}
			else structerror("SegWgStruct: circuit: nz");
		}
	}
	SegWgStruct v;
	v = h.rotate();
	Circuit circ(v.nz, h.nz);
	circ.lambda = h.wgvec[0].lambda; 
	circ.special = h.wgvec[0].special; 
	circ.hx = v.hz;
	circ.hz = h.hz;
	Waveguide wg;
	double x;
	for(int s=0; s<=h.nz+1; ++s)
	{
		wg = h.wgvec[s];
		for(int l=0; l<=v.nz+1; ++l)
		{
			if(l>=1 && l<=v.nz) x = 0.5*(v.hz(l-1)+v.hz(l)); 
			else 
			{
				if(l==0) x = v.hz(0)-1.0;
				else     x = v.hz(v.nz)+1.0;
			}
			circ.n(l, s) = wg.n(wg.layeridx(x)); 
		}
	}
	return circ;
}

/* modelling of periodically corrugated waveguides:
   extend this periodically over the display window zbeg, zend */
SegWgStruct SegWgStruct::unwrap(double zbeg, double zend) const
{
	if(zbeg >= zend) structerror("SegWgStruct: unwrap: zbeg, zend");
	double L = hz(nz)-hz(0);
	double z0 = hz(0);
	while(zbeg <= z0) z0 -= L;
	double zn = hz(nz);
	while(zend >= zn) zn += L;
	int nnz = (int) floor((zn-z0)/L+0.5)*nz;
	SegWgStruct st(nnz);	
	int j=0;
	st(j) = (*this)(0);
	st.hz(j) = z0;
	++j;
	int k = 1;
	while(j<=nnz)
	{
		st(j) = (*this)(k);
		st.hz(j) = st.hz(j-1)+(hz(k)-hz(k-1));
		++j;
		++k; if(k > nz) k = 1;
	}
	st(nnz+1) = (*this)(nz+1);
	st.optimize();
	return st;
}
/* modelling of periodically corrugated waveguides:
   remove outermost equal slices */
SegWgStruct SegWgStruct::reduce() const
{
	if(nz < 3) structerror("SegWgStruct: reduce: nz<3");
        if((*this)(0).equal((*this)(1)) != 1) structerror("SegWgStruct: reduce: (0)!=(1)");
        if((*this)(nz+1).equal((*this)(nz)) != 1) structerror("SegWgStruct: reduce: (nz)!=(nz+1)"); 
	SegWgStruct st(nz-2);	
	for(int l=1; l<=nz; ++l) st(l-1) = (*this)(l);
	for(int l=1; l<=nz-1; ++l) st.hz(l-1) = (*this).hz(l);
	return st;
}

/* refractive index profile -> MATLAB .m file  
   ext0, ext1: filename id characters 
   xbeg, xend, zbeg, zend: display window in the x-z-plane */
void SegWgStruct::plot(double xbeg, double xend, double zbeg, double zend,
                       char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "n__.m";
	name[1] = ext0;
	name[2] = ext1;
	fprintf(stderr, "n [%g, %g] x [%g, %g] >> %s\n",  
	               xbeg, xend, zbeg, zend, name);
	if(special()) fprintf(stderr, "  (special: eps<->n !)\n");
	dat = fopen(name, "w+");
	mlout_title(dat, name, "2D refractive index profile");
	name[3] = 0;
	mlout_gengeoxz(dat, (*this), xbeg, xend, zbeg, zend);
	mlout_windowxz(dat, xbeg, xend, zbeg, zend);
	mlout_gengeorefindpatches(dat, (*this));
	mlout_refindoverview(name, dat); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}
void SegWgStruct::plot(char ext0, char ext1) const
{
	Interval wx, wz;
	defaultwindow(wx, wz);
	plot(wx.x0, wx.x1, wz.x0, wz.x1, ext0, ext1); 
	return;
}

/* special() != 0: n stores permittivity values in place of refractive
   indices, consistency checks for ref.ind. and perm. are disabled;
   a misuse of the Circuit container for problems with negative permittivity */
int SegWgStruct::special() const
{
	int sp = wgvec[0].special;
	for(int i=1; i<=nz+1; ++i)
	{
		if(wgvec[i].special !=sp) 
			structerror("SegWgStruct: special() not unique");
	}
	return sp;
}
/* ------------------------------------------------------------------- */

/* rectangular 2D optical circuit,
   geometry & permittivity profile, vacuum wavelength;
   (so far) meant only as a help to define a SegWgStruct */

/* initialize */
Circuit::Circuit()
{
	(*this) = Circuit(1, 1);
}
Circuit::Circuit(int vnx, int vnz)
{
	if(vnx < 1) structerror("Circuit: vnx < 1");
	if(vnz < 1) structerror("Circuit: vnz < 1");
	nx = vnx;
	hx = Dvector(nx+1);
	nz = vnz;
	hz = Dvector(nz+1);
	n = Dmatrix(nx+2, nz+2);
	special = 0;
}

/* refractive index profile -> MATLAB .m file  
   ext0, ext1: filename id characters 
   xbeg, xend, zbeg, zend: display window in the x-z-plane */
void Circuit::plot(double xbeg, double xend, double zbeg, double zend,
                   char ext0, char ext1) const
{
	SegWgStruct s(*this);
	s.plot(xbeg, xend, zbeg, zend, ext0, ext1); 
	return;
}
void Circuit::plot(char ext0, char ext1) const
{
	SegWgStruct s(*this);
	s.plot(ext0, ext1); 
	return;
}

/* set Rectangle index (l, m) corresponding to position (x,z) */ 
void Circuit::rectidx(double x, double z, int& l, int& m) const
{
	l = 0;
	while(l<=nx && x>hx(l)) ++l;
	m = 0;
	while(m<=nz && z>hz(m)) ++m;
	return;
}

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

/* for the VEIMS solver:
   cross section (2D) of a 3D rectangular waveguide,
   defined as a sequence of (1D) Waveguides */

/* initialize */
SegWgCrs::SegWgCrs()
{
	ny = -1;
	wgvec = new Waveguide[ny+2];
	Waveguide wg(1);
	wg.n(0) =  1.0;
	wg.n(1) =  1.0;
	wg.n(2) =  1.0;
	wg.hx(0) = -1.0;
	wg.hx(1) =  1.0;
	wg.lambda = 1.0;
	wg.special = 0;
	for(int i=0; i<=ny+1; ++i) wgvec[i] = wg;
	hy = Dvector(ny+1);
	for(int i=0; i<=ny; ++i) hy(i) = ((double) i);
}
SegWgCrs::SegWgCrs(int vny)
{
	if(vny<-1) structerror("SeWgCrs: init < -1"); 
	ny = vny;
	wgvec = new Waveguide[ny+2];
	Waveguide wg(1);
	wg.n(0) =  1.0;
	wg.n(1) =  1.0;
	wg.n(2) =  1.0;
	wg.hx(0) = -1.0;
	wg.hx(1) =  1.0;
	wg.lambda = 1.0;
	wg.special = 0;
	for(int i=0; i<=ny+1; ++i) wgvec[i] = wg;
	hy = Dvector(ny+1);
	for(int i=0; i<=ny; ++i) hy(i) = ((double) i);
}
SegWgCrs::SegWgCrs(WgCrs rc)
{
	ny = rc.ny;
	hy = rc.hy;
	wgvec = new Waveguide[ny+2];
	for(int s=0; s<=ny+1; ++s)
	{
		Waveguide wg(rc.nx);
		wg.hx = rc.hx;
		wg.n = rc.n.col(s);
		wg.lambda = rc.lambda;
		wg.special = 0;
		wg.optimize();
		wgvec[s] = wg;
	}
}

/* destroy */
SegWgCrs::~SegWgCrs()
{
	if(wgvec != NULL) delete[] wgvec;
	wgvec = NULL;
	ny = 0; 
	hy.strip();
}

/* copy constructor */
SegWgCrs::SegWgCrs(const SegWgCrs& swg)
{
	ny = swg.ny;
	hy = swg.hy;
	wgvec = new Waveguide[ny+2];
	for(int i=0; i<=ny+1; ++i) wgvec[i] = swg.wgvec[i]; 
	
}
	
/* assignment */
SegWgCrs& SegWgCrs::operator=(const SegWgCrs& swg)
{
	if(this != &swg)
	{
		if(wgvec != NULL) delete[] wgvec;
		wgvec = NULL;
		ny = swg.ny; 
		hy = swg.hy; 
		wgvec = new Waveguide[ny+2];
		for(int i=0; i<=ny+1; ++i) wgvec[i] = swg.wgvec[i]; 
	}
	return *this;
}

/* delete all segment entries */
void SegWgCrs::clear()
{
	ny = -1; 
	hy.strip();
	if(wgvec != NULL) delete[] wgvec; 
	wgvec = NULL;
}

/* set all segments to wg */
void SegWgCrs::init(Waveguide wg)
{
	for(int i=0; i<=ny+1; ++i) wgvec[i] = wg; 
}

/* subscripting 
Waveguide& SegWgCrs::operator() (int s)
{ 	
	if(s < 0) structerror("SegWgCrs: access ( < 0)"); 
	if(s > ny+1) structerror("SegWgCrs: access ( > ny+1)");
	return wgvec[s]; 
}
Waveguide SegWgCrs::operator() (int s) const 
{
	if(s < 0) structerror("SegWgCrs: access ( < 0)"); 
	if(s > ny+1) structerror("SegWgCrs: access ( > ny+1)");
	return wgvec[s]; 
}
*/

/* input from FILE dat */
void SegWgCrs::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		structerror("SegWgCrs: read");
	}
	else
	{
		if(wgvec != NULL) delete[] wgvec;
		wgvec = NULL;
		ny = fgetint(dat);
		if(ny<-1) structerror("read: ny < -1");
		hy.read(dat);
		wgvec = new Waveguide[ny+2];
		for(int i=0; i<=ny+1; ++i) (wgvec[i]).read(dat);
	}
	return;
}

/* output to FILE dat */
void SegWgCrs::write(FILE *dat) const
{
	comment("SegWgCrs", dat);
	comment("ny", dat); fputint(ny, dat);
	comment("hy", dat); hy.write(dat); 
	for(int i=0; i<=ny+1; ++i) (wgvec[i]).write(dat);
	return;
}

/* vacuum wavelength: should be unique in all segments */
double SegWgCrs::lambda()
{
	return wgvec[0].lambda;
}


/* check consistency of the SegWgCrs object: error&exit, if not ok */
void SegWgCrs::consistency() const
{
	double l;
	l = wgvec[0].lambda;
	for(int i=1; i<=ny+1; ++i)
	{
		if(wgvec[i].lambda !=l) 
			structerror("SegWgCrs: lambda not unique");
	}
	int sp; 
	sp = wgvec[0].special;
	for(int i=1; i<=ny+1; ++i)
	{
		if(wgvec[i].special !=sp) 
			structerror("SegWgCrs: wg.special not unique");
	}
	for(int i=0; i<=ny+1; ++i) wgvec[i].consistency(); 
	return;
}

/* get segment index corresponding to position y */ 
int SegWgCrs::segidx(double y) const
{
	int s=0;
	while(s<=ny && y>hy(s)) ++s;
	return s;
}

/* get segment boundaries corresponding to index s */  
Interval SegWgCrs::segment(int s) const
{
	Interval i;
//	if(s < 0) structerror("SegWgCrs: segment( < 0)"); 
//	if(s > ny+1) structerror("SegWgCrs: segment( > ny+1)");
	if(s <= 0) i.x0 = -AWAY;
	else i.x0 = hy(s-1);
	if(s >= ny+1) i.x1 = AWAY;
	else i.x1 = hy(s);
	return i;
}

/* get segment boundaries corresponding to position y */ 
Interval SegWgCrs::segment(double y) const
{
	return segment(segidx(y));
}

/* get refractive index at position x, y */
double SegWgCrs::n(double x, double y) const
{
	int s = segidx(y);
	Waveguide wg = wgvec[s];
	double nv = wg.n(wg.layeridx(x));
	if(wg.special) return sgn_sqrt(nv);
	return nv;
}

/* get permittivity at position x, y */
double SegWgCrs::eps(double x, double y) const
{
	int s = segidx(y);
	Waveguide wg = wgvec[s];
	double nv = wg.n(wg.layeridx(x));
	if(wg.special) return nv;
	return nv*nv;
}

/* translate: hy -> hy+dy */
void SegWgCrs::ytranslate(double dy)
{
	for(int s=0; s<=ny; ++s) hy(s) += dy;
	return;
}

/* translate: hx -> hx+dx for all segments */
void SegWgCrs::xtranslate(double dx)
{
	for(int s=0; s<=ny+1; ++s) (wgvec[s]).translate(dx); 
	return;
}

/* construct the corresponding vertical segmented waveguide sequence */
SegWgCrs SegWgCrs::rotate() const
{
	consistency();
	int b=0;
	for(int s=0; s<=ny+1; ++s) b += wgvec[s].nx+1; 
	Dvector tbd(b);
	Dvector bd;
	b=0;
	for(int s=0; s<=ny+1; ++s)
	{
		tbd.setsubvector(wgvec[s].hx, b);		
		b += wgvec[s].nx+1;
	}
	tbd.sort();
	bd = tbd;
	if(bd.nel<2) structerror("SegWgCrs: rotated: (1) #x_bounds < 2.");
	tbd(0) = bd(0); b=1;
	for(int i=1; i<=bd.nel-1; ++i)
	{
		if(fabs(tbd(b-1)-bd(i)) > HDIST)
		{
			tbd(b) = bd(i);
			++b;
		}
	}
	bd = tbd.subvector(b, 0);
	if(bd.nel<2) structerror("SegWgCrs: rotated: (2) #x_bounds < 2.");

	int nny = bd.nel-1;	

	SegWgCrs nst(nny);
	nst.hy = bd;
	Waveguide twg(ny);
	twg.lambda = wgvec[0].lambda;
	twg.hx = hy;
	twg.n.init(1.0);
	nst.init(twg);
	
	Interval itv;
	double x, y, r;
	for(int s=0; s<=nst.ny+1; ++s)
	{	
		itv = nst.segment(s);
		x = 0.5*(itv.x0+itv.x1);
		for(int l=0; l<=twg.nx+1; ++l)
		{
			itv = twg.layer(l);
	                y = 0.5*(itv.x0+itv.x1);
			int idx = wgvec[segidx(y)].layeridx(x);
			r = wgvec[segidx(y)].n(idx);
			nst(s).n(l) = r;
		}
	}
	for(int s=0; s<=nst.ny+1; ++s) nst(s).optimize();
	return nst;
}

/* determine a default window around the "inner" structure,
   e.g. for purposes of plotting */
void SegWgCrs::defaultwindow(Interval& wx, Interval& wy) const
{
	wx = wy = Interval(-1.0, 1.0);
	if(ny < -1) return;
	Waveguide wg = wgvec[0];
	double x0 = wg.hx(0);
	double x1 = wg.hx(wg.nx);
	double wl = wg.lambda;
	double emin = wg.defaultepseffmin();
	double xt;
	for(int j=1; j<=ny+1; ++j)
	{
		wg = wgvec[j];
		xt = wg.hx(0);
		if(xt < x0) x0 = xt;
		xt = wg.hx(wg.nx);
		if(xt > x1) x1 = xt;
		if(wg.eps(0) < emin) emin = wg.eps(0);
		if(wg.eps(wg.nx+1) < emin) emin = wg.eps(wg.nx+1);
	}
	if(emin < 1.0) emin=1.0;
	double d = 1.5*wl/sqrt(emin);
	x0 = x0-d;
	x1 = x1+d;
	double y0 = -d;
	double y1 = d;
	if(ny >= 0) { y0 = hy(0)-d; y1 = hy(ny)+d; }
	wx = Interval(x0, x1);
	wy = Interval(y0, y1);
	return;
}

/* convert this to a WgCrs object */
WgCrs SegWgCrs::wgcrs() const
{
	SegWgCrs h;
	if(ny >= 1) h = (*this);
	else
	{
		h = SegWgCrs(1);
		if(ny == -1) 
		{
			Waveguide wg = wgvec[0];
			h(0) = wg;
			h(1) = wg;
			h(2) = wg;	
			double d;
			d = wg.hx(wg.nx) - wg.hx(0);
			h.hy(0) = -d/2.0;
			h.hy(1) =  d/2.0;
		}
		else 
		{
			if(ny == 0) 
			{
				Waveguide wg = wgvec[0];
				h(0) = wg;
				h(1) = wg;
				h(2) = wgvec[1];
				double d;
				d = wg.hx(wg.nx) - wg.hx(0);
				h.hy(0) = hy(0)-d;
				h.hy(1) = hy(0);
			}
			else structerror("SegWgCrs: circuit: ny");
		}
	}
	SegWgCrs v;
	v = h.rotate();
	WgCrs circ(v.ny, h.ny);
	circ.lambda = h.wgvec[0].lambda; 
	circ.hx = v.hy;
	circ.hy = h.hy;
	Waveguide wg;
	double x;
	for(int s=0; s<=h.ny+1; ++s)
	{
		wg = h.wgvec[s];
		for(int l=0; l<=v.ny+1; ++l)
		{
			if(l>=1 && l<=v.ny) x = 0.5*(v.hy(l-1)+v.hy(l)); 
			else 
			{
				if(l==0) x = v.hy(0)-1.0;
				else     x = v.hy(v.ny)+1.0;
			}
			circ.n(l, s) = wg.n(wg.layeridx(x)); 
		}
	}
	return circ;
}

/* refractive index profile -> MATLAB .m file  
   ext0, ext1: filename id characters 
   xbeg, xend, ybeg, yend: display window in the x-y-plane */
void SegWgCrs::plot(double xbeg, double xend, double ybeg, double yend,
                       char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "n__.m";
	name[1] = ext0;
	name[2] = ext1;
	fprintf(stderr, "n [%g, %g] x [%g, %g] >> %s\n",  
	               xbeg, xend, ybeg, yend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "2D waveguide cross section");
	name[3] = 0;
	mlout_gengeoxy(dat, (*this), Rect(xbeg, ybeg, xend, yend));
	mlout_windowxy(dat, xbeg, xend, ybeg, yend);
	mlout_gengeorefindpatches(dat, (*this));
	mlout_wgcrsoverview(name, dat); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}
void SegWgCrs::plot(Rect disp, char ext0, char ext1) const
{
	plot(disp.x0, disp.x1, disp.y0, disp.y1, ext0, ext1);
	return;
}
void SegWgCrs::plot(char ext0, char ext1) const
{
	Interval wx, wy;
	defaultwindow(wx, wy);
	plot(wx.x0, wx.x1, wy.x0, wy.x1, ext0, ext1); 
	return;
}

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

/* for the VEIMS solver:
   cross section (2D) of a 3D rectangular waveguide,
   geometry & permittivity profile, vacuum wavelength;
   (so far) meant as a help to define a SegWgCrs */

/* initialize */
WgCrs::WgCrs()
{
	(*this) = WgCrs(1, 1);
}
WgCrs::WgCrs(int vnx, int vny)
{
	if(vnx < 1) structerror("WgCrs: vnx < 1");
	if(vny < 1) structerror("WgCrs: vny < 1");
	nx = vnx;
	hx = Dvector(nx+1);
	ny = vny;
	hy = Dvector(ny+1);
	n = Dmatrix(nx+2, ny+2);
}
WgCrs::WgCrs(int vnx, Dvector vhx,
             int vny, Dvector vhy, 
             double l, Dmatrix vn)
{
	nx = vnx;
	ny = vny;
	hx = vhx;
	hy = vhy;
	lambda = l;
	n = vn;
}

/* free allocated memory */
void WgCrs::strip()
{
	nx = 0;
	ny = 0;
	hx.strip();
	hy.strip();
	n.strip();
}

/* set Rectangle index (l, m) corresponding to position (x,y) */ 
void WgCrs::rectidx(double x, double y, int& l, int& m) const
{
	l = 0;
	while(l<=nx && x>hx(l)) ++l;
	m = 0;
	while(m<=ny && y>hy(m)) ++m;
	return;
}

/* get rectangle boundaries corresponding to index (l,m) */  
Rect WgCrs::rectbounds(int l, int m) const
{
	Rect r;
	if(l <= 0) r.x0 = -AWAY;
	else r.x0 = hx(l-1);
	if(l >= nx+1) r.x1 = AWAY;
	else r.x1 = hx(l);
	if(m <= 0) r.y0 = -AWAY;
	else r.y0 = hy(m-1);
	if(m >= ny+1) r.y1 = AWAY;
	else r.y1 = hy(m);
	return r;
}

/* get rectangle boundaries corresponding to position (x,y) */ 
Rect WgCrs::rectbounds(double x, double y) const
{
	int l, m;
	rectidx(x, y, l, m);
	return rectbounds(l, m);
}

/* test neighbourhood of two rectangles */
int WgCrs::testconnect(int l0, int m0, int l1, int m1) const
{
	int ld, md;
	ld = l1-l0; if(ld<0) ld = -ld;	
	md = m1-m0; if(md<0) md = -md;	
	if(ld==1 && md==0) return 1;
	if(ld==0 && md==1) return 1;
	if(ld==0 && md==0) return 1;
	return 0;
}

/* test matching of boundary positions with another waveguide */
int WgCrs::bdmatch(WgCrs w) const
{
	int j, k, e;
	for(j=0; j<=w.nx; ++j)
	{
		e = 0;
		for(k=0; k<=nx; ++k)
		{
			if(fabs(hx(k)-w.hx(j))<=HDIST) e = 1; 
		}
		if(e != 1) return 1;	
	}
	for(j=0; j<=w.ny; ++j)
	{
		e = 0;
		for(k=0; k<=ny; ++k)
		{
			if(fabs(hy(k)-w.hy(j))<=HDIST) e = 1; 
		}
		if(e != 1) return 1;	
	}
	return 0;
}

/* permittivity on rectangle l, m */
double WgCrs::eps(int l, int m) const
{
	double e;
	e = n(l,m);
	return e*e;
}

/* permittivity at position x, y */
double WgCrs::eps(double x, double y) const
{
	double e;
	int l, m;
	rectidx(x, y, l, m);
	e = n(l,m);
	return e*e;
}


/* translate: hx -> hx+dx, hy -> hy+dy */
void WgCrs::translate(double dx, double dy)
{
	int l, m;
	for(l=0; l<=nx; ++l) hx(l) += dx;
	for(m=0; m<=ny; ++m) hy(m) += dy;
	return;
}

/* refractive index profile -> MATLAB .m file  
   ext0, ext1: filename id characters 
   xbeg, xend, ybeg, yend: display window in the x-y-plane */
void WgCrs::plot(double xbeg, double xend, double ybeg, double yend,
                   char ext0, char ext1) const
{
	SegWgCrs s(*this);
	s.plot(xbeg, xend, ybeg, yend, ext0, ext1); 
	return;
}
void WgCrs::plot(Rect disp, char ext0, char ext1) const
{
	SegWgCrs s(*this);
	s.plot(disp, ext0, ext1); 
	return;
}
void WgCrs::plot(char ext0, char ext1) const
{
	SegWgCrs s(*this);
	s.plot(ext0, ext1); 
	return;
}

/* Rect: rectangular section from the waveguide cross section */
/* initialize */
Rect::Rect()
{
	x0 = 0.0;
	y0 = 0.0;
	x1 = 1.0;
	y1 = 1.0;
}
Rect::Rect(double xa, double ya, double xb, double yb)
{
	if(xa <= xb) { x0 = xa; x1 = xb; }
	else         { x0 = xb; x1 = xa; }
	if(ya <= yb) { y0 = ya; y1 = yb; }
	else         { y0 = yb; y1 = ya; }
}

/* output to FILE dat */
void Rect::write(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "[%g,%g|%g,%g]", x0, y0, x1, y1);
	}
	else     
	{
		fputdouble(x0, dat);
		fputdouble(y0, dat);
		fputdouble(x1, dat);
		fputdouble(y1, dat);
	}
	return;
}

/* input from FILE dat */
void Rect::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		structerror("Rect: read");
	}
	else 
	{
		x0 = fgetdouble(dat);
		y0 = fgetdouble(dat);
		x1 = fgetdouble(dat);
		y1 = fgetdouble(dat);
	}
	return;
}
 
Rect XYplane(-AWAY, -AWAY, AWAY, AWAY);

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

/* more general specification of a dielectric structure
   as a stack of possibly overlapping geometric shapes with constant
   permittivity, individual shapes */

/* initialize: generic */
Overlay::Overlay()
{
	type = HLAY;
	eps  = 1.0;
	n = 1.0;
	ubd = 1.0;
	lbd = -1.0;
	xo = 0.0; 
	zo = 0.0;
}

/* initialize: horizontal (HLAY) or vertical layer (VLAY),
   refr.ind. nc, thickness d, center at x=ofs (HLAY) or z=ofs (VLAY) */
Overlay::Overlay(OvlType t, double nc, double d, double ofs)
{
	if(t != HLAY && t != VLAY) structerror("Overlay: type");
	type = t;
	eps  = nc*nc;
	n = nc;
	ubd = ofs+0.5*d;
	lbd = ofs-0.5*d;
	xo = ofs; 
	zo = ofs;
}

/* initialize: RING cavity, 
   refr.ind. nc, radius r, core width d, origin at ox, oz */
Overlay::Overlay(OvlType t, double nc, 
                 double r, double d, double ox, double oz)
{
	if(t != RING) structerror("Overlay: type");
	type = t;
	eps = nc*nc;
	n = nc;
	ubd = r;
	lbd = r-d;
	xo = ox;
	zo = oz;
}

/* initialize: DISK cavity, 
   refr.ind. nc, radius R, origin at ox, oz */
Overlay::Overlay(OvlType t, double nc, double r, double ox, double oz)
{
	if(t != DISK) structerror("Overlay: type");
	type = t;
        eps = nc*nc;
	n = nc;
        ubd = r;
        xo = ox;
        zo = oz;
}

/* initialize: rectangular region RECT,
   defined by Intervals rx, rz, refr.ind. nc */
Overlay::Overlay(OvlType t, double nc, Interval rx, Interval rz)
{
	if(t != RECT) structerror("Overlay: type");
	type = t;
        eps = nc*nc;
	n = nc;
	if(rx.x0 <= rx.x1) { lbd = rx.x0; ubd = rx.x1; }
	else               { lbd = rx.x1; ubd = rx.x0; }
	if(rz.x0 <= rz.x1) { xo  = rz.x0; zo  = rz.x1; }
	else               { xo  = rz.x1; zo  = rz.x0; }
}

/* decide whether a position falls inside the shape (1) or not (0) */
int Overlay::inside(double x, double z)
{
	double dx, dz, r;
	switch(type)
	{
	case HLAY:
		if(x > ubd) return 0;
		if(x < lbd) return 0;
		return 1;
		break;
	case VLAY:
		if(z > ubd) return 0;
		if(z < lbd) return 0;
		return 1;
		break;
	case RING:
		dx = x-xo;
		dz = z-zo;
		r = sqrt(dx*dx+dz*dz);
		if(r > ubd) return 0;
		if(r < lbd) return 0;
		return 1;
	case DISK:
		dx = x-xo;
		dz = z-zo;
		r = sqrt(dx*dx+dz*dz);
		if(r > ubd) return 0;
		return 1;
	case RECT:
		if(x > ubd) return 0;
		if(x < lbd) return 0;
		if(z > zo) return 0;
		if(z < xo) return 0;
		return 1;
		break;
	}
	return 0;
}

/* read from FILE dat */
void Overlay::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
        {
                structerror("Overlay: read");
        }
        comment("Overlay", dat);
        type = (OvlType)fgetint(dat);
        eps = fgetdouble(dat);
        n = fgetdouble(dat);
        xo = fgetdouble(dat);
        zo = fgetdouble(dat);
        lbd = fgetdouble(dat);
        ubd = fgetdouble(dat);
        return;
}

/* output to FILE dat */
void Overlay::write(FILE *dat) const
{
        comment("Overlay", dat);
        comment("type", dat); fputint(((int)type), dat);
        comment("eps", dat); fputdouble(eps, dat);
        comment("n", dat); fputdouble(n, dat);
        comment("xo", dat); fputdouble(xo, dat);
        comment("zo", dat); fputdouble(zo, dat);
        comment("lbd", dat); fputdouble(lbd, dat);
        comment("ubd", dat); fputdouble(ubd, dat);
        return;
}

/* the entire stack */

/* initialize */
OvlStruct::OvlStruct()
{
	no = 0;
	lambda = 1.0;
	bgeps = 1.0;
	bgn = 1.0;
	ovlvec = NULL;
}
/* initialize */
OvlStruct::OvlStruct(double b, double l)
{
	no = 0;
	bgeps = b;
	bgn = sqrt(fabs(b));
	lambda = l;
	ovlvec = NULL;
}

/* destroy */
OvlStruct::~OvlStruct()
{
	if(ovlvec != NULL) delete[] ovlvec;
	ovlvec = NULL;
	no = 0; 
}

/* copy constructor */
OvlStruct::OvlStruct(const OvlStruct& o)
{
	no = o.no;
	bgeps = o.bgeps;
	bgn = o.bgn;
	lambda = o.lambda;
	if(no >= 1)
	{
		ovlvec = new Overlay[no];
		for(int i=0; i<=no-1; ++i) ovlvec[i] = o.ovlvec[i]; 
	}
	
}
	
/* assignment */
OvlStruct& OvlStruct::operator=(const OvlStruct& o)
{
	if(this != &o)
	{
		if(ovlvec != NULL) delete[] ovlvec;
		ovlvec = NULL;
		no = o.no; 
		bgeps = o.bgeps; 
		bgn = o.bgn; 
		lambda = o.lambda; 
		if(no >= 1)
		{
			ovlvec = new Overlay[no];
			for(int i=0; i<=no-1; ++i) ovlvec[i] = o.ovlvec[i]; 
		}

	}
	return *this;
}

/* delete all entries */
void OvlStruct::clear()
{
	no = 0; 
	if(ovlvec != NULL) delete[] ovlvec; 
	ovlvec = NULL;
}

/* place an Overlay on the top of the stack */
void OvlStruct::place(Overlay o)
{
	if(no<0) no = 0;
	Overlay *nov;
        nov = new Overlay[no+1];
        Overlay* np = nov;
	if(no >= 1)
	{
		Overlay* mp = ovlvec;
		for(int i=0; i<=no-1; ++i) *np++ = *mp++;
	}
        *np = o;
        if(ovlvec != NULL) delete[] ovlvec;
        ovlvec = nov;
        ++no;
        return;
}

/* input from FILE dat */
void OvlStruct::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		structerror("OvlStruct: read");
	}
	else
	{
		if(ovlvec != NULL) delete[] ovlvec;
		ovlvec = NULL;
		bgeps = fgetdouble(dat);
		bgn = fgetdouble(dat);
		lambda = fgetdouble(dat);
		no = fgetint(dat);
		if(no<0) structerror("OvlStruct: read: no < 0");
		if(no >= 1)
		{
			ovlvec = new Overlay[no];
			for(int i=0; i<=no-1; ++i) (ovlvec[i]).read(dat); 
		}
	}
	return;
}

/* output to FILE dat */
void OvlStruct::write(FILE *dat) const
{
	if(dat==stdout || dat==stderr)
	{
		fprintf(dat, 
		"OvlStruct, shapes: %d, n = %g, eps = %g, lambda = %g\n",
		        no, bgn, bgeps, lambda);
		for(int i=0; i<=no-1; ++i)
		{
		fprintf(dat, "[%d] type %d, eps = %g, n = %g\n", 
		i, ovlvec[i].type, ovlvec[i].eps, ovlvec[i].n);
		fprintf(dat, "     xo = %g, zo = %g, lbd = %g, ubd = %g\n", 
		ovlvec[i].xo, ovlvec[i].zo, ovlvec[i].lbd, ovlvec[i].ubd);
		}
		return;
	}
	comment("OvlStruct", dat);
	comment("bgeps", dat); fputdouble(bgeps, dat);
	comment("bgn", dat); fputdouble(bgn, dat);
	comment("lambda", dat); fputdouble(lambda, dat);
	comment("no", dat); fputint(no, dat);
	if(no >= 1)
	{
		for(int i=0; i<=no-1; ++i) (ovlvec[i]).write(dat); 
	}
	return;
}

/* permittivity at position x, z */
double OvlStruct::eps(double x, double z) const
{
	if(no >= 1)
	{
		for(int i=no-1; i>=0; --i)
		{
			if((ovlvec[i]).inside(x, z) != 0) return ovlvec[i].eps;
		}
	}
	return bgeps;
}

/* refractive index at position x, z */
double OvlStruct::n(double x, double z) const
{
	if(no >= 1)
	{
		for(int i=no-1; i>=0; --i)
		{
			if(ovlvec[i].inside(x, z) != 0) return ovlvec[i].n;
		}
	}
	return bgn;
}

/* for purely rectangular structures: convert this to a Circuit object */
Circuit OvlStruct::circuit() const
{	
	int nx, nz;
	Dvector hx, hz;	
	Circuit c;
	if(no == 0)
	{
		nx = 1; 
		hx = Dvector(nx+1);
		hx(0) = -0.5; hx(1) = 0.5;
		nz = 1;
		hz = Dvector(nz+1);
		hz(0) = -0.5; hz(1) = 0.5;
		c = Circuit(nx, nz);
		c.hx = hx;
		c.hz = hz;
		c.lambda = lambda;
		c.n = Dmatrix(nx+2, nz+2);
		c.n.init(bgn);
		c.special = 0;
		return c;
	}
	nx = 0; 
	nz = 0;
	hx.strip();
	hz.strip();
	for(int i=0; i<=no-1; ++i)
	{
		switch(ovlvec[i].type)
		{
		case HLAY:  
			hx.append(ovlvec[i].ubd);
			hx.append(ovlvec[i].lbd);
			break;
		case VLAY:  
			hz.append(ovlvec[i].ubd);
			hz.append(ovlvec[i].lbd);
			break;
		case RECT: 
			hx.append(ovlvec[i].ubd);
			hx.append(ovlvec[i].lbd);
			hz.append(ovlvec[i].zo);
			hz.append(ovlvec[i].xo);
			break;
		default:
			structerror(
				"OvlStruct: circuit: non rectangular Overlay");
				return c;
			break;
		}
	}
	Dvector thx = hx;
	hx.strip(); 
	thx.sort();
	for(int j=0; j<=thx.nel-1; ++j)
	{
		double x = thx(j);	 
		if(fabs(x) < AWAY/10.0)
		{
			if(hx.nel <= 0) hx.append(x);
			else if(fabs(x-hx(hx.nel-1)) > HDIST) hx.append(x);
		}
	}
	if(hx.nel <= 0)
	{
		hx = Dvector(2);
		hx(0) = -0.5;
		hx(1) =  0.5;
	}
	if(hx.nel == 1)
	{
		thx = hx;
		hx = Dvector(2);
		hx(0) = thx(0);
		hx(1) = thx(0)+1.0;
	}
	Dvector thz = hz;
	hz.strip(); 
	thz.sort();
	for(int j=0; j<=thz.nel-1; ++j)
	{
		double z = thz(j);	 
		if(fabs(z) < AWAY/10.0)
		{
			if(hz.nel <= 0) hz.append(z);
			else if(fabs(z-hz(hz.nel-1)) > HDIST) hz.append(z);
		}
	}
	if(hz.nel <= 0)
	{
		hz = Dvector(2);
		hz(0) = -0.5;
		hz(1) =  0.5;
	}
	if(hz.nel == 1)
	{
		thz = hz;
		hz = Dvector(2);
		hz(0) = thz(0);
		hz(1) = thz(0)+1.0;
	}
	nx = hx.nel-1;
	nz = hz.nel-1;
	c = Circuit(nx, nz);
	c.hx = hx;
	c.hz = hz;
	for(int i=0; i<=nx+1; ++i)
	{
		double x0, x1;
		if(i == 0   ) x0 = -AWAY; else x0 = hx(i-1); 
		if(i == nx+1) x1 =  AWAY; else x1 = hx(i); 
		for(int j=0; j<=nz+1; ++j)
		{
			double z0, z1;
			if(j == 0   ) z0 = -AWAY; else z0 = hz(j-1); 
			if(j == nz+1) z1 =  AWAY; else z1 = hz(j); 
			c.n(i, j) = n(0.5*(x0+x1), 0.5*(z0+z1));
		}
	}
	c.lambda = lambda;
	c.special = 0;
	return c;
}

/* for purely rectangular structures: convert this to a SegWgStruct object */
SegWgStruct OvlStruct::segwgstruct() const
{	
	Circuit c = circuit();
	SegWgStruct s(c);
	return s;
}