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

/*
 * fim.cpp
 * guided modes of circular multi-step-index optical fibers 
 */

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

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

/* Cylinder functions & derivatives, integer order o, real argument a, 
   wrapper functions: warn at computation failure */
int FIMCylEvalError = 0; // error indicator 
double FIMCylSPV = 1.0e-15; // substituted in case of zero a 

double BfJ(int o, double a)
{
	FIMCylEvalError = 0;
	if(a < 0.0) a = 0.0;
	if(o >= 0) return bessel_J(o, a);
	if((-o)%2 == 1) return -bessel_J(-o, a);
	return bessel_J(-o, a);
}

double BfY(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) { FIMCylEvalError = 1; fprintf(stderr, "Y_(%d)(%g) failed.\n", o, a); return 0.0; }
	if(o >= 0) return bessel_Y(o, a);
	if((-o)%2 == 1) return -bessel_Y(-o, a);
	return bessel_Y(-o, a);
}

double BfI(int o, double a)
{
	FIMCylEvalError = 0;
	if(a < 0.0) a = 0.0;
	if(o >= 0) return bessel_I(o, a);
	return bessel_I(-o, a);
}

double BfK(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) { FIMCylEvalError = 1; fprintf(stderr, "K_(%d)(%g) failed.\n", o, a); return 0.0; }
	if(o >= 0) return bessel_K(o, a);
	return bessel_K(-o, a);
}

double dBfJ(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) a = FIMCylSPV;
	if(o == 0) return -bessel_J1(a);
	if(o < 0)
	{
		double s = 1.0;
		int oa = abs(o);
		if(oa%2 == 1) s = -1.0; 
		return s*(bessel_J(oa-1, a)-bessel_J(oa, a)*oa/a);
	}
	return bessel_J(o-1, a)-bessel_J(o, a)*o/a;
}

double dBfY(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) { FIMCylEvalError = 1; fprintf(stderr, "dY_(%d)(%g) failed.\n", o, a); return 0.0; }
	if(o == 0) return -bessel_Y1(a);
	if(o < 0)
	{
		double s = 1.0;
		int oa = abs(o);
		if(oa%2 == 1) s = -1.0; 
		return s*(bessel_Y(oa-1, a)-bessel_Y(oa, a)*oa/a);
	}
	return bessel_Y(o-1, a)-bessel_Y(o, a)*o/a;
}

double dBfI(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) a = FIMCylSPV;
	if(o == 0) return bessel_I1(a);
	if(o < 0)
	{
		int oa = abs(o);
		return bessel_I(oa+1, a)+bessel_I(oa, a)*oa/a;
	}
	return bessel_I(o+1, a)+bessel_I(o, a)*o/a;
}

double dBfK(int o, double a)
{
	FIMCylEvalError = 0;
	if(a <= 0.0) { FIMCylEvalError = 1; fprintf(stderr, "dK_(%d)(%g) failed.\n", o, a); return 0.0; }
	if(o == 0) return -bessel_K1(a);
	if(o < 0)
	{
		int oa = abs(o);
		return -bessel_K(oa+1, a)+bessel_K(oa, a)*oa/a;
	}
	return -bessel_K(o+1, a)+bessel_K(o, a)*o/a;
}

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

