/*
 * WMM - Wave-matching method for mode analysis of dielectric waveguides
 * https://wmm.computational-photonics.eu/
 */

/*
 * mode.cpp
 * general mode definitions 
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmvis.h"
#include"mode.h"

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

/* store num values of component cp between (x0, y0) and (x1, y1) in a vector 
   cp: EX - HZ, SZ */
Dvector Mode::fieldvec(Fcomp cp, int num, 
                       double x0, double y0, double x1, double y1) const
{
	double x, y;
	double dx, dy;
	int j;
	if(num <= 0) modeerror("fieldvec: num <= 0");
	Dvector f(num);
	if(num == 1) 
	{
		x = 0.5*(x0+x1);
		y = 0.5*(x0+x1);
		f(0) = field(cp, x, y);
		return f;
	}
	dx = (x1-x0)/(num-1);
	dy = (y1-y0)/(num-1);
	for(j=0; j<=num-1; ++j)
	{
		x = x0+j*dx;	
		y = y0+j*dy;	
		f(j) = field(cp, x, y);
	}
	return f;
}

/* evaluate component cp on a rectangular npx x npy mesh on area disp 
   cp: EX - HZ, SZ 
   foa: ORG, MOD, SQR */
Dmatrix Mode::fieldmat(Rect disp, int npx, int npy, Fcomp cp, Afo foa) const
{
	double x, y;
	double dx, dy;
	int i, j;
	if(npx <= 1) modeerror("fieldmat: npx <= 1");
	if(npy <= 1) modeerror("fieldmat: npy <= 1");
	Dmatrix f(npx, npy);
	double ft;
	dx = (disp.x1-disp.x0)/(npx-1);
	dy = (disp.y1-disp.y0)/(npy-1);
	for(i=0; i<=npx-1; ++i)
	{
		x = disp.x0+i*dx;	
		for(j=0; j<=npy-1; ++j)
		{
			y = disp.y0+j*dy;	
			ft = field(cp, x, y);
			switch(foa)
			{
			case MOD: 
				ft = fabs(ft); 
				break;
			case SQR: 
				ft = ft*ft; 
				break;
			case ORG: 
				break;
			case REP: 
				if(cp == EZ || cp == HZ) ft = 0.0;
				break;
			case IMP: 
				if(cp != EZ && cp != HZ) ft = 0.0;
				break;
			}
			f(i, j) = ft;
		}
	}
	return f;
} 

/* z component of the Poyntingvector, integrated over the entire x-y-plane */
/* TE part */
double Mode::tepower() const
{
	return -0.5*recint(HX, EY, wg.xyplane);
}
/* TM part */
double Mode::tmpower() const
{
	return 0.5*recint(EX, HY, wg.xyplane);
}
/* total power */
double Mode::totpower() const
{
	return tepower()+tmpower();
}

/* polarization character of a hybrid mode */
/* mod: 'E': (int E_y^2 dxdy)/(int E_x^2+E_y^2 dxdy) */
/*      'H': (int H_x^2 dxdy)/(int H_x^2+H_y^2 dxdy) */
/*      'S': (int S_z^TE dxdy)/(int S_z dxdy) */
double Mode::polcharacter(char mod) const
{
	double te, tm;
	switch(mod)
	{
	case 'E':
		te = recint(EY, EY, wg.xyplane); 
		tm = recint(EX, EX, wg.xyplane); 
		break;
	case 'H':
		te = recint(HX, HX, wg.xyplane); 
		tm = recint(HY, HY, wg.xyplane); 
		break;
	default:
		te = tepower();
		tm = tmpower();
		break;
	}
	return te/(te+tm);
}

