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

/*
 * bend.cpp
 * 2-D modes of bent slab waveguides 
 */

#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"bend.h"

/* error message */
void benderror(const char *msg)
{
	fprintf(stderr, "\nbend: %s.\n", msg);
	return;
}

/* Cylinder functions & derivatives, Complex order o, real argument a, 
   wrapper functions: indicate computation failure */
int BDMCylEvalError = 0;

#ifndef BF_CIRC
Complex BfJ(Complex o, double a)
{
	int info = 0;
	BDMCylEvalError = 0;
	Complex c = bessel_J(o, Complex(a, 0.0), info);
	if(info != 0) 
	{
		BDMCylEvalError = 1;
		fprintf(stderr, "J_(%g+i%g)(%g) -> info: %d, failed.\n", o.re, o.im, a, info); 
	}
	return c;
}
Complex BfY(Complex o, double a)
{
	int info = 0;
	BDMCylEvalError = 0;
	Complex c = bessel_Y(o, Complex(a, 0.0), info);
	if(info != 0) 
	{
		BDMCylEvalError = 1;
		fprintf(stderr, "Y_(%g+i%g)(%g) -> info: %d, failed.\n", o.re, o.im, a, info); 
	}
	return c;
}
Complex dBfJ(Complex o, double a)
{
	int info = 0;
	BDMCylEvalError = 0;
	Complex c = d_bessel_J(o, Complex(a, 0.0), info);
	if(info != 0) 
	{
		BDMCylEvalError = 1;
		fprintf(stderr, "dJ_(%g+i%g)(%g) -> info: %d, failed.\n", o.re, o.im, a, info); 
	}
	return c;
}
Complex dBfY(Complex o, double a)
{
	int info = 0;
	BDMCylEvalError = 0;
	Complex c = d_bessel_Y(o, Complex(a, 0.0), info);
	if(info != 0) 
	{
		BDMCylEvalError = 1;
		fprintf(stderr, "dY_(%g+i%g)(%g) -> info: %d, failed.\n", o.re, o.im, a, info); 
	}
	return c;
}
#else 
Complex BfJ(Complex o, double a)
{
	// if(o.re < 0.0) { o.neg(); Complex pio = o*PI; return cos(pio)*BfJ(o, a)-sin(pio)*BfY(o, a); }
	int nz, ierr;
	BDMCylEvalError = 0;
	Complex c = BesselJ(o, Complex(a, 0.0), 1, &nz, &ierr);
	if(ierr != 0) 
	{	
		BDMCylEvalError = 1;
		fprintf(stderr, "J_(%g+i%g)(%g) -> nz: %d ierr: %d, failed.\n", o.re, o.im, a, nz, ierr); 
	}
	return c;
}
Complex BfY(Complex o, double a)
{
	// if(o.re < 0.0) { o.neg(); Complex pio = o*PI; return sin(pio)*BfJ(o, a)+cos(pio)*BfY(o, a); }
	int nz, ierr;
	BDMCylEvalError = 0;
	Complex c = BesselY(o, Complex(a, 0.0), 1, &nz, &ierr);
	if(ierr != 0) 
	{	
		BDMCylEvalError = 1;
		fprintf(stderr, "Y_(%g+i%g)(%g) -> nz: %d ierr: %d, failed.\n", o.re, o.im, a, nz, ierr); 
	}
	return c;
}
Complex dBfJ(Complex o, double a)
{
	// if(o.re < 0.0) { o.neg(); Complex pio = o*PI; return cos(pio)*dBfJ(o, a)-sin(pio)*dBfY(o, a); }
	int nz, ierr;
	BDMCylEvalError = 0;
	Complex c = DBesselJ(o, Complex(a, 0.0), 1, &nz, &ierr);
	if(ierr != 0) 
	{	
		BDMCylEvalError = 1;
		fprintf(stderr, "dJ_(%g+i%g)(%g) -> nz: %d ierr: %d, failed.\n", o.re, o.im, a, nz, ierr); 
	}
	return c;
}
Complex dBfY(Complex o, double a)
{
	// if(o.re < 0.0) { o.neg(); Complex pio = o*PI; return sin(pio)*dBfJ(o, a)+cos(pio)*dBfY(o, a); }
	int nz, ierr;
	BDMCylEvalError = 0;
	Complex c = DBesselY(o, Complex(a, 0.0), 1, &nz, &ierr);
	if(ierr != 0) 
	{	
		BDMCylEvalError = 1;
		fprintf(stderr, "dY_(%g+i%g)(%g) -> nz: %d ierr: %d, failed.\n", o.re, o.im, a, nz, ierr); 
	}
	return c;
}
#endif

/* ---------------------------------------------------------------------- */
// transverse resonance evaluation: limit for invertibility
double BDMsP_TrcLimInvert = 1.0e-16;
// transverse resonance evaluation: error indicator 
int BDMsP_TrcFail = 0;

/* determine suitable initial amplitude */
void BDMode::tr_init(Complex tg)
{
	Complex m0 = tg*R; 
	double r0 = R+wg.hx(wg.nx); 
	double rho0 = k0*wg.n(wg.nx+1)*r0;
	Complex J0 = BfJ(m0, rho0); if(BDMCylEvalError) BDMsP_TrcFail = 1;
	Complex Y0 = BfY(m0, rho0); if(BDMCylEvalError) BDMsP_TrcFail = 1;
	if(J0.sqabs() < 1.0e+140 && Y0.sqabs() < 1.0e+140)
		tr_Ainit = CC1/(J0 - CCI*Y0); 
	else 
		tr_Ainit = 1.0e-70; 
	return;
}