int FIMsP_TrcFail = 0;
/* evaluate transverse resonance condition, set cA, cB, cC, cD accordingly */
/* l= 0, TE */
double FIMode::travresTE()
{
	FIMsP_TrcFail = 0;
	double r;
	double rho;
	double lon;
	double kpa;
	double bfc, bfd;
	double dfc, dfd;
	double em, ep;
	Dmatrix T;
	if(fillc >= 2)
	{
		cA.init(0.0);
		cB.init(0.0);
		cC(wg.nx+1) = 0.0;
		cD(wg.nx+1) = cDe;
		lon = wg.n(wg.nx+1);
		kappa(wg.nx+1) = k0*sqrt(neff*neff-lon*lon);
	}
	T.unit(2);
	for(int l=wg.nx; l>=0; --l)
	{	
		r = R+wg.hx(l); 
		em = wg.n(l)  *wg.n(l);
		ep = wg.n(l+1)*wg.n(l+1);
		
		lon = wg.n(l);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfc =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfd =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfd = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfc =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfd =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfd = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		Dmatrix O(2, 2);
		O(0, 0) = bfc; 
		O(0, 1) = bfd;
		O(1, 0) = (k0*k0*ep-beta*beta)*dfc;
		O(1, 1) = (k0*k0*ep-beta*beta)*dfd;

		if(fillc >= 2)
		{
			kappa(l) = kpa;
		}
		
		lon = wg.n(l+1);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfc =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfd =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfd = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfc =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfd =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfd = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		Dmatrix P(2, 2);
		P(0, 0) = bfc; 
		P(0, 1) = bfd;
		P(1, 0) = (k0*k0*em-beta*beta)*dfc;
		P(1, 1) = (k0*k0*em-beta*beta)*dfd;

		O.inverse();
		Dmatrix L;
		L = mult(O, P);
		T = mult(L, T);

		if(fillc >= 2)
		{
			Dvector vm(2), vp(2);
			vp(0) = cC(l+1); vp(1) = cD(l+1); 
			vm = mult(L, vp);
			cC(l) = vm(0); cD(l) = vm(1); 
// fprintf(stderr, "[l = %d]\n", l);
// fprintf(stderr, "   A: %.10g\n", cA(l));
// fprintf(stderr, "   B: %.10g\n", cB(l));
// fprintf(stderr, "   C: %.10g\n", cC(l));
// fprintf(stderr, "   D: %.10g\n", cD(l));
		}
	}
	double tr = T(1, 1);

	if(fillc <= 0) return tr;

	if(fillc == 1)
	{
		cBe = 0.0;
		cDe = 1.0; 
// fprintf(stderr, "Be: %.10g\n", cBe);
// fprintf(stderr, "De: %.10g\n", cDe);
		return tr;
	}
	// fillc >= 2
	cD(0) = 0.0;
	return tr;
}
/* l= 0, TM */
double FIMode::travresTM()
{
	FIMsP_TrcFail = 0;
	double r;
	double rho;
	double lon;
	double kpa;
	double bfa, bfb;
	double dfa, dfb;
	double em, ep;
	Dmatrix T;
	if(fillc >= 2)
	{
		cA(wg.nx+1) = 0.0;
		cB(wg.nx+1) = cBe;
		cC.init(0.0);
		cD.init(0.0);
		lon = wg.n(wg.nx+1);
		kappa(wg.nx+1) = k0*sqrt(neff*neff-lon*lon);
	}
	T.unit(2);
	for(int l=wg.nx; l>=0; --l)
	{	
		r = R+wg.hx(l); 
		em = wg.n(l)  *wg.n(l);
		ep = wg.n(l+1)*wg.n(l+1);
		
		lon = wg.n(l);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfa =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfa = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfa =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfa = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		Dmatrix O(2, 2);
		O(0, 0) = bfa; 
		O(0, 1) = bfb; 
		O(1, 0) = em*(k0*k0*ep-beta*beta)*dfa;
		O(1, 1) = em*(k0*k0*ep-beta*beta)*dfb;

		if(fillc >= 2)
		{
			kappa(l) = kpa;
		}
		
		lon = wg.n(l+1);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfa =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfa = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfa =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfa = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
		}
		Dmatrix P(2, 2);
		P(0, 0) = bfa; 
		P(0, 1) = bfb; 
		P(1, 0) = ep*(k0*k0*em-beta*beta)*dfa;
		P(1, 1) = ep*(k0*k0*em-beta*beta)*dfb;

		O.inverse();
		Dmatrix L;
		L = mult(O, P);
		T = mult(L, T);

		if(fillc >= 2)
		{
			Dvector vm(2), vp(2);
			vp(0) = cA(l+1); vp(1) = cB(l+1); 
			vm = mult(L, vp);
			cA(l) = vm(0); cB(l) = vm(1);
// fprintf(stderr, "[l = %d]\n", l);
// fprintf(stderr, "   A: %.10g\n", cA(l));
// fprintf(stderr, "   B: %.10g\n", cB(l));
// fprintf(stderr, "   C: %.10g\n", cC(l));
// fprintf(stderr, "   D: %.10g\n", cD(l));
		}
	}
	double tr = T(1, 1);

	if(fillc <= 0) return tr;

	if(fillc == 1)
	{
		cBe = 1.0;
		cDe = 0.0; 
// fprintf(stderr, "Be: %.10g\n", cBe);
// fprintf(stderr, "De: %.10g\n", cDe);
		return tr;
	}
	// fillc >= 2
	cB(0) = 0.0;
	return tr;
}
/* l>= 1, hybrid modes */
double FIMode::travres()
{
if(aol == 0)
{
	if(pid == 'e') return travresTE();
	if(pid == 'm') return travresTM();
	fimerror("travres: aol=0, confusing pid");
	return 0.0;
}
	FIMsP_TrcFail = 0;
	double r;
	double rho;
	double lon;
	double kpa;
	double bfa, bfb, bfc, bfd;
	double dfa, dfb, dfc, dfd;
	double em, ep;
	Dmatrix T;
	if(fillc >= 2)
	{
		cA(wg.nx+1) = 0.0;
		cB(wg.nx+1) = cBe;
		cC(wg.nx+1) = 0.0;
		cD(wg.nx+1) = cDe;
		lon = wg.n(wg.nx+1);
		kappa(wg.nx+1) = k0*sqrt(neff*neff-lon*lon);
	}
	T.unit(4);
	for(int l=wg.nx; l>=0; --l)
	{	
		r = R+wg.hx(l); 
		em = wg.n(l)  *wg.n(l);
		ep = wg.n(l+1)*wg.n(l+1);
		
		lon = wg.n(l);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfa =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfc = bfa;
			bfd = bfb;
			dfa = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = dfa;
			dfd = dfb;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfa =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfc = bfa;
			bfd = bfb;
			dfa = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = dfa;
			dfd = dfb;
		}
		Dmatrix O(4, 4);
		O(0, 0) = bfa; 
		O(0, 1) = bfb; 
		O(0, 2) = 0.0; 
		O(0, 3) = 0.0;
		O(1, 0) = 0.0; 
		O(1, 1) = 0.0; 
		O(1, 2) = bfc; 
		O(1, 3) = bfd;
		O(2, 0) = em*(k0*k0*ep-beta*beta)*dfa;
		O(2, 1) = em*(k0*k0*ep-beta*beta)*dfb;
		O(2, 2) = -ommu0*ep*beta*aol/r*bfc;
		O(2, 3) = -ommu0*ep*beta*aol/r*bfd;
		O(3, 0) = -omep0*ep*beta*aol/r*bfa;
		O(3, 1) = -omep0*ep*beta*aol/r*bfb;
		O(3, 2) = (k0*k0*ep-beta*beta)*dfc;
		O(3, 3) = (k0*k0*ep-beta*beta)*dfd;

		if(fillc >= 2)
		{
			kappa(l) = kpa;
		}
		
		lon = wg.n(l+1);
		if(lon > neff)
		{
			kpa = k0*sqrt(lon*lon - neff*neff);
			rho = kpa*r;
			bfa =      BfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfc = bfa;
			bfd = bfb;
			dfa = kpa*dBfJ(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfY(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = dfa;
			dfd = dfb;
		}
		else 
		{
			kpa = k0*sqrt(neff*neff-lon*lon);
			rho = kpa*r;
			bfa =      BfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfb =      BfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			bfc = bfa;
			bfd = bfb;
			dfa = kpa*dBfI(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfb = kpa*dBfK(aol, rho); if(FIMCylEvalError) FIMsP_TrcFail = 1;
			dfc = dfa;
			dfd = dfb;
		}
		Dmatrix P(4, 4);
		P(0, 0) = bfa; 
		P(0, 1) = bfb; 
		P(0, 2) = 0.0; 
		P(0, 3) = 0.0;
		P(1, 0) = 0.0; 
		P(1, 1) = 0.0; 
		P(1, 2) = bfc; 
		P(1, 3) = bfd;
		P(2, 0) = ep*(k0*k0*em-beta*beta)*dfa;
		P(2, 1) = ep*(k0*k0*em-beta*beta)*dfb;
		P(2, 2) = -ommu0*em*beta*aol/r*bfc;
		P(2, 3) = -ommu0*em*beta*aol/r*bfd;
		P(3, 0) = -omep0*em*beta*aol/r*bfa;
		P(3, 1) = -omep0*em*beta*aol/r*bfb;
		P(3, 2) = (k0*k0*em-beta*beta)*dfc;
		P(3, 3) = (k0*k0*em-beta*beta)*dfd;

		O.inverse();
		Dmatrix L;
		L = mult(O, P);
		T = mult(L, T);

		if(fillc >= 2)
		{
			Dvector vm(4), vp(4);
			vp(0) = cA(l+1); vp(1) = cB(l+1); vp(2) = cC(l+1); vp(3) = cD(l+1);
			vm = mult(L, vp);
			cA(l) = vm(0); cB(l) = vm(1); cC(l) = vm(2); cD(l) = vm(3);
// fprintf(stderr, "[l = %d]\n", l);
// fprintf(stderr, "   A: %.10g\n", cA(l));
// fprintf(stderr, "   B: %.10g\n", cB(l));
// fprintf(stderr, "   C: %.10g\n", cC(l));
// fprintf(stderr, "   D: %.10g\n", cD(l));
		}
	}
	double tr = T(1, 1)*T(3, 3) - T(3, 1)*T(1, 3);

	if(fillc <= 0) return tr;

	if(fillc == 1)
	{
		Dmatrix M(2, 2);
		M(0, 0) = T(1, 1); M(0, 1) = T(1, 3); 
		M(1, 0) = T(3, 1); M(1, 1) = T(3, 3); 

		double v0, v1;
		double d = M(0, 0)*M(0, 0) + 4.0*M(0, 1)*M(1, 0) - 2.0*M(0, 0)*M(1, 1) + M(1, 1)*M(1, 1);
		if(d > 0)
		{
			d = sqrt(d);
			double l0, l1;
			l0 = 0.5*(M(0, 0)+M(1, 1)-d);
			l1 = 0.5*(M(0, 0)+M(1, 1)+d);
			if(fabs(l0) < fabs(l1)) { v0 = -M(0, 1); v1 = M(0, 0)-l0; }
			else                    { v0 = -M(0, 1); v1 = M(0, 0)-l1; }
		}
		else 
		{
			if(fabs(M(1, 1)*M(1, 0)) > fabs(M(0, 1)*M(0, 0))) { v0 = -M(1, 1); v1 = M(1, 0); }
			else                                              { v0 = -M(0, 1); v1 = M(0, 0); }
		}

		double nrm = 1.0/sqrt(v0*v0+v1*v1);
		cBe = v0*nrm;
		cDe = v1*nrm; 
// fprintf(stderr, "Be: %.10g\n", cBe);
// fprintf(stderr, "De: %.10g\n", cDe);
		return tr;
	}
	// fillc >= 2
	cB(0) = 0.0;
	cD(0) = 0.0;
	return tr;
}
/* ...
   n: trial value for the effective index, 
   f1: fill in profile coefficients */
double FIMode::travres(double n, int f)
{
	neff = n;
	beta = n*k0;
	if(f == 0)
	{ 
		fillc = 0;
		return travres();
	}
	fillc = 1;
	travres();
	fillc = 2;
	return travres();
}

/* minimum effective index distance from local refractive indices */
double FIMsP_MinNDiff = 1.0e-4;

/* convergence to TRC roots: effective index accuray */
double FIMsP_NeffAcc = 1.0e-10;

/* minimum effective index distance for modes of the same angular order
   to be considered different */
double FIMsP_NeffSep = 1.0e-8;


/* trap a root in the transverse resonance condition (bisection) */
bool rootin(double a, double b)
{
	if(a > 0.0 && b < 0.0) return true;
	if(a < 0.0 && b > 0.0) return true;
	if(a == 0.0) return true;
	if(b == 0.0) return true;
	return false;
}
double FIMode::converge(double n0, double n1)
{
	double t0, t1, t;
	double n;
	if(n0 >= n1) { n=n0; n0=n1; n1=n; }
	t0 = travres(n0, 0);
	t1 = travres(n1, 0);
	while(fabs(n1-n0) > FIMsP_NeffAcc)
	{
		n = 0.5*(n0+n1);
		t = travres(n, 0);
// fprintf(stderr, "n %g:%g:%g   t %g:%g:%g\n", n0, n, n1, t0, t, t1);
		     if(rootin(t0, t)) n1 = n;
		else if(rootin(t, t1)) n0 = n;
		else fimerror("converge: confused");
	}
	return 0.5*(n0+n1);
}

/* radial profile & derivatives, principal components cp ~~ EZ, HZ */
void FIMode::profile(double r, double &e, double &de, double &h, double &dh) const 
{
	int l = wg.layeridx(r-R);
	if(disc_prof == 0 || r <= d_ri || r >= d_re)
	{
		double rho = kappa(l)*r;
		 e = 0.0;
		de = 0.0;
		 h = 0.0;
		dh = 0.0;
		if(wg.n(l) >= neff)
		{
			if(l > wg.nx) return;
			double  J =  BfJ(aol, rho);
			double dJ = dBfJ(aol, rho);
			 e += cA(l)* J;
			de += cA(l)*dJ;
			 h += cC(l)* J;
			dh += cC(l)*dJ;
			if(l >= 1)
			{
				double  Y =  BfY(aol, rho);
				double dY = dBfY(aol, rho);
				 e += cB(l)* Y;
				de += cB(l)*dY;
				 h += cD(l)* Y;
				dh += cD(l)*dY;
			}
		}
		else 
		{
			if(l <= wg.nx)
			{
				double  I =  BfI(aol, rho);
				double dI = dBfI(aol, rho);
				 e += cA(l)* I;
				de += cA(l)*dI;
				 h += cC(l)* I;
				dh += cC(l)*dI;
			}
			if(l >= 1)
			{
				double  K =  BfK(aol, rho);
				double dK = dBfK(aol, rho);
				 e += cB(l)* K;
				de += cB(l)*dK;
				 h += cD(l)* K;
				dh += cD(l)*dK;
			}
		}
		de = de*kappa(l);
		dh = dh*kappa(l);
		return;
	}
	int klo = 0;
	int khi = n_rad(l)-1;
	while(khi-klo > 1) 
	{
		int k = (khi+klo) >> 1;
		if(v_rad(l, k) > r) khi = k;
		else                klo = k;
	}
	double s = v_rad(l, khi)-v_rad(l, klo);
	double a = (v_rad(l, khi)-r)/s;
	double b = (r-v_rad(l, klo))/s;
	 e = v_prfE(l, klo)*a + v_prfE(l, khi)*b + ((a*a*a-a)*d_prfE(l, klo)+(b*b*b-b)*d_prfE(l, khi))*(s*s)/6.0;
	de = v_dprE(l, klo)*a + v_dprE(l, khi)*b + ((a*a*a-a)*d_dprE(l, klo)+(b*b*b-b)*d_dprE(l, khi))*(s*s)/6.0;
	 h = v_prfH(l, klo)*a + v_prfH(l, khi)*b + ((a*a*a-a)*d_prfH(l, klo)+(b*b*b-b)*d_prfH(l, khi))*(s*s)/6.0;
	dh = v_dprH(l, klo)*a + v_dprH(l, khi)*b + ((a*a*a-a)*d_dprH(l, klo)+(b*b*b-b)*d_dprH(l, khi))*(s*s)/6.0;
	return;
}

/* initialize discrete profiles 
   on radial interval [ir < wg.hx(.) < re] with nr steps */
void FIMode::discretize(double ri, double re, int nr)
{
	disc_prof = 0;
	n_rad = Ivector(wg.nx+2);	
	Dvector ra(wg.nx+2);
	Dvector rb(wg.nx+2);
	Dvector ds(wg.nx+2);
	Ivector ns(wg.nx+2);

	if(ri > R+wg.hx(0)) ri = (wg.hx(0)+R)-wg.lambda; 
	if(ri < 0.0) ri = 0.0;
	if(re < R+wg.hx(wg.nx)) re = (wg.hx(wg.nx)+R)+wg.lambda;
	
	d_ri = ri;
	d_re = re;

	double dr = (re-ri)/((double) nr);

	for(int l=0; l<=wg.nx+1; ++l)
	{
		if(l == 0) ra(l) = ri; 
		else       ra(l) = wg.hx(l-1)+R;
		if(l == wg.nx+1) rb(l) = re; 
		else             rb(l) = wg.hx(l)+R;
		ns(l) = (int)(ceil((rb(l)-ra(l))/dr));
		if(ns(l) < 2) ns(l) = 2;
		ds(l) = (rb(l)-ra(l))/((double) ns(l)); 
		n_rad(l) = ns(l)+3;
	}		
	int maxn = n_rad.max();
	v_rad = Dmatrix(wg.nx+2, maxn);
	v_prfE = Dmatrix(wg.nx+2, maxn);
	v_dprE = Dmatrix(wg.nx+2, maxn);
	d_prfE = Dmatrix(wg.nx+2, maxn);
	d_dprE = Dmatrix(wg.nx+2, maxn);
	v_prfH = Dmatrix(wg.nx+2, maxn);
	v_dprH = Dmatrix(wg.nx+2, maxn);
	d_prfH = Dmatrix(wg.nx+2, maxn);
	d_dprH = Dmatrix(wg.nx+2, maxn);
	for(int l=0; l<=wg.nx+1; ++l)
	{
		int j=0;
		v_rad(l, j) = ra(l); ++j;
		v_rad(l, j) = ra(l)+ds(l)*0.02; ++j;
		v_rad(l, j) = ra(l)+ds(l); ++j;
		for(int m=2; m<=ns(l)-1; ++m)
		{
			v_rad(l, j) = v_rad(l, j-1)+ds(l); ++j;
		}
		v_rad(l, j) = rb(l)-ds(l)*0.02; ++j;
		v_rad(l, j) = rb(l); ++j;
		if(j != n_rad(l)) fimerror("discretize: confused");
	}
	for(int l=0; l<=wg.nx+1; ++l)
	{
		for(int m=0; m<=n_rad(l)-1; ++m)
		{
			double tdr = 0.0;
			if(m == 0)          tdr =  HDIST;
			if(m == n_rad(l)-1) tdr = -HDIST;
			double e, de, h, dh;
			profile(v_rad(l, m)+tdr, e, de, h, dh);
			v_prfE(l, m) = e;
			v_dprE(l, m) = de;
			v_prfH(l, m) = h;
			v_dprH(l, m) = dh;
		}
	}
	for(int l=0; l<=wg.nx+1; ++l)
	{
		int n;
		Dvector x;
		Dvector y;
		Dvector y2;
		n = n_rad(l);
		Dvector u(n);
		y2 = Dvector(maxn);

		double sig;
		double p;
		x = v_rad.row(l).subvector(n, 0);

		y = v_prfE.row(l).subvector(n, 0);
		y2.init(0.0);
		u(0) = 0.0;
		for(int i=1; i<=n-2; ++i) 
		{
			sig = (x(i)-x(i-1))/(x(i+1)-x(i-1));
			p = sig*y2(i-1)+2.0;
			y2(i) = (sig-1.0)/p;
			u(i) = (y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1))/(x(i)-x(i-1));
			u(i) = (6.0*u(i)/(x(i+1)-x(i-1))-sig*u(i-1))/p;
		}
		y2(n-1) = 0.0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_prfE.row(y2, l);

		y = v_dprE.row(l).subvector(n, 0);
		y2.init(0.0);
		u(0) = 0.0;
		for(int i=1; i<=n-2; ++i) 
		{
			sig = (x(i)-x(i-1))/(x(i+1)-x(i-1));
			p = sig*y2(i-1)+2.0;
			y2(i) = (sig-1.0)/p;
			u(i) = (y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1))/(x(i)-x(i-1));
			u(i) = (6.0*u(i)/(x(i+1)-x(i-1))-sig*u(i-1))/p;
		}
		y2(n-1) = 0.0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_dprE.row(y2, l);
	}
	for(int l=0; l<=wg.nx+1; ++l)
	{
		int n;
		Dvector x;
		Dvector y;
		Dvector y2;
		n = n_rad(l);
		Dvector u(n);
		y2 = Dvector(maxn);

		double sig;
		double p;
		x = v_rad.row(l).subvector(n, 0);

		y = v_prfH.row(l).subvector(n, 0);
		y2.init(0.0);
		u(0) = 0.0;
		for(int i=1; i<=n-2; ++i) 
		{
			sig = (x(i)-x(i-1))/(x(i+1)-x(i-1));
			p = sig*y2(i-1)+2.0;
			y2(i) = (sig-1.0)/p;
			u(i) = (y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1))/(x(i)-x(i-1));
			u(i) = (6.0*u(i)/(x(i+1)-x(i-1))-sig*u(i-1))/p;
		}
		y2(n-1) = 0.0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_prfH.row(y2, l);

		y = v_dprH.row(l).subvector(n, 0);
		y2.init(0.0);
		u(0) = 0.0;
		for(int i=1; i<=n-2; ++i) 
		{
			sig = (x(i)-x(i-1))/(x(i+1)-x(i-1));
			p = sig*y2(i-1)+2.0;
			y2(i) = (sig-1.0)/p;
			u(i) = (y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1))/(x(i)-x(i-1));
			u(i) = (6.0*u(i)/(x(i+1)-x(i-1))-sig*u(i-1))/p;
		}
		y2(n-1) = 0.0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_dprH.row(y2, l);
	}
	disc_prof = 1;
	return;
}