/* set maximum values for E, H, Sz (only a coarse approximation !) */
void Mode::setfieldmax()
{
	int l, m;
	double x, y, w, s;
	int np = 10;
	Rect r;
	minE = 1.0e+300;
	maxE = -1.0e+300;
	minH = 1.0e+300;
	maxH = -1.0e+300;
	minS = 1.0e+300;
	maxS = -1.0e+300;
	for(l=1; l<=wg.nx; ++l)
	{
		for(m=1; m<=wg.ny;++m)
		{
			r = wg.rectbounds(l, m);
			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(EX, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EY, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EZ, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(EX, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EY, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
				w = field(EZ, x, y); 
				if(w > maxE) maxE = w;
				if(w < minE) minE = w;
			}

			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(HX, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HY, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HZ, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(HX, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HY, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
				w = field(HZ, x, y); 
				if(w > maxH) maxH = w;
				if(w < minH) minH = w;
			}

			s = (r.x1-r.x0)/np;
			y = (r.y1+r.y0)/2.0;
			for(x=r.x0; x<=r.x1; x+=s)
			{
				w = field(SZ, x, y); 
				if(w > maxS) maxS = w;
				if(w < minS) minS = w;
			}
			s = (r.y1-r.y0)/np;
			x = (r.x1+r.x0)/2.0;
			for(y=r.y0; y<=r.y1; y+=s)
			{
				w = field(SZ, x, y); 
				if(w > maxS) maxS = w;
				if(w < minS) minS = w;
			}
		}
	}
	return;
}


/* coordinates of the maximum mode intensity on rectangle r, 
   average 1/e^2-radius in x- and y-direction 
   only a coarse approximation ! */
double Mode::maxintensity(Rect r, double& xm, double& ym, 
				  double& rx, double& ry) const
{
	double x, y, w, sx, sy;
	int np = 50;
	double max;
	double a, b;
	max = -1.0e+300;
	xm = 0.0;
	ym = 0.0;

	sx = (r.x1-r.x0)/np;
	sy = (r.y1-r.y0)/np;
	for(x=r.x0; x<=r.x1; x+=sx)
	{
		for(y=r.y0; y<=r.y1; y+=sy)
		{
			w = fabs(field(SZ, x, y)); 
			if(w > max)
			{
				max = w;
				xm = x;
				ym = y;
			}
		}
	}

	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		x -= sx;
		w = fabs(field(SZ, x, y)); 
	}
	a = x;
	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		x += sx;
		w = fabs(field(SZ, x, y)); 
	}
	b = x;
	rx = 0.5*(b-a);

	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		y -= sy;
		w = fabs(field(SZ, x, y)); 
	}
	a = y;
	x = xm; y = ym; w = fabs(field(SZ, x, y)); 
	while(w > max*0.1353)
	{
		y += sy;
		w = fabs(field(SZ, x, y)); 
	}
	b = y;
	ry = 0.5*(b-a);
	 
	return max;
}

/* normalize mode to totpower() == nrm */
void Mode::normalize(double nrm)
{
	ampfac = 1.0;
	ampfac = sqrt(nrm)/sqrt(totpower());
	setfieldmax();
	return;
}

/* write field profile  
	cp: EX - HZ, SZ
	foa: MOD, ORG, SQR  
	fac: output factor (typically 1.0 or -1.0)
	disp: output region on the x-y-plane 
	ext0, ext1: filename id characters
	npx, npy: number of points in output mesh
	adj: 1: adjust contour levels to cp
	adj: 0: adjust contour levels to minE, maxEm, minH ... 
*/
void Mode::writeprofile(Fcomp cp, Afo foa, double fac, Rect disp,
		   char ext0, char ext1, int npx, int npy, int adj) const
{
	FILE *dat;
	int xi, yi, l, m;
	double x, y, w, s;
	double minz, maxz;
	double mi, ma;
	double maxf;
	char name[13] = "t_______.hvl";
	Dmatrix f(npx+1,npy+1);
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = '-'; break;
	}
	switch(foa)
	{
		case MOD: name[3] = 'm'; break;
		case ORG: name[3] = 'f'; break;
		case SQR: name[3] = 'q'; break;
		case REP: name[3] = 'r'; break;
		case IMP: name[3] = 'i'; break;
	}
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	name[6] = ext0;
	name[7] = ext1;
	minz = 1.0e+300;
	maxz = -1.0e+300;
	maxf = -1.0e+300;
	for(xi=0; xi<=npx; ++xi)
	{
		x = disp.x0+(disp.x1-disp.x0)/npx*xi;
		for(yi=0; yi<=npy; ++yi)
		{
			y = disp.y0+(disp.y1-disp.y0)/npy*yi;
			w = field(cp, x, y);
			if(fabs(w) > maxf) maxf = fabs(w);
			w = fac*w;
			if(foa == MOD) w = fabs(w);
			if(foa == SQR) w = w*w;
			if(foa == IMP) w = 0.0;
			f(xi,yi) = w;
			if(w < minz) minz = w;
			if(w > maxz) maxz = w;
		}
	}
	fprintf(stderr, "   max|%c%c| = %g;  *(%g) >> %s\n",
		fldchr(cp), cpchr(cp), maxf, fac, name);
	if(adj == 0) 
	{
		mi = 0.0;
		ma = 1.0;
		switch(fldchr(cp))
		{
		case 'E': mi = minE; ma = maxE; break; 
		case 'H': mi = minH; ma = maxH; break; 
		case 'S': mi = minS; ma = maxS; break; 
		}
		switch(foa)
		{
		case MOD: 
			ma = fabs(ma);
			mi = fabs(ma);
			minz = 0.0;
			if(ma > mi) maxz = ma;
			else        maxz = mi;
			break;
		case ORG: case REP:
			maxz = ma;
			minz = mi;
			break;
		case SQR: 
			maxz = ma*ma; 
			if(maxz < mi*mi) maxz = mi*mi;	
			minz = 0.0;
			break;
		case IMP: 
			maxz = 1.0;
			minz = 0.0;
			break;
		}
	}
	dat = fopen(name, "w+");
	fprintf(dat, "c\n");
	fprintf(dat, "z\n");
	fprintf(dat, "%g\n", disp.y0);
	fprintf(dat, "%g\n", disp.x0);
	fprintf(dat, "%g\n", (disp.y1-disp.y0)/npy);
	fprintf(dat, "%g\n", (disp.x1-disp.x0)/npx);
	fprintf(dat, "%d\n", npy+1);
	fprintf(dat, "%d\n", npx+1);
	fprintf(dat, "v\n");
	fprintf(dat, "19\n");
	s = (maxz-minz)/20.0;
	for(xi=1; xi<=19; ++xi) fprintf(dat, "%g\n", minz+xi*s); 
	for(xi=0; xi<=npx; ++xi)
	{
		x = disp.x0+(disp.x1-disp.x0)/npx*xi;
		for(yi=0; yi<=npy; ++yi)
		{
			y = disp.y0+(disp.y1-disp.y0)/npy*yi;
			w = f(xi, yi);
			fprintf(dat, "%g\n", w);
		}
	}

	fprintf(dat, "%g\n", minz);
	fprintf(dat, "%g\n", maxz);

	fprintf(dat, "%d\n", wg.nx);
	fprintf(dat, "%d\n", wg.ny);
	for(l=0; l<=wg.nx; ++l)
		fprintf(dat, "%g\n", wg.hx(l));
	for(m=0; m<=wg.ny; ++m)
		fprintf(dat, "%g\n", wg.hy(m));
	for(l=0; l<=wg.nx+1; ++l)
	{
		for(m=0; m<=wg.ny+1; ++m)
			fprintf(dat, "%g\n", wg.n(l, m));
	}
	fprintf(dat, "%g\n", disp.x0);
	fprintf(dat, "%g\n", disp.x1);
	fprintf(dat, "%g\n", disp.y0);
	fprintf(dat, "%g\n", disp.y1);

	fclose(dat);
	return;
}