/* evaluate transverse resonance condition, set cA, cB accordingly */
Complex BDMode::travres()
{
	BDMsP_TrcFail = 0;
	double r, rho;
	Complex m = c_m();
	Complex tr;
	cA(wg.nx+1) = tr_Ainit; 
	cB(wg.nx+1) = CCmI*cA(wg.nx+1);
	if(pol == TE) 
	{
		for(int l=wg.nx; l>=0; --l)
		{	
			r = R+wg.hx(l); 
			rho = k0*wg.n(l)*r;
			Complex m_a =  BfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_b =  BfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_c = dBfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_d = dBfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			m_c = m_c*wg.n(l);
			m_d = m_d*wg.n(l);
			rho = k0*wg.n(l+1)*r;
			Complex m_e =  BfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_f =  BfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_g = dBfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			Complex m_h = dBfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
			m_g = m_g*wg.n(l+1);
			m_h = m_h*wg.n(l+1);
			Complex m_D = m_a*m_d-m_b*m_c;
			Complex m_U = m_e*cA(l+1)+m_f*cB(l+1);
			Complex m_V = m_g*cA(l+1)+m_h*cB(l+1);
			m_D.inv();
			if(m_D.isnotOK()) BDMsP_TrcFail = 1;
			cA(l) = m_D*(m_d*m_U-m_b*m_V);
			cB(l) = m_D*(m_a*m_V-m_c*m_U);
			if(cA(l).isnotOK()) BDMsP_TrcFail = 1;
			if(cB(l).isnotOK()) BDMsP_TrcFail = 1;
		}
		tr = cB(0);
		tr_error = tr.abs();
		if(!std::isfinite(tr_error)) BDMsP_TrcFail = 1;
		return tr;
	}
	for(int l=wg.nx; l>=0; --l)
	{	
		r = R+wg.hx(l); 
		rho = k0*wg.n(l)*r;
		Complex m_a =  BfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_b =  BfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_c = dBfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_d = dBfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		m_c = m_c/wg.n(l);
		m_d = m_d/wg.n(l);
		rho = k0*wg.n(l+1)*r;
		Complex m_e =  BfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_f =  BfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_g = dBfJ(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		Complex m_h = dBfY(m, rho); if(BDMCylEvalError) BDMsP_TrcFail = 1;
		m_g = m_g/wg.n(l+1);
		m_h = m_h/wg.n(l+1);
		Complex m_D = m_a*m_d-m_b*m_c;
		Complex m_U = m_e*cA(l+1)+m_f*cB(l+1);
		Complex m_V = m_g*cA(l+1)+m_h*cB(l+1);
		m_D.inv();
		if(m_D.isnotOK()) BDMsP_TrcFail = 1;
		cA(l) = m_D*(m_d*m_U-m_b*m_V);
		cB(l) = m_D*(m_a*m_V-m_c*m_U);
		if(cA(l).isnotOK()) BDMsP_TrcFail = 1;
		if(cB(l).isnotOK()) BDMsP_TrcFail = 1;
	}
	tr = cB(0);
	tr_error = tr.abs();
	if(!std::isfinite(tr_error)) BDMsP_TrcFail = 1;
	return tr;
}
Complex BDMode::travres(Complex g)
{
	beta = g.re;
	alpha = -g.im;
	return travres();
}

/* default profile crop radii: */
/* contributions for radii below 
   (radius of the innermost interface)-BDMp_R_CROP_IN*(vacuum wavelength) 
   will be neglected */
double BDMp_R_CROP_IN = 5.0;
/* contributions for radii above  
   (radius of the outermost interface)+BDMp_R_CROP_OUT*(vacuum wavelength)
   will be neglected */
double BDMp_R_CROP_OUT = 7.0;
/* relative field levels below BDMp_R_CROP_LVL will be neglected */
double BDMp_R_CROP_LVL = 1.0e-3;

/* mode profile evaluation: neglect fields at radii below 
   min(R, vac.wavelength)*BDMp_R_tiny */
double BDMp_R_TINY = 1.0e-2;

/* mode profile classification: relative node detection level */ 
double BDMp_NodeLVL = 1.0e-10;

/* radial profile, principal component */
Complex BDMode::profile(double r) const
{
	if(r < BDMp_R_TINY*R) return CC0;
	if(r < BDMp_R_TINY*wg.lambda) return CC0;
	int l = wg.layeridx(r-R);
	if(disc_prof == 0 || r <= d_ri || r >= d_re)
	{
		double rho = k0*wg.n(l)*r;
		Complex p;
			 p  = cA(l)*BfJ(c_m(), rho);
		if(l>=1) p += cB(l)*BfY(c_m(), rho);
		return p*nrm;
	}
	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 h = v_rad(l, khi)-v_rad(l, klo);
	double a = (v_rad(l, khi)-r)/h;
	double b = (r-v_rad(l, klo))/h;
	return v_prf(l, klo)*a + v_prf(l, khi)*b + ((a*a*a-a)*d_prf(l, klo)+(b*b*b-b)*d_prf(l, khi))*(h*h)/6.0;
}

/* ... and its derivative */
Complex BDMode::d_profile(double r) const
{
	if(r < BDMp_R_TINY*R) return CC0;
	if(r < BDMp_R_TINY*wg.lambda) return CC0;
	int l = wg.layeridx(r-R);
	if(disc_prof == 0 || r <= d_ri || r >= d_re)
	{
		double rho = k0*wg.n(l)*r;
		Complex p;
			 p  = cA(l)*dBfJ(c_m(), rho);
		if(l>=1) p += cB(l)*dBfY(c_m(), rho);
		return p*k0*wg.n(l)*nrm;
	}
	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 h = v_rad(l, khi)-v_rad(l, klo);
	double a = (v_rad(l, khi)-r)/h;
	double b = (r-v_rad(l, klo))/h;
	return v_dpr(l, klo)*a + v_dpr(l, khi)*b + ((a*a*a-a)*d_dpr(l, klo)+(b*b*b-b)*d_dpr(l, khi))*(h*h)/6.0;
}

/* initialize discrete profiles 
   on radial interval [ir < wg.hx(.) < re] with nr steps */
void BDMode::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)-BDMp_R_CROP_IN*wg.lambda; 
	if(ri < R*BDMp_R_TINY) ri = R*BDMp_R_TINY;
	if(re < R+wg.hx(wg.nx)) re = (wg.hx(wg.nx)+R)+BDMp_R_CROP_OUT*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_prf = Cmatrix(wg.nx+2, maxn);
	v_dpr = Cmatrix(wg.nx+2, maxn);
	d_prf = Cmatrix(wg.nx+2, maxn);
	d_dpr = Cmatrix(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)) benderror("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;
			v_prf(l, m) =   profile(v_rad(l, m)+tdr);
			v_dpr(l, m) = d_profile(v_rad(l, m)+tdr);
		}
	}
	for(int l=0; l<=wg.nx+1; ++l)
	{
		int n;
		Dvector x;
		Cvector y;
		Cvector y2;
		n = n_rad(l);
		Cvector u(n);
		y2 = Cvector(maxn);

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

		y = v_prf.row(l).subvector(n, 0);
		y2.init(CC0);
		u(0) = CC0;
		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) = CC0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_prf.row(y2, l);

		y = v_dpr.row(l).subvector(n, 0);
		y2.init(CC0);
		u(0) = CC0;
		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) = CC0;
		for(int k=n-2; k>=1; --k) y2(k) = y2(k)*y2(k+1)+u(k);
		d_dpr.row(y2, l);
	}
	disc_prof = 1;
	return;
}