/* FIM field */
Complex FIMode::field(Fcomp cp, double x, double y) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dy;
	double cost, sint, r, t;
	double e, de, h, dh;
	Complex efac;
	dx = x-x0;
	dy = y-y0;
	r = sqrt(dx*dx+dy*dy);
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	cost = dx/r; 
	sint = dy/r; 
	t = atan2(dy, dx); if(t<0) t+=2.0*PI;
	profile(r, e, de, h, dh);
	e *= s; de *= s;
	efac = expi(-((double)aol)*t*s);
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor = 0.0;
	if(aol > 0) lor = ((double)aol)*s/r;

	switch(cp)
	{
	case EZ:
		return CC1*e*efac;
		break;
	case HZ:
		return CCmI*h*efac;
		break;
	case ER:
		return CCI*mookq*(-ommu0*lor*h+beta*de)*efac;
		break;
	case HR:
		return CC1*mookq*(-omep0*nq*lor*e+beta*dh)*efac;
		break;
	case ET:
		return CC1*mookq*(-ommu0*dh+beta*lor*e)*efac;
		break;
	case HT:
		return CCI*mookq*(omep0*nq*de-beta*lor*h)*efac;
		break;
	case EX: 
		return mookq*(CCI*(-ommu0*lor*h+beta*de)*cost-(-ommu0*dh+beta*lor*e)*sint)*efac;
		break;
	case EY:
		return mookq*(CCI*(-ommu0*lor*h+beta*de)*sint+(-ommu0*dh+beta*lor*e)*cost)*efac;
		break;
	case HX:
		return mookq*((-omep0*nq*lor*e+beta*dh)*cost-CCI*(omep0*nq*de-beta*lor*h)*sint)*efac;
		break;
	case HY:
		return mookq*((-omep0*nq*lor*e+beta*dh)*sint+CCI*(omep0*nq*de-beta*lor*h)*cost)*efac;
		break;
	default:
		fimerror("field: cp not implemented");
		break;
	}
	return CC0;
}
Complex FIMode::field(Fcomp cp, double r) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	double e, de, h, dh;
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	profile(r, e, de, h, dh);
	e *= s; de *= s;
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor = 0.0;
	if(aol > 0 ) lor = ((double)aol)*s/r;

	switch(cp)
	{
	case EZ:
		return CC1*e;
		break;
	case HZ:
		return CCmI*h;
		break;
	case ER:
		return CCI*mookq*(-ommu0*lor*h+beta*de);
		break;
	case HR:
		return CC1*mookq*(-omep0*nq*lor*e+beta*dh);
		break;
	case ET:
		return CC1*mookq*(-ommu0*dh+beta*lor*e);
		break;
	case HT:
		return CCI*mookq*(omep0*nq*de-beta*lor*h);
		break;
	default:
		fimerror("field: cp not implemented");
		break;
	}
	return CC0;
}
void FIMode::emfield(double x, double y, 
	             Complex &ex, Complex &ey, Complex &ez, 
	             Complex &hx, Complex &hy, Complex &hz) const
{
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dy;
	double cost, sint, r, t;
	double pe, de, ph, dh;
	Complex efac;
	dx = x-x0;
	dy = y-y0;
	r = sqrt(dx*dx+dy*dy);
	ex = CC0; ey = CC0; ez = CC0; 
	hx = CC0; hy = CC0; hz = CC0; 
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	cost = dx/r; 
	sint = dy/r; 
	t = atan2(dy, dx); if(t<0) t+=2.0*PI;
	profile(r, pe, de, ph, dh);
	pe *= s; de *= s;
	efac = expi(-((double)aol)*t*s);
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor = 0.0;
	if(aol > 0) lor = ((double)aol)*s/r;
	ex = mookq*(CCI*(-ommu0*lor*ph+beta*de)*cost-(-ommu0*dh+beta*lor*pe)*sint)*efac;
	ey = mookq*(CCI*(-ommu0*lor*ph+beta*de)*sint+(-ommu0*dh+beta*lor*pe)*cost)*efac;
	ez = pe*efac;
	hx = mookq*((-omep0*nq*lor*pe+beta*dh)*cost-CCI*(omep0*nq*de-beta*lor*ph)*sint)*efac;
	hy = mookq*((-omep0*nq*lor*pe+beta*dh)*sint+CCI*(omep0*nq*de-beta*lor*ph)*cost)*efac;
	hz = CCmI*ph*efac;
	return;
}
void FIMode::emfield_r(double x, double y, 
	               Complex &er, Complex &et, Complex &ez, 
	               Complex &hr, Complex &ht, Complex &hz) const
{
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dy;
	double r, t;
	double pe, de, ph, dh;
	Complex efac;
	dx = x-x0;
	dy = y-y0;
	r = sqrt(dx*dx+dy*dy);
	er = CC0; et = CC0; ez = CC0; 
	hr = CC0; ht = CC0; hz = CC0; 
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	t = atan2(dy, dx); if(t<0) t+=2.0*PI;
	profile(r, pe, de, ph, dh);
	pe *= s; de *= s;
	efac = expi(-((double)aol)*t*s);
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor = 0.0;
	if(aol > 0) lor = ((double)aol)*s/r;

	er = CCI*mookq*(-ommu0*lor*ph+beta*de)*efac;
	et = CC1*mookq*(-ommu0*dh+beta*lor*pe)*efac;
	ez = CC1*pe*efac;
	hr = CC1*mookq*(-omep0*nq*lor*pe+beta*dh)*efac;
	ht = CCI*mookq*(omep0*nq*de-beta*lor*ph)*efac;
	hz = CCmI*ph*efac;
	return;
}