/* write profile section along line (x0,y0) -> (x1,y1) 
	cp: EX - HZ, SZ
	ext0, ext1: filename id characters
	np: number of output points
*/
void Mode::writesec(Fcomp cp, char ext0, char ext1, int np,
	            double x0, double y0, double x1, double y1) const
{
	FILE *dat;
	int i;
	double x, y, dx, dy;
	char name[13] = "t______.xyf";
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = '-'; break;
	}
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, ">> %s\n", name);
	dat = fopen(name, "w+");
	dx = x1-x0; dy = y1-y0;
	for(i=0; i<=np; ++i)
	{
		x = x0+i*dx/np;
		y = y0+i*dy/np;
		if(dx == 0.0) fprintf(dat, "%g  %g\n", y, field(cp, x, y));
		else
		{
			if(dy == 0.0) fprintf(dat, "%g  %g\n", 
						    x, field(cp, x, y));
			else
			{
			   fprintf(dat, "%g  %g\n", i*sqrt(dx*dx+dy*dy)/np,
			   				field(cp, x, y));
			}
		}
	}
	fclose(dat);
	return;
}

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

/* permittivity perturbation: propagation constant shift */
Complex Mode::phaseshift(Perturbation p) const
{
	Complex pt;
	double nrm;

	Complex db;
	db = CC0;

	// xx
	pt = p.e(0, 0); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EX, p.r); 	

	// yy
	pt = p.e(1, 1); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EY, p.r); 	

	// zz
	pt = p.e(2, 2); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EZ, EZ, p.r); 	

	// xy, yx
	pt = p.e(0, 1) + p.e(1, 0); 
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EY, p.r); 	

	// xz, zx
	pt = p.e(0, 2) - p.e(2, 0); 
	pt = CCI*pt;
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EZ, p.r); 	

	// yz, zy
	pt = p.e(1, 2) - p.e(2, 1); 
	pt = CCI*pt;
	if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EZ, p.r); 	
	
	nrm  = recint(HY, EX, wg.xyplane);
	nrm -= recint(HX, EY, wg.xyplane);
	nrm *= 2.0;
	db = db*(val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm);
	return db; 
}

/* geometry variations: propagation constant shift due to
   moving the horizontal boundary between rectangles [l,m] and [l+1,m]
   ( moving hx(l) for hy(m-1) < y < hy(m) ) */