/* BDM field */
Complex BDMode::field(Fcomp cp, double x, double z) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dz;
	double cosp, sinp, r, p;
	Complex prf, d_prf, efac;
	dx = x-x0;
	dz = z-z0;
	r = sqrt(dx*dx+dz*dz);
	if(r < BDMp_R_TINY*R) return CC0;
	if(r < BDMp_R_TINY*wg.lambda) return CC0;
	cosp = dx/r; 
	sinp = dz/r; 
	p = atan2(dz, dx); // if(p<0) p+=2.0*PI;
	prf   =   profile(r);
	d_prf = d_profile(r);
	efac = expmi(c_m()*p*s);
	double nq = wg.eps(r-R);

	Complex dxPy, dzPy;
	
	if(pol == TE) 
	{
		switch(cp)
		{
		case EX: case EZ: case HY: case ER: case EP:
			return CC0;
			break;
		case EY:
			return prf*efac;
			break;
		case HX:
			dzPy = (d_prf*sinp - CCI*c_m()*s/r*cosp*prf)*efac;
			return CCmI*dzPy*invommu0;
			break;
		case HZ:
			dxPy = (d_prf*cosp + CCI*c_m()*s/r*sinp*prf)*efac;
			return CCI*dxPy*invommu0;
			break;
		case HR: 
			return CCm1*invommu0*c_m()*s/r*prf*efac;
			break;
		case HP: 
			return CCI*invommu0*d_prf*efac;
			break;
		default:
			benderror("field: cp not implemented");
			break;
		}
		return CC0;
	}
	switch(cp)
	{
	case HX: case HZ: case EY: case HR: case HP:
		return CC0;
		break;
	case HY:
		return prf*efac;
		break;
	case EX:
		dzPy = (d_prf*sinp - CCI*c_m()*s/r*cosp*prf)*efac;
		return CCI*dzPy*invomep0/nq;
		break;
	case EZ:
		dxPy = (d_prf*cosp + CCI*c_m()*s/r*sinp*prf)*efac;
		return CCmI*dxPy*invomep0/nq;
		break;
	case ER: 
		return CC1*invomep0*c_m()*s/r*prf*efac/nq;
		break;
	case EP: 
		return CCmI*invomep0*d_prf*efac/nq;
		break;
	default:
		benderror("field: cp not implemented");
		break;
	}
	return CC0;
}
Complex BDMode::field(Fcomp cp, double r) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	if(r < BDMp_R_TINY*R) return CC0;
	if(r < BDMp_R_TINY*wg.lambda) return CC0;
	Complex prf, d_prf;
	prf   =   profile(r);
	d_prf = d_profile(r);
	double nq = wg.eps(r-R);

	Complex dxPy, dzPy;
	
	if(pol == TE) 
	{
		switch(cp)
		{
		case EX: case EZ: case HY: case ER: case EP:
			return CC0;
			break;
		case EY:
			return prf;
			break;
		case HX:
			dzPy = CCmI*c_m()*s/r*prf;
			return CCmI*dzPy*invommu0;
			break;
		case HZ:
			dxPy = d_prf;
			return CCI*dxPy*invommu0;
			break;
		case HR: 
			return CCm1*invommu0*c_m()*s/r*prf;
			break;
		case HP: 
			return CCI*invommu0*d_prf;
			break;
		default:
			benderror("field: cp not implemented");
			break;
		}
		return CC0;
	}
	switch(cp)
	{
	case HX: case HZ: case EY: case HR: case HP:
		return CC0;
		break;
	case HY:
		return prf;
		break;
	case EX:
		dzPy = CCmI*c_m()*s/r*prf;
		return CCI*dzPy*invomep0/nq;
		break;
	case EZ:
		dxPy = d_prf;
		return CCmI*dxPy*invomep0/nq;
		break;
	case ER: 
		return CC1*invomep0*c_m()*s/r*prf/nq;
		break;
	case EP: 
		return CCmI*invomep0*d_prf/nq;
		break;
	default:
		benderror("field: cp not implemented");
		break;
	}
	return CC0;
}
void BDMode::emfield(double x, double z, 
	             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, dz;
	double cosp, sinp, r, p;
	Complex prf, d_prf, efac;
	dx = x-x0;
	dz = z-z0;
	r = sqrt(dx*dx+dz*dz);
	if(r < BDMp_R_TINY*R || r < BDMp_R_TINY*wg.lambda)
	{
		ex = CC0; ey = CC0; ez = CC0; 
		hx = CC0; hy = CC0; hz = CC0;
		return;
	} 
	cosp = dx/r; 
	sinp = dz/r; 
	p = atan2(dz, dx); // if(p<0) p+=2.0*PI;
	prf   =   profile(r);
	d_prf = d_profile(r);
	efac = expmi(c_m()*s*p);
	double nq = wg.eps(r-R);

	Complex dxPy, dzPy;
	
	if(pol == TE) 
	{
		ex = CC0;
		ey = prf*efac; 
		ez = CC0;
		dzPy = (d_prf*sinp - CCI*c_m()*s/r*cosp*prf)*efac;
		hx = CCmI*dzPy*invommu0;
		dxPy = (d_prf*cosp + CCI*c_m()*s/r*sinp*prf)*efac;
		hz = CCI*dxPy*invommu0;
		return;
	}
	dzPy = (d_prf*sinp - CCI*c_m()*s/r*cosp*prf)*efac;
	ex = CCI*dzPy*invomep0/nq;
	ey = CC0;
	dxPy = (d_prf*cosp + CCI*c_m()*s/r*sinp*prf)*efac;
	ez = CCmI*dxPy*invomep0/nq;
	hx = CC0;
	hy = prf*efac;
	hz = CC0;
	return;
}