/* directional interference: mix with a*(reversed mode) */
Complex FIMode::field_ri(Fcomp cp, double x, double y, Complex a) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dy;
	double cost, sint, r, t;
	double e, de, h, dh;
	double er, der, hr, dhr;

	Complex efac;
	dx = x-x0;
	dy = y-y0;
	r = sqrt(dx*dx+dy*dy);
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	cost = dx/r; 
	sint = dy/r; 
	t = atan2(dy, dx); if(t<0) t+=2.0*PI;

	profile(r, e, de, h, dh);
	e = e*s;
	de = de*s;
	er = -e;
	der = -de;
	hr = h; 
	dhr = dh;

	efac = expi(-((double)aol)*s*t);
	double sm = -s;
	Complex efacm = expi(-((double)aol)*sm*t);
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor  = 0.0;
	double lorm = 0.0;
	if(aol > 0) { lor = ((double)aol)*s/r; lorm = ((double)aol)*sm/r; }

	switch(cp)
	{
	case EZ:
		return CC1*(e*efac + a*er*efacm);
		break;
	case HZ:
		return CCmI*(h*efac + a*hr*efacm);
		break;
	case ER:
		return CCI*mookq*((-ommu0*lor*h+beta*de)*efac + a*(-ommu0*lorm*hr+beta*der)*efacm);
		break;
	case HR:
		return CC1*mookq*((-omep0*nq*lor*e+beta*dh)*efac + a*(-omep0*nq*lorm*er+beta*dhr)*efacm);
		break;
	case ET:
		return CC1*mookq*((-ommu0*dh+beta*lor*e)*efac + a*(-ommu0*dhr+beta*lorm*er)*efacm);
		break;
	case HT:
		return CCI*mookq*((omep0*nq*de-beta*lor*h)*efac + a*(omep0*nq*der-beta*lorm*hr)*efacm);
		break;
	case EX: 
		return mookq*(CCI*(-ommu0*lor*h+beta*de)*cost-(-ommu0*dh+beta*lor*e)*sint)*efac
		  + a*(mookq*(CCI*(-ommu0*lorm*hr+beta*der)*cost-(-ommu0*dhr+beta*lorm*er)*sint)*efacm);
		break;
	case EY:
		return mookq*(CCI*(-ommu0*lor*h+beta*de)*sint+(-ommu0*dh+beta*lor*e)*cost)*efac
		  + a*(mookq*(CCI*(-ommu0*lorm*hr+beta*der)*sint+(-ommu0*dhr+beta*lorm*er)*cost)*efacm);
		break;
	case HX:
		return mookq*((-omep0*nq*lor*e+beta*dh)*cost-CCI*(omep0*nq*de-beta*lor*h)*sint)*efac
		  + a*(mookq*((-omep0*nq*lorm*er+beta*dhr)*cost-CCI*(omep0*nq*der-beta*lorm*hr)*sint)*efacm);
		break;
	case HY:
		return mookq*((-omep0*nq*lor*e+beta*dh)*sint+CCI*(omep0*nq*de-beta*lor*h)*cost)*efac
		  + a*(mookq*((-omep0*nq*lorm*er+beta*dhr)*sint+CCI*(omep0*nq*der-beta*lorm*hr)*cost)*efacm);
		break;
	default:
		fimerror("field: cp not implemented");
		break;
	}
	return CC0;
}
void FIMode::emfield_ri(double x, double y, Complex a, 
	               Complex &er, Complex &et, Complex &ez, 
	               Complex &hr, Complex &ht, Complex &hz) const
{
	if(a.abs() < 1.0e-10)
	{
		emfield_r(x, y, er, et, ez, hr, ht, hz);
		return;
	}
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dy;
	double r, t;
	double pe, de, ph, dh;
	double per, der, phr, dhr;
	Complex efac;
	dx = x-x0;
	dy = y-y0;
	r = sqrt(dx*dx+dy*dy);
	er = CC0; et = CC0; ez = CC0; 
	hr = CC0; ht = CC0; hz = CC0; 
	if(r < FIMp_R_TINY*wg.lambda) r = FIMp_R_TINY*wg.lambda;
	t = atan2(dy, dx); if(t<0) t+=2.0*PI;
	profile(r, pe, de, ph, dh);
	pe = pe*s;
	de = de*s;
	per = -pe;
	der = -de;
	phr = ph; 
	dhr = dh;
	efac = expi(-((double)aol)*s*t);
	double sm = -s;
	Complex efacm = expi(-((double)aol)*sm*t);
	double nq = wg.eps(r-R);
	double mookq = -1.0/(k0*k0*(nq-neff*neff));
	double lor  = 0.0;
	double lorm = 0.0;
	if(aol > 0) { lor = ((double)aol)*s/r; lorm = ((double)aol)*sm/r; }

	er = CCI*(mookq*(-ommu0*lor*ph+beta*de)*efac + a*(mookq*(-ommu0*lorm*phr+beta*der)*efacm));
	et = mookq*(-ommu0*dh+beta*lor*pe)*efac + a*(mookq*(-ommu0*dhr+beta*lorm*per)*efacm);
	ez = pe*efac + a*(per*efacm);
	hr = mookq*(-omep0*nq*lor*pe+beta*dh)*efac + a*(mookq*(-omep0*nq*lorm*per+beta*dhr)*efacm);
	ht = CCI*(mookq*(omep0*nq*de-beta*lor*ph)*efac + a*(mookq*(omep0*nq*der-beta*lorm*phr)*efacm));
	hz = CCmI*(ph*efac + a*phr*efacm);
	return;
}