double Mode::horgeovar(int l, int m) const
{
	double pt;
	double nrm;
	double y0, y1;
	if(l<0) return 0.0;
	if(m<0) return 0.0;
	if(l>wg.nx+1) return 0.0;
	if(m>wg.ny) return 0.0;
	if(m <= 0) y0 = wg.xyplane.y0;
	else       y0 = wg.hy(m-1);
	if(m >= wg.ny+1) y1 = wg.xyplane.y1;
	else             y1 = wg.hy(m);

/*
	if(pol == QTM)
	{
		pt = 0.5*(horlineint(l+1, m, HY, HY, wg.hx(l), y0, y1)
		         +horlineint(l, m,   HY, HY, wg.hx(l), y0, y1))
		   /wg.eps(l,m)/wg.eps(l+1,m); 
		pt += 0.5*(horlineint(l+1, m, EZ, EZ, wg.hx(l), y0, y1)
		          +horlineint(l, m,   EZ, EZ, wg.hx(l), y0, y1))
			   /val_invomep0(wg.lambda)/val_invomep0(wg.lambda) 
			   /beta/beta;
		pt *= (wg.eps(l,m)-wg.eps(l,m+1));
		nrm = 0.0;
		for(int li=0; li<=wg.nx+1; ++li)
		{
			for(int mi=0; mi<=wg.ny+1; ++mi)
				nrm += recintlm(li, mi, HY, HY)/wg.eps(li,mi);
		}
		pt *= beta/2.0/nrm;
		return pt;
	}
*/
/*
	if(pol == QTE)
	{
		pt = 0.5*(horlineint(l+1, m, EY, EY, wg.hx(l), y0, y1)
		         +horlineint(l, m,   EY, EY, wg.hx(l), y0, y1));
		pt *= (wg.eps(l,m)-wg.eps(l,m+1));
		nrm  = recint(EY, EY, wg.xyplane);
		pt *= val_k0(wg.lambda)*val_k0(wg.lambda)/2.0/beta/nrm;
		return pt;
	}
*/

	pt  = 0.5*horlineint(l+1, m, EX, EX, wg.hx(l), y0, y1)
	         *wg.eps(l+1,m)/wg.eps(l,m); 
	pt += 0.5*horlineint(l, m,   EX, EX, wg.hx(l), y0, y1)
	         *wg.eps(l,m)/wg.eps(l+1,m); 

	pt += 0.5*(horlineint(l+1, m, EY, EY, wg.hx(l), y0, y1) 
	          +horlineint(l, m,   EY, EY, wg.hx(l), y0, y1)); 

	pt += 0.5*(horlineint(l+1, m, EZ, EZ, wg.hx(l), y0, y1) 
	          +horlineint(l, m,   EZ, EZ, wg.hx(l), y0, y1)); 

	pt *= (wg.eps(l,m)-wg.eps(l+1,m));

	nrm  = recint(HY, EX, wg.xyplane);
	nrm -= recint(HX, EY, wg.xyplane);
	nrm *= 2.0;

	pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm;

	return pt;
}

/* geometry variations: propagation constant shift due to
   moving the vertical boundary between rectangles [l,m] and [l,m+1]
   ( moving hy(m) for hx(l-1) < x < hx(l) ) */
double Mode::vergeovar(int l, int m) const
{
	double pt;
	double nrm;
	double x0, x1;
	if(m<0) return 0.0;
	if(l<0) return 0.0;
	if(m>wg.ny) return 0.0;
	if(l>wg.nx+1) return 0.0;
	if(l == 0) x0 = wg.xyplane.x0;
	else       x0 = wg.hx(l-1);
	if(l == wg.nx+1) x1 = wg.xyplane.x1;
	else             x1 = wg.hx(l);

	if(pol == QTM)
	{
		pt = 0.5*(verlineint(l, m+1, HY, HY, x0, x1, wg.hy(m))
		         +verlineint(l, m,   HY, HY, x0, x1, wg.hy(m)));
		pt *= (wg.eps(l,m)-wg.eps(l,m+1));
		nrm = recint(HY, HY, wg.xyplane);
		pt *= val_k0(wg.lambda)*val_k0(wg.lambda)/2.0/beta/nrm;
		return pt;
	}

	pt  = 0.5*(verlineint(l, m+1, EX, EX, x0, x1, wg.hy(m)) 
	          +verlineint(l, m,   EX, EX, x0, x1, wg.hy(m))); 
	pt += 0.5*verlineint(l, m+1, EY, EY, x0, x1, wg.hy(m))
	         *wg.eps(l,m+1)/wg.eps(l,m); 
	pt += 0.5*verlineint(l, m,   EY, EY, x0, x1, wg.hy(m))
	         *wg.eps(l,m)/wg.eps(l,m+1); 
	pt += 0.5*(verlineint(l, m+1, EZ, EZ, x0, x1, wg.hy(m)) 
	          +verlineint(l, m,   EZ, EZ, x0, x1, wg.hy(m))); 
	pt *= (wg.eps(l,m)-wg.eps(l,m+1));

	nrm  = recint(HY, EX, wg.xyplane);
	nrm -= recint(HX, EY, wg.xyplane);
	nrm *= 2.0;

	pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm;

	return pt;
}

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