void BDMode::emfield_r(double x, double z, 
	               Complex &er, Complex &ey, Complex &ep, 
	               Complex &hr, Complex &hy, Complex &hp) const
{
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx = x-x0;
	double dz = z-z0;
	double r = sqrt(dx*dx+dz*dz);
	if(r < BDMp_R_TINY*R || r < BDMp_R_TINY*wg.lambda)
	{
		er = CC0; ey = CC0; ep = CC0; 
		hr = CC0; hy = CC0; hp = CC0;
		return;
	} 
	double p = atan2(dz, dx); // if(p<0) p+=2.0*PI;
	Complex prf   =   profile(r);
	Complex d_prf = d_profile(r);
	Complex efac = expmi(c_m()*p*s);
	double nq = wg.eps(r-R);
	
	if(pol == TE) 
	{
		er = CC0;
		ey = prf*efac;
		ep = CC0;
		hr = CCm1*invommu0*c_m()*s/r*prf*efac; 
		hy = CC0;
		hp = CCI*invommu0*d_prf*efac;
		return;
	}
	er = CC1*invomep0*c_m()*s/r*prf/nq*efac;
	ey = CC0;
	ep = CCmI*invomep0*d_prf/nq*efac;
	hr = CC0;
	hy = prf*efac;
	hp = CC0;
	return;
}

/* ... values on a rectangular mesh */
Cmatrix BDMode::field(Fcomp cp, Interval dwx, int npx, 
                                Interval dwz, int npz) const
{
	if(npx <= 1 || npz <= 1) benderror("discretized field: npx, npz"); 
	Cmatrix fld(npx, npz);
	double dx = dwx.len()/((double)(npx-1));
	double dz = dwz.len()/((double)(npz-1));
	double x, z;

	z = dwz.x0;
	for(int j=0; j<=npz-1; ++j)
	{
		x = dwx.x0;
		for(int i=0; i<=npx-1; ++i) 
		{
			fld(i, j) = field(cp, x, z);
			x += dx;
		}
		z += dz;
	}
	return fld;
}