/* ... values on a rectangular mesh */
Cmatrix FIMode::field(Fcomp cp, Rect disp, int npx, int npy) const
{
	if(npx <= 1 || npy <= 1) fimerror("discretized field: npx, npy"); 
	Cmatrix fld(npx, npy);
	double dx = (disp.x1-disp.x0)/((double)(npx-1));
	double dy = (disp.y1-disp.y0)/((double)(npy-1));
	double x, y;

	y = disp.y0;
	for(int j=0; j<=npy-1; ++j)
	{
		x = disp.x0;
		for(int i=0; i<=npx-1; ++i) 
		{
			fld(i, j) = field(cp, x, y);
			x += dx;
		}
		y += dy;
	}
	return fld;
}

/* set the mode id string, for mode count o */
void FIMode::classify(int o)
{
	ord = o;
	int l = abs(aol);
	int i;
	i = 0;
	if(aol == 0)
	{
		ids[i++] = 'T';
		if(pid == 'e') ids[i++] = 'E';
		else           ids[i++] = 'M';
	}
	else
	{
		ids[i++] = 'O';
		ids[i++] = 'A';
		ids[i++] = 'M';
	}
	ids[i++] = '(';
        if(dir == BACK && l >= 1) ids[i++] = '-';
        if(l >=1000)     ids[i++] = dig1000(l);
        if(l >=100)      ids[i++] = dig100(l);
        if(l >=10)       ids[i++] = dig10(l);
        ids[i++] = dig1(l);
	ids[i++] = ',';
        if(ord >=1000)  ids[i++] = dig1000(ord);
        if(ord >=100)   ids[i++] = dig100(ord);
        if(ord >=10)    ids[i++] = dig10(ord);
        ids[i++] = dig1(ord);
	ids[i++] = ')';
        while(i<=19) ids[i++] = 0;

	i = 0;
	if(aol == 0)
	{
		sids[i++] = 'T';
		if(pid == 'e') sids[i++] = 'E';
		else           sids[i++] = 'M';
	}
	else
	{
		sids[i++] = 'O';
		sids[i++] = 'A';
		sids[i++] = 'M';
	}
	sids[i++] = '_';
        if(dir == BACK && l >= 1) sids[i++] = '-';
        if(l >=1000)     sids[i++] = dig1000(l);
        if(l >=100)      sids[i++] = dig100(l);
        if(l >=10)       sids[i++] = dig10(l);
        sids[i++] = dig1(l);
	sids[i++] = '_';
        if(ord >=1000)  sids[i++] = dig1000(ord);
        if(ord >=100)   sids[i++] = dig100(ord);
        if(ord >=10)    sids[i++] = dig10(ord);
        sids[i++] = dig1(ord);
        while(i<=19) sids[i++] = 0;
        return;
}