/* output to file */
void Mode::write(FILE *dat)
{
	writerdpid(dat);
	comment("MODE", dat);	
	comment("pol", dat); fputint(pol, dat);
	comment("sym", dat); fputint(sym, dat);
	comment("beta", dat); fputdouble(beta, dat);
	comment("k0", dat); fputdouble(k0, dat);
	comment("neff", dat); fputdouble(neff, dat);
	comment("npcB", dat); fputdouble(npcB, dat);
	comment("minE", dat); fputdouble(minE, dat);
	comment("maxE", dat); fputdouble(maxE, dat);
	comment("minH", dat); fputdouble(minH, dat);
	comment("maxH", dat); fputdouble(maxH, dat);
	comment("minS", dat); fputdouble(minS, dat);
	comment("maxS", dat); fputdouble(maxS, dat);
	comment("ampfac", dat); fputdouble(ampfac, dat);
	wg.write(dat);
	writerdp(dat);
	return;
}

/* input from file */
void Mode::read(FILE *dat)
{
	readrdpid(dat);
	pol = Polarization(fgetint(dat));
	sym = Symmetry(fgetint(dat));
	beta = fgetdouble(dat);
	k0 = fgetdouble(dat);
	neff = fgetdouble(dat);
	npcB = fgetdouble(dat);
	minE = fgetdouble(dat);
	maxE = fgetdouble(dat);
	minH = fgetdouble(dat);
	maxH = fgetdouble(dat);
	minS = fgetdouble(dat);
	maxS = fgetdouble(dat);
	ampfac = fgetdouble(dat);
	wg.read(dat);
	readrdp(dat);
	return;
}

/* output to default file */
void Mode::write_def(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	name[0] = 't';
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = '-'; break;
	}
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 'm';
	name[7] = 'o';
	name[8] = 'd';
	name[9] = 0;

	fprintf(stderr, "beta: %.10g,  neff: %.10g,  npcB: %.10g  >> %s\n",
			 beta, neff, npcB, name); 
	dat = fopen(name, "w+");
	write(dat);
	fclose(dat);
	return;
}

/* input from default file */
void Mode::read_def(Polarization p, Symmetry s, char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	name[0] = 't';
	switch(p)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(s)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = '-'; break;
	}
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 'm';
	name[7] = 'o';
	name[8] = 'd';
	name[9] = 0;

	fprintf(stderr, "<< %s : ", name);
	
	dat = fopen(name, "r");
	read(dat);
	fclose(dat);

	fprintf(stderr, "beta: %.10g,  neff: %.10g,  npcB: %.10g\n",
			 beta, neff, npcB); 
	return;
}

/* overlap 0.5*intint mode_EX*inif_HY - mode_EY*inif_HX dxdy with 
   an initial field, with real HX and HY components, integrals are 
   evaluated by gaussian quadrature formulas, separately on the
   fundamental rectangles, inside the domain rectangle r */

/* helper functions */
double Mode::ovleval(double (*inif)(Fcomp cp, double x, double y), 
		     double x, double y) const
{
	double fex, fhy, fey, fhx;
	fhy = inif(HY, x, y);
	fex = 0; if(fhy != 0.0) fex = field(EX, x, y);
	fhx = inif(HX, x, y);
	fey = 0; if(fhx != 0.0) fey = field(EY, x, y);
 	return 0.5*(fex*fhy-fey*fhx); 
}

/* based on: Numerical Recipes in C --- The Art of Scientific Computing
   Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */ 