/* determine radial mode order, set the mode id string */
#define BDMp_NumRDStepspWl 50.0;
void BDMode::classify()
{
	disc_prof = 0;
	v_rad.strip();
	v_prf.strip();
	v_dpr.strip();
	d_prf.strip();
	d_dpr.strip();

	double rm = R+wg.hx(wg.nx);
	double re = (wg.hx(wg.nx)+R)+BDMp_R_CROP_OUT*(wg.lambda);
	double dr = fabs(wg.lambda/wg.n.max())/BDMp_NumRDStepspWl;
	double pvm = profile(R).abs();
	double rb = R*BDMp_R_TINY; 
	while(rb < R+wg.hx(0)-10.0*dr && profile(rb).abs() < pvm*BDMp_R_CROP_LVL) rb += 5.0*dr;
	if(rb > R+wg.hx(0)) rb = (wg.hx(0)+R)-BDMp_R_CROP_IN*wg.lambda; 
	if(rb < R*BDMp_R_TINY) rb = R*BDMp_R_TINY;

	int npi = ((int)ceil((re-rb)/dr)); 
	Dvector pfl(npi); 
	Dvector ffl(npi);
	Dvector rps(npi); 
	double ph = profile(R).arg();	
	Complex phcf = expmi(ph);
	cRi = rb;
	cRe = re;
	double r0 = rb;
	double ipfmax = 0.0;
	double fa;
	for(int i=0; i<=npi-1; ++i)
	{
		Complex f = profile(r0)*phcf; 
		rps(i) = r0;
		ffl(i) = f.re; 
		fa = f.abs();
		pfl(i) = fa;
		if(r0 < rm && ipfmax < fa) ipfmax = fa;	
		r0 += dr;
	}

	ord = 0;
//	int ili = 0;
	double zl = ipfmax*BDMp_NodeLVL;
	for(int i=0; i<=npi-2; ++i)
	{
		if(rps(i+1)<=rm+HDIST && ffl(i)*ffl(i+1)<-zl) ++ord; 
//		if(rps(i)<=rm+dr+HDIST) ili=i;
	}
//	if(ili==0 || pfl(ili+1) > pfl(ili)+zl) ord = -1;

	int i=0;
        ids[i++] = 'T';
        if(pol == TE) ids[i++] = 'E';
        else          ids[i++] = 'M';
	if(ord >= 0)
	{
		if(ord >=100) ids[i++] = dig100(ord);
		if(ord >=10)  ids[i++] = dig10(ord);
		ids[i++] = dig1(ord);
	}
	else                  ids[i++] = 'x';
	if(dir == BACK) ids[i++] = '-';
        while(i<=15) ids[i++] = 0;
	
	double cl = pfl.max()*BDMp_R_CROP_LVL;
	i = 0;
	while(pfl(i) < cl) ++i;
	cRi = rb+((double) i)*dr;
	i = pfl.nel-1;
	while(pfl(i) < cl) --i;
	cRe = rb+((double) i)*dr;
        return;
}

/* power transported in angular direction 
   on infinite radial interval at angle theta */
double BDMode::power(double theta) const
{
	double s = 1.0; if(dir == BACK) s = -1.0;
	double p;
	if(fabs(c_neff().im) > 1.0e-8)
	{
		Complex f = profile(R);
		double rho = k0*wg.n(wg.nx+1)*R;
		f = f/(BfJ(c_m(), rho)-CCI*BfY(c_m(), rho));
		p = f.sqabs()/2.0/alpha/R/PI;
		if(pol == TE) p *= val_invommu0(wg.lambda);
		else          p *= val_invomep0(wg.lambda)/wg.eps(wg.nx+1);	
		p *= exp(alpha*R*(PI-2.0*s*theta));
		return s*p;
	}
	
        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};

	p =0.0;
	double ri = cRi;
	if(ri > R+wg.hx(0)) ri = (wg.hx(0)+R)-BDMp_R_CROP_IN*wg.lambda; 
	if(ri < R*BDMp_R_TINY) ri = R*BDMp_R_TINY;
	double re = cRe;
	if(re < R+wg.hx(wg.nx)) re = (wg.hx(wg.nx)+R)+BDMp_R_CROP_OUT*wg.lambda; 
	double ds0 = (re-ri)/70.0;

	for(int l=0; l<=wg.nx+1; ++l)
	{
		double epsl = wg.eps(l);
		double ra, rb;
		if(l == 0) ra = ri; 
		else       ra = wg.hx(l-1)+R;
		if(l == wg.nx+1) rb = re; 
		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 = profile(rx).sqabs()/rx;
				if(pol == TM) fp = fp/epsl;
				rx = xm-dx;
				double fm = profile(rx).sqabs()/rx;
				if(pol == TM) fm = fm/epsl;
				s += gipw[j]*(fp+fm);
			}
			p += s*xr;
			a = b;
		}
		
	}		
	p *= beta*R/2.0;
	if(pol == TE) p *= val_invommu0(wg.lambda);
	else          p *= val_invomep0(wg.lambda);	
	p *= exp(-2.0*alpha*R*theta*s);
	return s*p;
}

/* normalize mode to unit angular power at angular origin */ 
void BDMode::normalize()
{
	nrm = CC1;
	double p = power(0.0);
	nrm = Complex(1.0/sqrt(fabs(p)), 0.0);
	double ph = profile(R).arg();	
	nrm = nrm*expmi(ph);
	return;
}

/* translate the BDM: shift the origin by dx, dz */
void BDMode::translate(double dx, double dz)
{
	x0 += dx; 
	z0 += dz; 
	for(int o=0; o<=bend.no-1; ++o)	
	{
		switch(bend(o).type) 
		{
		case RING: case DISK:
			bend(o).xo += dx;
			bend(o).zo += dz;
			break;
		default:
			benderror("translate: bend.type");
			break;
		}
	}
	return;
}

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

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