/* local axial Poynting vector component */
double FIMode::pfSz(double r) const
{ 
	Complex er, et, hr, ht;
	er = field(ER, r);
	et = field(ET, r);
	hr = field(HR, r);
	ht = field(HT, r);
	return 0.5*(er.conj()*ht-et.conj()*hr).re;
}

/* power transported in axial direction */
double FIMode::power() const
{
        static double gipx[] = {0.1488743389, 0.4333953941, 0.6794095682, 0.8650633666, 0.9739065285};
        static double gipw[] = {0.2955242247, 0.2692667193, 0.2190863625, 0.1494513491, 0.0666713443};

	double p = 0.0;
	double ds0 = wg.lambda/wg.n.max()/100.0;
	for(int l=0; l<=wg.nx+1; ++l)
	{
		double ra, rb;
		if(l == 0) ra = cRi; 
		else       ra = wg.hx(l-1)+R;
		if(l == wg.nx+1) rb = cRe; 
		else             rb = wg.hx(l)+R;

		int ns = (int) ceil((rb-ra)/ds0);
		if(ns < 10) ns = 10;
		double ds = (rb-ra)/((double) ns);
		double a, b;
		a = ra; 
		for(int i=0; i<=ns-1; ++i)
		{
			b = a+ds;
			double xm = 0.5*(b+a);
			double xr = 0.5*(b-a);
			double s = 0.0;
			for(int j=0; j<=4; j++) 
			{
				double dx = xr*gipx[j];
				double rx;
				rx = xm+dx;
				double fp = rx*pfSz(rx);
				rx = xm-dx;
				double fm = rx*pfSz(rx);
				s += gipw[j]*(fp+fm);
			}
			p += s*xr;
			a = b;
		}
	}		
	p *= 2.0*PI;
	return p;
}

/* normalize mode to unit power */ 
void FIMode::normalize()
{
	double p;
	double ri = 0.0;
	double re = wg.hx(wg.nx)+R;
	p = fabs(pfSz(re));
	double lwl;
	lwl = wg.lambda/wg.n(wg.nx+1);
	re += lwl;
	int i = 0;
	while(i<=1000 && fabs(pfSz(re)) > p/1.0e+6) { re += lwl; ++i; }
	if(i >= 1000) fimerror("normalize: convergence");
	cRi = ri;
	cRe = re;

	p = power();
	p = 1.0/sqrt(fabs(p));
	cBe *= p;	
	cDe *= p;	
	fillc = 2; tr_error = fabs(travres());
	return;
}


/* test all interface conditions, 0: reject, 1: acceptable */
double FIMsP_IFCerr = 1.0e-10;
int FIMode::check()
{
	double err = 0.0;
	for(int i=0; i<=wg.nx; ++i)	
	{
		double r = wg.hx(i)+R;
		double dr = 1.0e-10;
		Complex f0, f1;
		f0 = field(HR, r-dr); f1 = field(HR, r+dr); err += (f0-f1).sqabs();
		f0 = field(HT, r-dr); f1 = field(HT, r+dr); err += (f0-f1).sqabs();
		f0 = field(HZ, r-dr); f1 = field(HZ, r+dr); err += (f0-f1).sqabs();
	}
// fprintf(stderr, "\n --- err: %g - ", err);
	if(err >= FIMsP_IFCerr) return 0;
	return 1;
}

/* local product of field components, conjugate of first field */
Complex FIMode::lfp(Fcomp cpc, Fcomp cpn, double r) const
{ 
	return field(cpc, r).conj()*field(cpn, r);
}

/* numerically integrate product of field components (conjugate of first field)
   over full circle x radial layer l */