double Mode::ovlpix(double (*inif)(Fcomp cp, double x, double y), 
                    double x0, double x1, double y, int numx) const
{
	int j;
	double xr,xm,dx,s;
	static double x[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	double ovl = 0.0; 
	double a, b, step;
	int i;
	step = (x1-x0)/((double) numx);
	for(i=0; i<=numx-1; ++i)
	{
		a = x0+((double) i)*step;
		b = a+step;
		xm=0.5*(b+a);
		xr=0.5*(b-a);
		s=0.0;
		for (j=1;j<=5;j++) {
			dx=xr*x[j];
			s += w[j]*(ovleval(inif, xm+dx, y)
				  +ovleval(inif, xm-dx, y));
		}
		ovl += s*xr;
	}
	return ovl;
}

double Mode::ovlpint(double (*inif)(Fcomp cp, double x, double y), Rect r,
		     int numx, int numy) const
{
	double ovl = 0.0;
	double a, b, step;
	int i;
	int j;
	double yr,ym,dy,s;
	static double y[]={0.0,0.1488743389,0.4333953941,
		0.6794095682,0.8650633666,0.9739065285};
	static double w[]={0.0,0.2955242247,0.2692667193,
		0.2190863625,0.1494513491,0.0666713443};
	step = (r.y1-r.y0)/((double) numy);
	for(i=0; i<=numy-1; ++i)
	{
		a = r.y0+((double) i)*step;
		b = a+step;
		ym=0.5*(b+a);
		yr=0.5*(b-a);
		s=0.0;
		for (j=1;j<=5;j++) {
			dy=yr*y[j];
			s += w[j]*(ovlpix(inif, r.x0, r.x1, ym+dy, numx)
				  +ovlpix(inif, r.x0, r.x1, ym-dy, numx));
		}
		ovl += s*yr;
	}
	return ovl;
}

#define HDIST 1.0e-8
double Mode::ovlp(double (*inif)(Fcomp cp, double x, double y), 
		  Rect r, int numx, int numy) const
{
	double ovl = 0.0; 
	double xb, yl;
	double xt, yr;
	int l, m, dum;
	Rect rp;
	if(r.x0 > r.x1) { xb=r.x1; r.x1=r.x0; r.x0=xb; }
	if(r.y0 > r.y1) { yl=r.y1; r.y1=r.y0; r.y0=yl; }
	yl = r.y0;
	while(fabs(r.y1-yl) > HDIST)
	{
		wg.rectidx(0.0, yl+HDIST/2.0, dum, m);
		yr = wg.rectbounds(0,m).y1;
		if(r.y1<yr) yr = r.y1;
		xb = r.x0;	
		while(fabs(r.x1-xb) > HDIST)
		{
			wg.rectidx(xb+HDIST/2.0, 0.0, l, dum);
			xt = wg.rectbounds(l,0).x1;
			if(r.x1<xt) xt = r.x1;
			rp.x0 = xb; rp.y0 = yl; rp.x1 = xt; rp.y1 = yr;
			ovl += ovlpint(inif, rp, numx, numy); 
			xb = xt;
		}
		yl = yr;
	}
	return ovl;
}

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

/* write single component to MATLAB .m file  
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR  
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters
   pltype: 'C': contour plot 
           'S': surface plot
           'I': intensity image 
           'N': field + mesh only, no plot commands (default) */
void Mode::mfile(Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		 char ext0, char ext1, char pltype) const
{
	FILE *dat;
	char name[13] = "t________.m";
	double minf, maxf;
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = 'n'; break;
	}
	switch(foa)
	{
		case MOD: name[3] = 'm'; break;
		case ORG: name[3] = 'f'; break;
		case SQR: name[3] = 'q'; break;
		case REP: name[3] = 'r'; break;
		case IMP: name[3] = 'i'; break;
	}
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	name[6] = ext0;
	name[7] = ext1;
	name[8] = pltype;
	fprintf(stderr, "%c%c >> %s\n", fldchr(cp), cpchr(cp), name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM mode profile");
	maxf = 1.0;
	minf = -1.0;
	switch(fldchr(cp))
	{
	case 'E': 
		switch(foa)
		{
		case ORG: case REP:
			minf = minE;
			maxf = maxE;
			break;
		case MOD:  
			maxf = fabs(maxE);
			if(fabs(minE)>maxf) maxf=fabs(minE);
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxE*maxE;
			if(minE*minE>maxf) maxf=minE*minE;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;
			break;
		}
		break;
	case 'H': 
		switch(foa)
		{
		case ORG: case REP:
			minf = minH;
			maxf = maxH;
			break;
		case MOD:  
			maxf = fabs(maxH);
			if(fabs(minH)>maxf) maxf=fabs(minH);
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxH*maxH;
			if(minH*minH>maxf) maxf=minH*minH;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;
			break;
		}
		break;
	case 'S': 
		switch(foa)
		{
		case ORG: case REP: case MOD:
			minf = minS;
			maxf = maxS;
			break;
		case SQR:  // not reasonable
			maxf = maxS*maxS;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;
			break;
		}
		break;
	}
	mlout_geo(dat, wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;		
	fld = fieldmat(disp, npx, npy, cp, foa);
	mlout_fld(dat, npx, npy, cp, fld);
	name[9] = 0; 
	switch(pltype)
	{
		case 'C': mlout_contour(name, dat, cp, foa); break;
		case 'S': mlout_surface(name, dat, cp, foa); break;
		case 'I': mlout_image(name, dat, cp, foa); break;
		default: break;
	}
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}


/* write all components to MATLAB .m file, surface plots  
   ext0, ext1: filename id characters
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh */
void Mode::acmfile(char ext0, char ext1, Rect disp, int npx, int npy) const
{
	FILE *dat;
	Dmatrix fld;		
	char name[13] = "t____A.m";
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = 'n'; break;
	}
	name[3] = ext0;
	name[4] = ext1;
	fprintf(stderr, "Ex - Hz >> %s\n", name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM mode profile");
	name[6] = 0;  
	mlout_geo(dat, wg, 0.0, 1.0);
	mlout_meshxy(dat, disp, npx, npy);
	fld = fieldmat(disp, npx, npy, EX, ORG);
	mlout_fld(dat, npx, npy, EX, fld);
	fld = fieldmat(disp, npx, npy, EY, ORG);
	mlout_fld(dat, npx, npy, EY, fld);
	fld = fieldmat(disp, npx, npy, EZ, ORG);
	mlout_fld(dat, npx, npy, EZ, fld);
	fld = fieldmat(disp, npx, npy, HX, ORG);
	mlout_fld(dat, npx, npy, HX, fld);
	fld = fieldmat(disp, npx, npy, HY, ORG);
	mlout_fld(dat, npx, npy, HY, fld);
	fld = fieldmat(disp, npx, npy, HZ, ORG);
	mlout_fld(dat, npx, npy, HZ, fld);

	mlout_acmfile(name, dat, minE, maxE, minH, maxH);
	mlout_print(dat, name, 'e');
	
	fclose(dat);
	return;
}

/* write profile section along line (x0,y0) -> (x1,y1) to MATLAB .m file
	cp: EX - HZ, SZ
	ext0, ext1: filename id characters
	np: number of output points
	pltype: 'L': geometry information + plot commands (default)
	        'V': field, mesh, and plot command, to be included into *L.m 
   --- at present implemented for horizontal or vertical lines only ---
*/
#define HDISTF 1.0e-8
void Mode::secmfile(Fcomp cp, char ext0, char ext1, int np, char pltype,
	            double x0, double y0, double x1, double y1) const
{
	FILE *dat;
	int i, j;
	double x, y, dx, dy, t;
	double xbd, ybd;
	double epsold, epsnew;
	double minf, maxf, f;
	int nsec, nbd;
	char name[13] = "t_______.m";
	char ori;
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = '-'; break;
	}
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	if(pltype != 'V') pltype = 'L';
	name[7] = pltype;

	dx = fabs(x1-x0);
	dy = fabs(y1-y0);
	
	if(dx > 1.0e-10 && dy > 1.0e-10)
	{
		fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> file \n", 
	                        fldchr(cp), cpchr(cp), x0, y0, x1, y1);
		fprintf(stderr, 
                " not implemented for tilted lines, use Mode.writesec().\n"); 
		return;
	}
	if(dx < 1.0e-10 && dy < 1.0e-10)
	{
		fprintf(stderr, "Mode.secmfile(): Nothing to do !\n");
		return;
	}
	if(dx < 1.0e-10) ori = 'h';
	else             ori = 'v';
	
	if(ori == 'v') { if(x0 > x1) { t = x0; x0 = x1; x1 = t; }}
	else           { if(y0 > y1) { t = y0; y0 = y1; y1 = t; }}

	fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> %s\n", 
	                fldchr(cp), cpchr(cp), x0, y0, x1, y1, name);

	Dmatrix fld(wg.nx+2+wg.ny+2, np+1);
	Dmatrix pos(wg.nx+2+wg.ny+2, np+1);
	Ivector nump(wg.nx+2+wg.ny+2);
	Dvector of(np+1);
	Dvector op(np+1);
	nsec = 0;
	Dvector bd(wg.nx+1+wg.ny+1);
	Dvector ri(wg.nx+2+wg.ny+2);
	nbd  = 0;

	dx = x1-x0; dy = y1-y0;
	nump(nsec) = 0;
	x = x0;
	y = y0;
	epsold = wg.eps(x, y);  	
	ri(nbd) = sqrt(epsold);
	if(ori=='h') pos(nsec, nump(nsec)) = y;  
	else         pos(nsec, nump(nsec)) = x;  
	fld(nsec, nump(nsec)) = field(cp, x, y);  
	++nump(nsec);
	for(i=1; i<=np; ++i)
	{
		x = x0+i*dx/np;
		y = y0+i*dy/np;
		epsnew = wg.eps(x, y);
		if(fabs(epsnew-epsold)<1.0e-10)
		{
			if(ori=='h') pos(nsec, nump(nsec)) = y;  
			else         pos(nsec, nump(nsec)) = x;  
			fld(nsec, nump(nsec)) = field(cp, x, y);  
			++nump(nsec);
		}
		else
		{
			if(ori=='h') { xbd = x0; 
                                       ybd = (wg.rectbounds(x, y)).y0; }
			else         { xbd = (wg.rectbounds(x, y)).x0; 
                                       ybd = y0; }
			if(ori=='h') pos(nsec, nump(nsec)) = ybd;  
			else         pos(nsec, nump(nsec)) = xbd;  
			fld(nsec, nump(nsec)) = 	
				field(cp, xbd-dx*HDISTF, ybd-dy*HDISTF);  
			++nump(nsec);
			if(ori=='h') bd(nbd) = ybd; 
			else         bd(nbd) = xbd; 
			++nbd;
			epsold = epsnew;
			ri(nbd) = sqrt(epsold);
			++nsec;
			nump(nsec) = 0;
			if(ori=='h') pos(nsec, nump(nsec)) = ybd;  
			else         pos(nsec, nump(nsec)) = xbd;  
			fld(nsec, nump(nsec)) = 	
				field(cp, xbd+dx*HDISTF, ybd+dy*HDISTF);  
			++nump(nsec);
			if(ori=='h') pos(nsec, nump(nsec)) = y;  
			else         pos(nsec, nump(nsec)) = x;  
			fld(nsec, nump(nsec)) = field(cp, x, y);  
			++nump(nsec);
		}
	}
	++nsec;

	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM mode field section");
	name[8] = 0;
	minf = -1.0;
	maxf = 1.0;
	if(pltype == 'L')
	{
		switch(fldchr(cp))
		{
			case 'E': minf = minE; maxf = maxE; break;
			case 'H': minf = minH; maxf = maxH; break;
			case 'S': minf = minS; maxf = maxS; break;
		}
		mlout_geo(dat, wg, minf, maxf); 
		mlout_Lsecgeo(dat, x0, y0, x1, y1, nbd, bd, ri);
	}
	minf = fld(0, 0);
	maxf = fld(0, 0);
	
	for(j=0; j<=nsec-1; ++j)
	{
		for(i=0; i<=nump(j)-1; ++i)
		{
			f = fld(j, i);
			if(f < minf) minf = f;
			if(f > maxf) maxf = f;
			of(i) = f;
			op(i) = pos(j, i);
		}
		if(pltype == 'L')
			mlout_sec1D(dat, cp, dig10(j), dig1(j), ' ', ' ', 
                                    nump(j), of, op); 
		else
			mlout_sec1D(dat, cp, dig10(j), dig1(j), ext0, ext1, 
                                    nump(j), of, op); 
	}
	if(pltype == 'L') 
	{
		mlout_Lsecplot(name, dat, ori, minf, maxf, cp, nbd, nsec); 
		mlout_print(dat, name, 'e');
	}
	else mlout_Vsecplot(dat, cp, nsec, ext0, ext1); 
	fclose(dat);
	return;
}