/* BDM field profile plots, image of component cp 
   cp: one of EX, EY, EZ, HX, HY, HZ, ER, EP, HR, HP
   foa: MOD, SQR, REP, IMP, ORG = REP
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters for the m.file */
void BDMode::plot(Fcomp cp, Afo foa,
                  double xbeg, double xend, double zbeg, double zend, 
                  int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "pl_____I.m";
	name[2] = afochr(foa);
        name[3] = fldchr(cp);
        name[4] = cpchr(cp);      
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "plot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "BDM field");
	mlout_gengeoxz(dat, bend, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix cfld = field(cp, Interval(xbeg, xend), npx, 
	                         Interval(zbeg, zend), npz);
	Dmatrix fld = realize(cfld, foa); 
	mlout_fld(dat, npx, npz, cp, fld);
	name[8] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* export full field data into a viewer m-file
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters */
void BDMode::viewer(double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "fld____A.m";
	name[3] = 't';
	name[4] = polchr(pol);
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   xbeg, npx, xend, zbeg, npz, zend, name);
	dat = fopen(name, "w+");
	name[8] = 0; 
	mlout_viewertop(dat, name, pol, wg.lambda);
	mlout_gengeoxz(dat, bend, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix f;
	Fcomp cp;
	if(pol == TE)
	{
		cp = EX;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EZ;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HX;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HY;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HZ;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);

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

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



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

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

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

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

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

}

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

}

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

/* add a mode */
void BDModeArray::add(BDMode m)
{
	BDMode *avec;
	avec = new BDMode[num+1];
	BDMode* ap = avec;
	BDMode* 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 BDModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		benderror("BDModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	BDMode *avec;
	avec = new BDMode[num-1];
	BDMode* ap = avec;
	BDMode* 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 BDModeArray nma */ 
void BDModeArray::merge(BDModeArray& ma)
{
	BDMode *avec;
	avec = new BDMode[num+ma.num];
	BDMode* ap = avec;
	BDMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	BDMode* 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;
}

/* sort the array by phase propagation constants, highest first */
void BDModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	BDMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		maxi = j; 
		maxb = (mvec[j]).beta;
		for(k=j+1; k<=num-1; ++k)
		{
			if(maxb<(mvec[k]).beta)
			{
				maxb = (mvec[k]).beta;
				maxi = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t;
	}
	return;
}

/* sort the array by mode order, lowest first, invalid last */
void BDModeArray::order()
{
	int j, k, mini;
	int mino;
	BDMode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		mini = j; 
		mino = (mvec[j]).ord;
		for(k=j+1; k<=num-1; ++k)
		{
			if((mvec[k]).ord >= 0 && ((mvec[k]).ord < mino || mino <0))
			{
				mini = k; 
				mino = (mvec[k]).ord;
			}
		}
		t = mvec[j]; mvec[j] = mvec[mini]; mvec[mini] = t;
	}
	for(j=0; j<=num-2; ++j)
	{
		mini = j; 
		mino = (mvec[j]).ord;
		if(mino < 0)	
		{
			for(k=j+1; k<=num-1; ++k)
			{
//				if(mvec[k].ord < 0 && mvec[k].tr_error < mvec[j].tr_error)
				if(mvec[k].ord < 0 && mvec[k].beta > mvec[j].beta)
				{
					mini = k; 
					mino = (mvec[k]).ord;
				}
			}
			t = mvec[j]; mvec[j] = mvec[mini]; mvec[mini] = t;
		}
	}
	return;
}

/* ------------------------------------------------------------------------
   BDM mode analysis procedures
   ------------------------------------------------------------------------ */

/* Heuristics for mode finding: 
   mode order identification by solving an equivalent straight problem 
   with Dirichlet boundary conditions */

// staircase approximation of the equiv. st. problem, number of slices
int    BDMsP_NumSlices = 100;
// radial distance of interior boundary to innermost radial interface, in vacuum wavelengths
double BDMsP_IntBd = 2.5;
// radial distance of exterior boundary to outermost radial interface, in vacuum wavelengths
double BDMsP_ExtBd = 2.5;
// number of spectral terms examined
int BDMsP_NSpT = 99;
// minimum effective index (real part) considered
double BDMsP_MinNeff = 0.1;

/* complex root finding */

// initial attenuation, -log10(upper level) */
double BDMsP_IniMaxAtt = 2;
// initial attenuation, -log10(lower level) */
double BDMsP_IniMinAtt = 14;
// secant method, terminal distance of root prototypes
double BDMsP_RootAcc = 1.0e-13;
// number of trials for convergence with modified initial conditions
int BDMsP_LimNumIts = 30;
// maximum growth of trial value distance permitted per iteration step
double BDMsP_ConvGrowth = 2.0;
// minimum difference in absolute effective mode index for two modes
// to be considered different
double BDMsP_MinNdiff = 1.0e-4;
// maximum acceptable level of absolute TRC violation 
double BDMsP_MaxTRC = 1.0e-3;


/* --- solver procedure ---------------------------------------------------*/ 

/* for the bend with radial layering w, radius r, centered at (xc, zc),  
   try to find pol-polarized BDMs of lowest order, 
   consider ntre * ntim trial values for complex effective mode indices,
   lorm=1: consider modes of lowest radial order only,  
   quiet = 0, 1, 2: all, less, no  log output  */
int bdmsolve(Waveguide w, double r, double xc, double zc, Polarization pol, 
             int lorm, int ntre, int ntim, BDModeArray& bdma, int quiet)
{
// quiet = 0; 
	if(quiet < 2)
	{
fprintf(stderr, "\n------------- BDMsolve --------------------------------------------- '18 ---\n");
fprintf(stderr, "T%c  lambda: %g  ", polCHR(pol), w.lambda);
fprintf(stderr, "R: %g   (x0, z0): (%g, %g)\n", r, xc, zc);
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");
	}

	Dvector tbv;
	double db;

	if(ntre < 0)
	{

		ModeArray ma, mas;
		mas.clear();

		double ib, eb, ai, ae;
		ib = r+w.hx(0)-BDMsP_IntBd*w.lambda; if(ib < r/10.0) ib = r/10.0; 
		eb = r+w.hx(w.nx)+BDMsP_ExtBd*w.lambda; 
		ai = r*log(ib/r);
		ae = r*log(eb/r);
		int nx = BDMsP_NumSlices;	

		Waveguide eqwg(nx);
		eqwg.special = 1;
		double da = (ae-ai)/((double) nx+2);
		for(int i=0; i<=nx; ++i) eqwg.hx(i) = ai+da + ((double) i)*da;
		double maxe = 0.0;
		for(int l=0; l<=w.nx; ++l)
		{
		
			double ta = r*log((r+w.hx(l))/r);
			double me = w.eps(l)*exp(2.0*ta/r);
			if(me > maxe) maxe = me;
		}
		for(int i=0; i<=nx+1; ++i) 
		{
			double a;
			if(i<=0) a = eqwg.hx(0)-da/2.0;
			else
			{
				if(i>=nx) a = eqwg.hx(nx)+da/2.0;
				else a = 0.5*(eqwg.hx(i)+eqwg.hx(i+1));
			}
			double p;
			p = r*exp(a/r);
			eqwg.n(i) = w.eps(p-r)*exp(2.0*a/r);
			if(eqwg.n(i) > maxe) eqwg.n(i) = maxe; 
		}
		eqwg.lambda = w.lambda;
		
		modespectrum(eqwg, pol, Interval(ai, ae), LIMD, BDMsP_NSpT, ma, 1);
		double wnm;
		if(lorm == 1) wnm = w.n.min()*val_k0(w.lambda);
		else          wnm = BDMsP_MinNeff*val_k0(w.lambda);
		wnm = wnm*wnm;
		for(int j=0; j<=ma.num-1; ++j) 
		{
                        if(ma(j).betaq >= wnm) mas.add(ma(j));
		}

		ma = mas; 
		mas.clear();

		ma.sort();
		db = 0.5;
		for(int i=0; i<=ma.num-2; ++i)
		{
			double d = fabs(ma(i).beta-ma(i+1).beta);
			if(d < db) db = d;
		}
		db = db/3.0;

if(quiet < 1)
{
		fprintf(stderr, "Eq.st.wg, LimDbc @ %.4g, %.4g -> %d candidates.\n", ib, eb, ma.num);
		fprintf(stderr, "----------------------------------------------------------------------------\n");
}

		tbv = Dvector(ma.num);
		for(int i=0; i<=ma.num-1; ++i) tbv(i) = ma(i).beta;
	}
	else
	{
		double mib;
                if(lorm == 1)
                {
                        mib = w.n.min()*val_k0(w.lambda);
                }
                else
                {
                        mib = BDMsP_MinNeff*val_k0(w.lambda);
                }
		double mab = w.n.max()*val_k0(w.lambda); 
		tbv = Dvector(ntre);
		for(int i=0; i<=ntre-1; ++i) 
			tbv(i) = mib + (mab-mib)/((double) (ntre-1))*((double) i);
		db = (mab-mib)/((double) (ntre-1))/3;
	}

	ModeArray mg;
	if(w.nx >= 1)
	{
		modeanalysis(w, pol, mg, 1);
		for(int j=0; j<=mg.num-1; ++j) tbv.append(mg(j).beta);
	}

	if(ntim <= 0) ntim = 10;
	double astep = (BDMsP_IniMinAtt-BDMsP_IniMaxAtt)/((double) ntim);

	if(quiet < 2)
	{	
	fprintf(stderr, "neff-lookup: %d x %d trials\n", tbv.nel, ntim);
	fprintf(stderr, "----------------------------------------------------------------------------\n");
	}


	int trcE;
	bdma.clear();

	BDMode bdm;
	bdm.pol = pol;
	bdm.R = r;
	bdm.x0 = xc;
	bdm.z0 = zc;
	bdm.wg = w;
	bdm.bend = OvlStruct(w.n(w.nx+1), w.lambda); 
	bdm.bend.clear();
	for(int l=w.nx; l>=0; --l)
	{
		if(fabs(w.n(l+1)-w.n(l)) > COMPTOL_N)
		{
			Overlay o(DISK, w.n(l), r+w.hx(l), xc, zc);
			bdm.bend.place(o);
		}
	}	
	bdm.k0 = val_k0(w.lambda);
	bdm.invommu0 = val_invommu0(w.lambda);
	bdm.invomep0 = val_invomep0(w.lambda);
	bdm.dir = FORW;
	bdm.cRi = (bdm.wg.hx(0)+bdm.R)-BDMp_R_CROP_IN*bdm.wg.lambda; 
	if(bdm.cRi < bdm.R*BDMp_R_TINY) bdm.cRi = bdm.R*BDMp_R_TINY;
	bdm.cRe = (bdm.wg.hx(bdm.wg.nx)+bdm.R)+BDMp_R_CROP_OUT*bdm.wg.lambda; 
	bdm.nrm = CC1;
	bdm.cA = Cvector(w.nx+2);
	bdm.cB = Cvector(w.nx+2);

	for(int oi=0; oi<=1; ++oi)
	{
	for(double ai=BDMsP_IniMaxAtt; ai<=BDMsP_IniMinAtt; ai += astep)
	{
		double maxatt = pow(10.0, -ai);
		double minatt = pow(10.0, -(ai+astep));

	for(int mi=0; mi<=tbv.nel-1; ++mi)
	{
		trcE = 0;

		double tb = tbv(mi);
if(quiet < 1)
{
		fprintf(stderr, "[neff = %.4g-i%.4g]:\n", tb/val_k0(w.lambda), 0.5*(maxatt+minatt)/val_k0(w.lambda));
}
		bdm.tr_init(Complex(tb, -0.5*(maxatt+minatt)));
		
		Complex z, z1, z2;
		Complex f1, f2, f;
		if(oi==0)
		{
			z1 = Complex(tb-db, -maxatt); 
			z2 = Complex(tb+db, -minatt); 
		}
		else 
		{
			z1 = Complex(tb+db, -maxatt); 
			z2 = Complex(tb-db, -minatt); 
		}

		f1 = bdm.travres(z1); if(BDMsP_TrcFail != 0) ++trcE;
		f2 = bdm.travres(z2); if(BDMsP_TrcFail != 0) ++trcE;
		if(f1.abs() < f2.abs())
		{
			f  = f2; z  = z2; 
			f2 = f1; z2 = z1; 
			f1 = f;  z1 = z; 
		}
if(quiet < 1) fprintf(stderr, "   i: "); 
		int nit = 0;
		while(trcE == 0 && ((z1-z2).abs() > BDMsP_RootAcc) && nit < BDMsP_LimNumIts)
		{
			z = (z2*f1-z1*f2)/(f1-f2);
if(quiet < 1) fprintf(stderr, "c"); 
			if(z.im >= 0.0) 
			{
				z.im = 0.5*z2.im;
if(quiet < 1) fprintf(stderr, "s"); 
			}
			if((z-z2).abs() > BDMsP_ConvGrowth*(z2-z1).abs())
			{
				z = z2 + BDMsP_ConvGrowth*((z2-z1).abs()/(z-z2).abs())*(z-z2);
if(quiet < 1) fprintf(stderr, "g"); 
			}
			f = bdm.travres(z); if(BDMsP_TrcFail != 0) ++trcE;
			f1 = f2; z1 = z2; 
			f2 = f;  z2 = z; 
			if(f1.abs() < f2.abs())
			{
if(quiet < 1) fprintf(stderr, "r"); 
				f  = f2; z  = z2; 
				f2 = f1; z2 = z1; 
				f1 = f;  z1 = z; 
			}
if(quiet < 1) fprintf(stderr, " "); 
			++nit;
		}
if(quiet < 1) fprintf(stderr, "\n   "); 
		if(trcE > 0)
		{
if(quiet < 1) fprintf(stderr, "F.\n"); 
		}
// if(quiet < 1) fprintf(stderr, "trcE: %d, |dz| = %g, Im(z) = %g, |TRC| = %g\n", trcE, (z1-z2).abs(), z2.im, f2.sqabs()); 
		if(trcE == 0 && ((z1-z2).abs() <= BDMsP_RootAcc) && (z2.im < 0.0) && f2.sqabs() < BDMsP_MaxTRC)
		{
			bdm.nrm = CC1;
			bdm.beta = z2.re; 
			bdm.alpha = -z2.im; 
			Complex q;
			q = bdm.travres(); if(BDMsP_TrcFail != 0) ++trcE;
			bdm.classify();
			bdm.normalize();
if(quiet < 1)
{
			Complex neff = bdm.c_neff();
			fprintf(stderr, "-> %s, neff=%.4g+i%.4g, errTRC: %.4g", bdm.ids, neff.re, neff.im, bdm.tr_error);
}
			if(trcE == 0)
			{
				bdma.add(bdm);
if(quiet < 1) fprintf(stderr, ".\n");
			}
			else
			{
if(quiet < 1) fprintf(stderr, " (-).\n");
			}
		}
		else
		{
if(quiet < 1) fprintf(stderr, "-> -.\n"); 
		}
	}
	}
	}
if(quiet < 1)
{
	fprintf(stderr, "----------------------------------------------------------------------------\n");
}
	bdma.sort();
        int done = 0;
	while(done == 0 && bdma.num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=bdma.num-1; ++j)
		{
			Complex tn = bdma(j).c_neff();
			for(int k=0; k<=bdma.num-1; ++k)
			{
				Complex nn = bdma(k).c_neff();
				if(j!=k && done == 1 && (tn-nn).abs() < BDMsP_MinNdiff)
				{
					if(bdma(j).tr_error <= bdma(k).tr_error) rmm = k;
					else                                     rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) bdma.remove(rmm);
	}
	bdma.sort();
        done = 0;
	while(done == 0 && bdma.num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=bdma.num-1; ++j)
		{
			int to = bdma(j).ord;
			for(int k=0; k<=bdma.num-1; ++k)
			{
				int no = bdma(k).ord;
				if(j!=k && done == 1 && to >= 0 && to == no)
				{
					if(bdma(j).tr_error <= bdma(k).tr_error) rmm = k;
					else                                     rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) bdma.remove(rmm);
	}
	bdma.order();
if(quiet < 2)
{
	for(int i=0; i<=bdma.num-1; ++i)
	{
		fprintf(stderr, "[%d] %s, neff=%g+i%g, beta=%g, alpha=%g\n", 
		i, bdma(i).ids, bdma(i).c_neff().re, bdma(i).c_neff().im, bdma(i).beta, bdma(i).alpha);
	}
	fprintf(stderr, "Ok.\n\n");
}
	return bdma.num;
}