Complex FIMode::numquad(Fcomp cpc, Fcomp cpn, int l) const
{
        static double gipx[] = {0.1488743389, 0.4333953941, 0.6794095682, 0.8650633666, 0.9739065285};
        static double gipw[] = {0.2955242247, 0.2692667193, 0.2190863625, 0.1494513491, 0.0666713443};

	if(l<0 || l>wg.nx+1) return CC0;
	Complex p = CC0;
	double ds0 = wg.lambda/wg.n.max()/100.0;
	double ra, rb;
	if(l == 0) ra = cRi; 
	else       ra = wg.hx(l-1)+R;
	if(l == wg.nx+1) rb = cRe; 
	else             rb = wg.hx(l)+R;

	int ns = (int) ceil((rb-ra)/ds0);
	if(ns < 10) ns = 10;
	double ds = (rb-ra)/((double) ns);
	double a, b;
	a = ra; 
	for(int i=0; i<=ns-1; ++i)
	{
		b = a+ds;
		double xm = 0.5*(b+a);
		double xr = 0.5*(b-a);
		Complex s = CC0;
		for(int j=0; j<=4; j++) 
		{
			double dx = xr*gipx[j];
			double rx;
			rx = xm+dx;
			Complex fp = rx*lfp(cpc, cpn, rx);
			rx = xm-dx;
			Complex fm = rx*lfp(cpc, cpn, rx);
			s += gipw[j]*(fp+fm);
		}
		p += s*xr;
		a = b;
	}
	p *= 2.0*PI;
	return p;
}

/* determine normalized effective permittivity B */
void FIMode::setnpcB()
{
	npcB = 1.0;
	double nmin = wg.n(wg.nx+1);
	double nmax = 0.0;
	for(int l=0; l<=wg.nx; ++l) if(wg.n(l) > nmax) nmax = wg.n(l);
	if(nmin >= nmax) return;
	npcB = (neff*neff-nmin*nmin)/(nmax*nmax-nmin*nmin);
	return;
}

/* translate the FIM: shift the origin by dx, dy */
void FIMode::translate(double dx, double dy)
{
	x0 += dx; 
	y0 += dy; 
	for(int o=0; o<=fib.no-1; ++o)	
	{
		switch(fib(o).type) 
		{
		case RING: case DISK:
			fib(o).xo += dx;
			fib(o).zo += dy;
			break;
		default:
			fimerror("translate: fib.type");
			break;
		}
	}
	return;
}

/* reverse the FIM: clockwise propagation <-> anticlockwise propagation */
void FIMode::reverse()
{
	if(dir == FORW) dir = BACK;
	else            dir = FORW;
	classify(ord);
	return;
}

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

/* FIM field profile plots, image of component cp (static, t=0)
   cp: one of EX, EY, EZ, HX, HY, HZ, ER, ET, HR, HT
   foa: MOD, SQR, REP, IMP, ORG = REP
   disp: x resp. y extensions of the plot window
   npx, npy: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void FIMode::plot(Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
                  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", 
	                disp.x0, disp.x1, disp.y0, disp.y1, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "FIM profile");
	mlout_gengeoxy(dat, fib, disp);
	mlout_meshxy(dat, disp, npx, npy);
	Cmatrix cfld = field(cp, disp, npx, npy);
	Dmatrix fld = realize(cfld, foa); 
	mlout_fld(dat, npx, npy, cp, fld);
	name[8] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage('y', 'x', 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;
}

/* export full field data into a viewer m-file
   disp: x resp. y extensions of the plot window
   npx, npy: number of plot points in the x and z directions
   ext0, ext1: filename id characters;
   field animations ignore the attenuation */
void FIMode::viewer(Rect disp, int npx, int npy, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "prf__A.m";
	name[3] = ext0;
	name[4] = ext1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   disp.x0, npx, disp.x1, disp.y0, npy, disp.y1, name);
	dat = fopen(name, "w+");
	name[6] = 0; 
	mlout_viewertopxy(dat, name, wg.lambda);
	mlout_gengeoxy(dat, fib, disp);
	mlout_meshxy(dat, disp, npx, npy);
	Cmatrix f;
	Fcomp cp;

	cp = EX;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);
	cp = EY;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);
	cp = EZ;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);
	cp = HX;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);
	cp = HY;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);
	cp = HZ;
	f = field(cp, disp, npx, npy); 
	mlout_fld(dat, npx, npy, cp, f.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, npx, npy, cp, f.im());
	mlout_fldtoim(dat, cp);

	Dmatrix n(npx, npy);
	double dx = (disp.x1-disp.x0)/((double)(npx-1));
	double dy = (disp.y1-disp.y0)/((double)(npy-1));
	double x, y;
	y = disp.y0;
	for(int j=0; j<=npy-1; ++j)
	{
		x = disp.x0;
		for(int i=0; i<=npx-1; ++i) { n(i, j) = fib.n(x, y); x += dx; }
		y += dy;
	}
	mlout_n(dat, npx, npy, n);
	mlout_fldviewerxy(dat, name);
	fclose(dat);
	return;
}

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

/* initialize */
FIModeArray::FIModeArray()
            : num(0)
{
	mvec = NULL;
}

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

/* copy constructor */
FIModeArray::FIModeArray(const FIModeArray& ma)
	    : num(ma.num)
{
	mvec = new FIMode[num];
	FIMode* ap = ma.mvec;
	FIMode* mp =    mvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}

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

}

/* delete all Mode entries */
void FIModeArray::clear()
{
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0;

}

/* subscripting */
FIMode& FIModeArray::operator() (int i)  
{ 
	if(i<0 || i>=num) fimerror("FIModeArray: () out of range");
	return mvec[i]; 
}
FIMode FIModeArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) fimerror("FIModeArray: () out of range");
	return mvec[i]; 
}

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

/* delete a mode entry */
void FIModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		fimerror("FIModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	FIMode *avec;
	avec = new FIMode[num-1];
	FIMode* ap = avec;
	FIMode* mp = mvec;
	for(int j=0; j<=i-1; ++j) *ap++ = *mp++;
	mp++;
	for(int j=i+1; j<=num-1; ++j) *ap++ = *mp++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	--num;
	return;
}

/* add an entire FIModeArray nma */ 
void FIModeArray::merge(FIModeArray& ma)
{
	FIMode *avec;
	avec = new FIMode[num+ma.num];
	FIMode* ap = avec;
	FIMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	FIMode* np = ma.mvec;
	for(int i=0; i<=ma.num-1; ++i) *ap++ = *np++;
	if(mvec != NULL) delete[] mvec;
	mvec = avec;
	num += ma.num;
	return;
}

int compare(const FIMode& a, const FIMode& b)
{	
	if(a.neff > b.neff) return -1;
	return 1;
}