/* write single component to MATLAB .m file, fancy style :-)  
   cp: EX - HZ, SZ
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters */
void Mode::fancymfile(Fcomp cp, Rect disp, int npx, int npy, 
                      char ext0, char ext1) const
{
	FILE *dat;
	int np, l, m;
	double x0, x1, xp, y0, y1, yp;
	int numc;
	char name[13] = "t______F.m";
	switch(pol)
	{
		case QTE: name[1] = 'e'; break;
		case QTM: name[1] = 'm'; break;
		case VEC: name[0] = 'v';
			  name[1] = 'c'; break;
	}
	switch(sym)
	{
		case SYM: name[2] = 's'; break;
		case ASY: name[2] = 'a'; break;
		case NOS: name[2] = 'n'; break;
	}
	name[3] = fldchr(cp);
	name[4] = cpchr(cp);
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "%c%c >> %s\n", fldchr(cp), cpchr(cp), name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM mode profile :-)");
	name[8] = 0; 

	switch(fldchr(cp))
	{
		case 'E': mlout_geo(dat, wg, minE, maxE); break;
		case 'H': mlout_geo(dat, wg, minH, maxH); break;
		case 'S': mlout_geo(dat, wg, minS, maxS); break;
	}
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = fieldmat(disp, npx, npy, cp, ORG);
	mlout_fld(dat, npx, npy, cp, fld);

	Dvector fv;
	numc = 0;
	for(l=0; l<=wg.nx+1; ++l)
	{
		if(l==0) x0 = disp.x0; 	
                else     x0 = wg.hx(l-1);
		if(l==wg.nx+1) x1 = disp.x1; 	
                else           x1 = wg.hx(l);
		for(m=0; m<=wg.ny; ++m)
		{
			if(fabs(wg.n(l,m)-wg.n(l,m+1)) > 1.0e-10)
			{
				yp = wg.hy(m);
				np = (int)(((double) npx)
                                     *(x1-x0)/(disp.x1-disp.x0)); 
				if(np >= 2)
				{
				fv = fieldvec(cp, np, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				fv = fieldvec(cp, np, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	for(m=0; m<=wg.ny+1; ++m)
	{
		if(m==0) y0 = disp.y0; 	
                else     y0 = wg.hy(m-1);
		if(m==wg.ny+1) y1 = disp.y1; 	
                else           y1 = wg.hy(m);
		for(l=0; l<=wg.nx; ++l)
		{
			if(fabs(wg.n(l,m)-wg.n(l+1,m)) > 1.0e-10)
			{
				xp = wg.hx(l);
				np = (int)(((double)npy)
                                     *(y1-y0)/(disp.y1-disp.y0)); 
				if(np >= 2)
				{
				fv = fieldvec(cp, np, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST); 
				mlout_sec(dat, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				fv = fieldvec(cp, np, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST);
				mlout_sec(dat, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST, 
					  np, cp, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	fv = fieldvec(cp, npx, disp.x0, disp.y0, disp.x1, disp.y0);
	mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, 
                  npx, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(cp, npy, disp.x1, disp.y0, disp.x1, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, 
                  npy, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(cp, npx, disp.x1, disp.y1, disp.x0, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, 
                  npx, cp, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(cp, npy, disp.x0, disp.y1, disp.x0, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, 
                  npy, cp, dig10(numc), dig1(numc), fv);
	++numc;

	mlout_fancy(name, dat, cp, numc);  	
	mlout_print(dat, name, 'p');

	fclose(dat);
	return;
}