/* sort the array */
void FIModeArray::sort()
{
	int j, k, mini;
	FIMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		mini = j; 
		for(k=j+1; k<=num-1; ++k)
		{
			if(compare(mvec[k], mvec[mini]) == -1)
			{
				mini = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[mini]; mvec[mini] = t;
	}
	return;
}

/* prune the array: remove modes with identical orders l, m */
void FIModeArray::prune()
{
        int done = 0;
	while(done == 0 && num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=num-1; ++j)
		{
			double nj = mvec[j].neff;
			int mj = mvec[j].aol;
			for(int k=0; k<=num-1; ++k)
			{
				double nk = mvec[k].neff;
				int mk = mvec[k].aol;
				if(j!=k && done == 1 && mj==mk && fabs(nj-nk) < FIMsP_NeffSep)
				{
					if(mvec[j].tr_error <= mvec[k].tr_error) rmm = k;
					else                                     rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) remove(rmm);
	}
	return;
}

/* ------------------------------------------------------------------------
   FIM mode analysis procedures
   ------------------------------------------------------------------------ */

/* --- solver procedures ---------------------------------------------------*/ 

/* small radius/wavelength, to avoid singluarity at the origin */
double FIMp_R_TINY = 1.0e-6;

/* for the fiber with radial layering w, radius r, centered at (xc, yc),  
   identify guided FIMs at vacuum wavelength w.lambda, of angular order l  
   | l==0: TE (p='e') or TM (p='m') 
   quiet = 0, 1, 2: all, less, no  log output  */
int fimsolveproc(Waveguide w, double r, double xc, double yc, int l, char p, FIModeArray& fima, int quiet)
{
	fima.clear();
if(quiet < 1)
{	
	if(l == 0)
	{
		if(p == 'e') fprintf(stderr, "<l=%d, TE>\n", l);
		else         fprintf(stderr, "<l=%d, TM>\n", l);
	}
	else fprintf(stderr, "<l=%d, HYB>\n", l);
}
	FIMode m;
	m.R = r;
	m.x0 = xc;
	m.y0 = yc;
	m.wg = w;
	m.fib = OvlStruct(w.n(w.nx+1), w.lambda);
	m.fib.clear();
        for(int li=w.nx; li>=0; --li)
        {
                if(fabs(w.n(li+1)-w.n(li)) > COMPTOL_N)
                {
                        Overlay o(DISK, w.n(li), r+w.hx(li), xc, yc);
                        m.fib.place(o);
                }
        }
	m.aol = l;
	m.pid = p;
	m.k0 = val_k0(w.lambda);
	m.omep0 = 1.0/val_invomep0(w.lambda);
	m.ommu0 = 1.0/val_invommu0(w.lambda);
	m.dir = FORW;
	m.cA = Dvector(w.nx+2);
	m.cB = Dvector(w.nx+2);
	m.cC = Dvector(w.nx+2);
	m.cD = Dvector(w.nx+2);
	m.kappa = Dvector(w.nx+2);
	m.cA.init(0.0);
	m.cB.init(0.0);
	m.cC.init(0.0);
	m.cD.init(0.0);
	m.kappa.init(1.0);
	m.ids[0] = 0;
	m.sids[0] = 0;
	m.ord = 0;

	double nmax = 0.0;
	for(int l=0; l<=w.nx; ++l) if(w.n(l) > nmax) nmax = w.n(l);
	double nmin = w.n(w.nx+1);
	if(nmin+FIMsP_MinNDiff >= nmax-FIMsP_MinNDiff)
	{
if(quiet < 1)
{	
		fprintf(stderr, "neff in [%g, %g] ?\n", nmin, nmax);
}
		return fima.num;
	}
		
	Dvector rilv;
	rilv.strip();
	rilv.append(nmin);
	rilv.append(nmax);
	for(int l=0; l<=w.nx; ++l)
	{
		double nx = w.n(l);
		bool done = false;
		for(int i=0; i<=rilv.nel-1; ++i)
		{
			if(fabs(rilv(i) - nx) <= FIMsP_MinNDiff) done = true;
		}
		if(done == false) rilv.append(nx);
	}
	rilv.sort(1);
	double n0, n1, no, n, nf, s, s0;
	double tro, tr;
	s0 = (nmax-nmin)/1000;
	for(int r=rilv.nel-2; r>=0; --r)
	{
		n0  = rilv(r)+FIMsP_MinNDiff;
		n1  = rilv(r+1)-FIMsP_MinNDiff;
		s = s0;
		if((n1-n0)/s < 100.0) s = (n1-n0)/100.0;
		no = n1;
		tro = m.travres(no, 0);
		n = no-s;
if(quiet < 1) fprintf(stderr, "    neff[%g:%g]: ", n0, n1); 
		while(n >= n0)
		{
			tr = m.travres(n, 0);
			if(rootin(tr, tro))
			{
				nf = m.converge(n, no);
				m.tr_error = fabs(m.travres(nf, 1));
				m.setnpcB();
				m.normalize();
				if(m.check() == 1)
				{
					m.classify(fima.num+1);
					fima.add(m);
if(quiet < 1) fprintf(stderr, "%s ", m.ids); 
				}
			}
			no = n; tro = tr;
			n -= s;
		}
if(quiet < 1) fprintf(stderr, "\n"); 
	}
	fima.prune();
	fima.sort();
	return fima.num;
/*
	m.tr_error = fabs(m.travres(nf, 1));
fprintf(stderr, "cA:\n"); m.cA.write(stderr);
fprintf(stderr, "cB:\n"); m.cB.write(stderr);
fprintf(stderr, "cC:\n"); m.cC.write(stderr);
fprintf(stderr, "cD:\n"); m.cD.write(stderr);

for(double r=0.01; r<=10; r+=0.001)
{
	double e, de, h, dh;
	m.profile(r, e, de, h, dh);
	apptoxyf("e", r, e);
	apptoxyf("de", r, de);
	apptoxyf("h", r, h);
	apptoxyf("dh", r, dh);
}
*/
}

/* ... find modes of all angular orders */
int fimsolve(Waveguide w, double r, double xc, double zc, FIModeArray& fima, int quiet)
{
	return fimsolve(w, r, xc, zc, 0, 100, fima, quiet);
}
/* ... find modes with angular orders between and including lmin and lmax */
int fimsolve(Waveguide w, double r, double xc, double zc, int lmin, int lmax, FIModeArray& fima, int quiet)
{	
	int l;
	lmin = abs(lmin);
	lmax = abs(lmax);
	if(lmax < lmin) { l = lmax; lmax = lmin; lmin = l; }
	if(quiet < 2)
	{
fprintf(stderr, "\n------------- FIMsolve --------------------------------------------- '19 ---\n");
fprintf(stderr, "R: %g   (x0, z0): (%g, %g)   lambda: %g\n", r, xc, zc, w.lambda);
fprintf(stderr,   "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", w.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=w.nx; ++j) fprintf(stderr, "%6.10g  ", w.hx(j));
fprintf(stderr, "\n");
if(w.special) fprintf(stderr, "   eps:\n");
else          fprintf(stderr, "     n:\n");
for(int i=0; i<=w.nx+1; ++i) fprintf(stderr, "%6.10g  ", w.n(i));
if(w.special) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "l: [%d, %d]\n", lmin, lmax);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	fima.clear();
	FIModeArray tfima;
	for(l=lmin; l<=lmax; ++l)
	{
		if(l == 0)
		{
			tfima.clear();
			fimsolveproc(w, r, xc, zc, l, 'e', tfima, quiet);
			fima.merge(tfima);
			tfima.clear();
			fimsolveproc(w, r, xc, zc, l, 'm', tfima, quiet);
			fima.merge(tfima);
		}
		else 
		{
			tfima.clear();
			fimsolveproc(w, r, xc, zc, l, 'v', tfima, quiet);
			if(tfima.num <= 0 && l >= 1) l=lmax;
			else fima.merge(tfima);
		}
	}
	fima.sort();
if(quiet < 1)
{
	fprintf(stderr, "\n----------------------------------------------------------------------------\n");
}
if(quiet < 2)
{
	if(fima.num<=0) fprintf(stderr, "No guided modes found.\n");
	else
	for(int j=0; j<=fima.num-1; ++j)
		fprintf(stderr, "%s: n_eff = %.8g, beta = %.8g, B = %.8g\n", fima(j).ids, fima(j).neff, fima(j).beta, fima(j).npcB);
}
	return fima.num;
}