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

/*
 * wgm.cpp
 * whispering-gallery modes of circular cavities in 2-D 
 */

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

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

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

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


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

/* determine suitable initial amplitude */
void WGMode::tr_init(Complex o)
{
	omega = o.re;
	alpha = o.im;
	double r0 = R+wg.hx(wg.nx);
	Complex rho0 = c_k()*wg.n(wg.nx+1)*r0;
	Complex J0 = BfJ(m, rho0); if(WGMCylEvalError) WGMsP_TrcFail = 1;
	Complex Y0 = BfY(m, rho0); if(WGMCylEvalError) WGMsP_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 WGMode::travres()
{
	WGMsP_TrcFail = 0;
	cA(wg.nx+1) = tr_Ainit; 
	cB(wg.nx+1) = CCmI*cA(wg.nx+1);
	double r;
	Complex rho;
	Complex tr;
	if(pol == TE) 
	{
		for(int l=wg.nx; l>=0; --l)
		{	
			r = R+wg.hx(l); 
			rho = c_k()*wg.n(l)*r;
			Complex m_a =  BfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_b =  BfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_c = dBfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_d = dBfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			m_c = m_c*wg.n(l);
			m_d = m_d*wg.n(l);
			rho = c_k()*wg.n(l+1)*r;
			Complex m_e =  BfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_f =  BfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_g = dBfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
			Complex m_h = dBfY(m, rho); if(WGMCylEvalError) WGMsP_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()) WGMsP_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()) WGMsP_TrcFail = 1;
			if(cB(l).isnotOK()) WGMsP_TrcFail = 1;
		}
		tr = cB(0);
		tr_error = tr.abs();
		if(!std::isfinite(tr_error)) WGMsP_TrcFail = 1;
		return tr;
	}
	for(int l=wg.nx; l>=0; --l)
	{	
		r = R+wg.hx(l); 
		rho = c_k()*wg.n(l)*r;
		Complex m_a =  BfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_b =  BfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_c = dBfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_d = dBfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		m_c = m_c/wg.n(l);
		m_d = m_d/wg.n(l);
		rho = c_k()*wg.n(l+1)*r;
		Complex m_e =  BfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_f =  BfY(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_g = dBfJ(m, rho); if(WGMCylEvalError) WGMsP_TrcFail = 1;
		Complex m_h = dBfY(m, rho); if(WGMCylEvalError) WGMsP_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()) WGMsP_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()) WGMsP_TrcFail = 1;
		if(cB(l).isnotOK()) WGMsP_TrcFail = 1;
	}
	tr = cB(0);
	tr_error = tr.abs();
	if(!std::isfinite(tr_error)) WGMsP_TrcFail = 1;
	return tr;
}
Complex WGMode::travres(Complex o)
{
	omega = o.re;
	alpha = o.im;
	return travres();
}

/* converge transverse resonance condition towards root;
   return 1: success; 0: failure */
int WGMode::converge(Complex z1, Complex z2, int quiet)
{
	Complex f1, f2, f;
	Complex z = (z1+z2)*0.5;
if(quiet < 1)
{
	fprintf(stderr, "[m=%d, om=%.4g+i%.4g]:\n", m, z.re, z.im);
}
	int trcE = 0;
	tr_init(z);

	f1 = travres(z1); if(WGMsP_TrcFail != 0) ++trcE;
	f2 = travres(z2); if(WGMsP_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() > WGMsP_RootAcc) && nit < WGMsP_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() > WGMsP_ConvGrowth*(z2-z1).abs())
		{
			z = z2 + WGMsP_ConvGrowth*((z2-z1).abs()/(z-z2).abs())*(z-z2);
if(quiet < 1) fprintf(stderr, "g"); 
		}
		f = travres(z); if(WGMsP_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() <= WGMsP_RootAcc) && (z2.im > 0.0) && f2.sqabs() < WGMsP_MaxTRC)
	{
		omega = z2.re; 
		alpha = z2.im;
		Complex q;
		tr_init(Complex(omega, alpha));
		q = travres(); if(WGMsP_TrcFail != 0) ++trcE;
		cav.lambda = val_lambda(omega);
		wg.lambda = val_lambda(omega);
		classify();
if(quiet < 1)
{
		fprintf(stderr, "-> %s, om=%.4g+i%.4g, errTRC: %.4g", ids, omega, alpha, tr_error);
}
		if(trcE == 0)
		{
if(quiet < 1) fprintf(stderr, ".\n");
			return 1;
		}
		else
		{
if(quiet < 1) fprintf(stderr, " (-).\n");
			return 0;
		}
	}
	else
	{
if(quiet < 1) fprintf(stderr, "-> -.\n"); 
		return 0;
	}
	return 0;
}

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

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

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

/* radial profile, principal component */
Complex WGMode::profile(double r) const
{
	if(r < WGMp_R_TINY*R) return CC0;
	if(r < WGMp_R_TINY*wg.lambda) return CC0;
	int l = wg.layeridx(r-R);
	if(disc_prof == 0 || r <= d_ri || r >= d_re)
	{
		Complex rho = c_k()*wg.n(l)*r;
		Complex p;
			 p  = cA(l)*BfJ(m, rho);
		if(l>=1) p += cB(l)*BfY(m, rho);
		return p;
	}
	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 WGMode::d_profile(double r) const
{
	if(r < WGMp_R_TINY*R) return CC0;
	if(r < WGMp_R_TINY*wg.lambda) return CC0;
	int l = wg.layeridx(r-R);
	if(disc_prof == 0 || r <= d_ri || r >= d_re)
	{
		Complex rho = c_k()*wg.n(l)*r;
		Complex p;
			 p  = cA(l)*dBfJ(m, rho);
		if(l>=1) p += cB(l)*dBfY(m, rho);
		return p*c_k()*wg.n(l);
	}
	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 WGMode::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)-WGMp_R_CROP_IN*wg.lambda; 
	if(ri < R*WGMp_R_TINY) ri = R*WGMp_R_TINY;
	if(re < R+wg.hx(wg.nx)) re = (wg.hx(wg.nx)+R)+WGMp_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)) wgmerror("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;
}

/* WGM field, time t=0 */
Complex WGMode::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 < WGMp_R_TINY*R || r < WGMp_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 = expi(-((double)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*((double)m)*s/r*cosp*prf)*efac;
			return CCmI*dzPy*invommu0();
			break;
		case HZ:
			dxPy = (d_prf*cosp + CCI*((double)m)*s/r*sinp*prf)*efac;
			return CCI*dxPy*invommu0();
			break;
		case HR: 
			return CCm1*invommu0()*((double)m)*s/r*prf*efac;
			break;
		case HP: 
			return CCI*invommu0()*d_prf*efac;
			break;
		default:
			wgmerror("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*((double)m)*s/r*cosp*prf)*efac;
		return CCI*dzPy*invomep0()/nq;
		break;
	case EZ:
		dxPy = (d_prf*cosp + CCI*((double)m)*s/r*sinp*prf)*efac;
		return CCmI*dxPy*invomep0()/nq;
		break;
	case ER: 
		return CC1*invomep0()*((double)m)*s/r*prf*efac/nq;
		break;
	case EP: 
		return CCmI*invomep0()*d_prf*efac/nq;
		break;
	default:
		wgmerror("field: cp not implemented");
		break;
	}
	return CC0;
}
Complex WGMode::field(Fcomp cp, double r) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	if(r < WGMp_R_TINY*R) return CC0;
	if(r < WGMp_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*((double)m)*s/r*prf;
			return CCmI*dzPy*invommu0();
			break;
		case HZ:
			dxPy = d_prf;
			return CCI*dxPy*invommu0();
			break;
		case HR: 
			return CCm1*invommu0()*((double)m)*s/r*prf;
			break;
		case HP: 
			return CCI*invommu0()*d_prf;
			break;
		default:
			wgmerror("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*((double)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()*((double)m)*s/r*prf/nq;
		break;
	case EP: 
		return CCmI*invomep0()*d_prf/nq;
		break;
	default:
		wgmerror("field: cp not implemented");
		break;
	}
	return CC0;
}
void WGMode::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 < WGMp_R_TINY*R || r < WGMp_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 = expi(-((double)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*((double)m)*s/r*cosp*prf)*efac;
		hx = CCmI*dzPy*invommu0();
		dxPy = (d_prf*cosp + CCI*((double)m)*s/r*sinp*prf)*efac;
		hz = CCI*dxPy*invommu0();
		return;
	}
	dzPy = (d_prf*sinp - CCI*((double)m)*s/r*cosp*prf)*efac;
	ex = CCI*dzPy*invomep0()/nq;
	ey = CC0;
	dxPy = (d_prf*cosp + CCI*((double)m)*s/r*sinp*prf)*efac;
	ez = CCmI*dxPy*invomep0()/nq;
	hx = CC0;
	hy = prf*efac;
	hz = CC0;
	return;
}
void WGMode::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 < WGMp_R_TINY*R || r < WGMp_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 = expi(-((double)m)*s*p);
	double nq = wg.eps(r-R);
	
	if(pol == TE) 
	{
		er = CC0;
		ey = prf*efac;
		ep = CC0;
		hr = CCm1*invommu0()*((double)m)*s/r*prf*efac; 
		hy = CC0;
		hp = CCI*invommu0()*d_prf*efac;
		return;
	}
	er = CC1*invomep0()*((double)m)*s/r*prf/nq*efac;
	ey = CC0;
	ep = CCmI*invomep0()*d_prf/nq*efac;
	hr = CC0;
	hy = prf*efac;
	hp = CC0;
	return;
}

/* directional interference: mix with a*(reversed mode) */
Complex WGMode::field_ri(Fcomp cp, double x, double z, Complex a) const
{ 
	double s = 1.0; if(dir == BACK) s = -1.0;
	double dx, dz;
	double r, p;
	Complex prf, d_prf;
	dx = x-x0;
	dz = z-z0;
	r = sqrt(dx*dx+dz*dz);
	if(r < WGMp_R_TINY*R || r < WGMp_R_TINY*wg.lambda) return CC0;
	p = atan2(dz, dx); if(p<0) p+=2.0*PI;
	prf   =   profile(r);
	d_prf = d_profile(r);
	double sm = -s;
	Complex efac  = expi(-((double)m)*s*p);
	Complex efacm = expi(-((double)m)*sm*p);
	double nq = wg.eps(r-R);

	if(pol == TE) 
	{
		switch(cp)
		{
		case HY: case ER: case EP:
			return CC0;
			break;
		case EY:
			return prf*(efac+a*efacm);
			break;
		case HR: 
			return CCm1*invommu0()
			       *prf*((double)m)/r*(s*efac+a*sm*efacm); 
			break;
		case HP: 
			return CCI*invommu0()*d_prf*(efac+a*efacm);
			break;
		default:
			wgmerror("field: cp not implemented");
			break;
		}
		return CC0;
	}
	switch(cp)
	{
	case EY: case HR: case HP:
		return CC0;
		break;
	case HY:
		return prf*(efac+a*efacm);
		break;
	case ER: 
		return CC1*invomep0()*((double)m)/r*prf/nq*(s*efac+a*sm*efacm);
		break;
	case EP: 
		return CCmI*invomep0()*d_prf/nq*(efac+a*efacm);
		break;
	default:
		wgmerror("field: cp not implemented");
		break;
	}
	return CC0;
}
void WGMode::emfield_ri(double x, double z, Complex a, 
	               Complex &er, Complex &ey, Complex &ep, 
	               Complex &hr, Complex &hy, Complex &hp) const
{
	if(a.abs() < 1.0e-10)
	{
		emfield_r(x, z, er, ey, ep, hr, hy, hp);
		return;
	}
	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 < WGMp_R_TINY*R || r < WGMp_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);
	double sm = -s;
	Complex efac  = expi(-((double)m)*s*p);
	Complex efacm = expi(-((double)m)*sm*p);
	double nq = wg.eps(r-R);
	
	if(pol == TE) 
	{
		er = CC0;
		ey = prf*(efac+a*efacm);
		ep = CC0;
		hr = CCm1*invommu0()*prf*((double)m)/r*(s*efac+a*sm*efacm); 
		hy = CC0;
		hp = CCI*invommu0()*d_prf*(efac+a*efacm);
		return;
	}
	er = CC1*invomep0()*((double)m)/r*prf/nq*(s*efac+a*sm*efacm);
	ey = CC0;
	ep = CCmI*invomep0()*d_prf/nq*(efac+a*efacm);
	hr = CC0;
	hy = prf*(efac+a*efacm);
	hp = CC0;
	return;
}

/* ... values on a rectangular mesh */
Cmatrix WGMode::field(Fcomp cp, Interval dwx, int npx, 
                                Interval dwz, int npz) const
{
	if(npx <= 1 || npz <= 1) wgmerror("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 WGMp_NumRDStepspWl 50.0;
void WGMode::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)+WGMp_R_CROP_OUT*(wg.lambda);
        double dr = fabs(wg.lambda/wg.n.max())/WGMp_NumRDStepspWl;
        double pvm = profile(R).abs();
        double rb = R*WGMp_R_TINY;
        while(rb < R+wg.hx(0)-10.0*dr && profile(rb).abs() < pvm*WGMp_R_CROP_LVL) rb += 5.0*dr;
        if(rb > R+wg.hx(0)) rb = (wg.hx(0)+R)-WGMp_R_CROP_IN*wg.lambda;
        if(rb < R*WGMp_R_TINY) rb = R*WGMp_R_TINY;

/*
	double rb = (wg.hx(0)+R)-WGMp_R_CROP_IN*(wg.lambda); 
	if(rb < R/10.0) rb = R/10.0;
	double re = (wg.hx(wg.nx)+R)+WGMp_R_CROP_OUT*(wg.lambda);
	double rm = R+wg.hx(wg.nx);
	double dr = fabs(wg.lambda/wg.n.max())/WGMp_NumRDStepspWl;
*/
	
	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*WGMp_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';
	ids[i++] = '(';
        if(ord >=100)      ids[i++] = dig100(ord);
        if(ord >=10)       ids[i++] = dig10(ord);
        ids[i++] = dig1(ord);
	ids[i++] = ',';
        if(dir == BACK)    ids[i++] = '-';
        if(abs(m) >=1000)  ids[i++] = dig1000(abs(m));
        if(abs(m) >=100)   ids[i++] = dig100(abs(m));
        if(abs(m) >=10)    ids[i++] = dig10(abs(m));
        ids[i++] = dig1(abs(m));
	ids[i++] = ')';
        while(i<=15)       ids[i++] = 0;

	double cl = pfl.max()*WGMp_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;
}

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

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

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

/* WGM field profile plots, image of component cp (static, t=0)
   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 WGMode::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, "WGM profile");
	mlout_gengeoxz(dat, cav, 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;
   field animations ignore the attenuation */
void WGMode::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, val_lambda(omega));
	mlout_gengeoxz(dat, cav, 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) = cav.n(x, z); x += dx; }
		z += dz;
	}
	mlout_n(dat, npx, npz, n);
	mlout_fldviewer(dat, name);
	fclose(dat);
	return;
}



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

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

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

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

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

}

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

}

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

/* add a mode */
void WGModeArray::add(WGMode m)
{
	WGMode *avec;
	avec = new WGMode[num+1];
	WGMode* ap = avec;
	WGMode* 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 WGModeArray::remove(int i)
{
	if(i<0 || i>=num) 
		wgmerror("WGModeArray: remove: argument out of range");
	if(num == 1)
	{
		delete[] mvec;
		mvec = NULL;
		num = 0;
		return;
	}
	WGMode *avec;
	avec = new WGMode[num-1];
	WGMode* ap = avec;
	WGMode* 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 WGModeArray nma */ 
void WGModeArray::merge(WGModeArray& ma)
{
	WGMode *avec;
	avec = new WGMode[num+ma.num];
	WGMode* ap = avec;
	WGMode* mp = mvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	WGMode* 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 WGMode& a, const WGMode& b)
{	
	if(a.alpha < 1.0e-14 && b.alpha < 1.0e-14)
	{
		if(a.ord < b.ord) return -1;
		if(a.ord > b.ord) return 1;
		if(a.m < b.m) return -1;
		return 1;
	}
	if(a.alpha < b.alpha) return -1;
	return 1;
}


/* sort the array */
void WGModeArray::sort()
{
	int j, k, mini;
	WGMode 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 order m, with closeby frequencies or with identical radial order */
void WGModeArray::prune()
{
        int done = 0;
	while(done == 0 && num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=num-1; ++j)
		{
			Complex to = mvec[j].c_omega();
			int mj = mvec[j].m;
			for(int k=0; k<=num-1; ++k)
			{
				Complex no = mvec[k].c_omega();
				int mk = mvec[k].m;
				if(j!=k && done == 1 && mj==mk && (to-no).abs() < WGMsP_MinOmDiff)
				{
					if(mvec[j].tr_error <= mvec[k].tr_error) rmm = k;
					else                                     rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) remove(rmm);
	}
        done = 0;
	while(done == 0 && num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=num-1; ++j)
		{
			int to = mvec[j].ord;
			int mj = mvec[j].m;
			for(int k=0; k<=num-1; ++k)
			{
				int no = mvec[k].ord;
				int mk = mvec[k].m;
				if(j!=k && done == 1 && to >= 0 && mj==mk && to == no)
				{
					if(mvec[j].tr_error <= mvec[k].tr_error) rmm = k;
					else                                     rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) remove(rmm);
	}
	return;
}

/* prune the array: among modes with identical radial order 
   select the one closest to resonance wavelength twl */
void WGModeArray::prune(double twl)
{
	int done = 0;
	while(done == 0 && num>=2)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=num-1; ++j)
		{
			int to = mvec[j].ord;
			for(int k=0; k<=num-1; ++k)
			{
				int no = mvec[k].ord;
				if(j!=k && done == 1 && to >= 0 && to == no)
				{
					if(fabs(mvec[j].lambda_r() - twl) 
					   <= fabs(mvec[k].lambda_r() - twl)) rmm = k;
					else                                  rmm = j;
					done = 0;
				}
			}
		}
		if(done == 0) remove(rmm);
	}
	return;
}

/* discard all modes with resonant wavelengths outside interval twi */
void WGModeArray::prune(Interval twi)
{
	int done = 0;
	while(done == 0 && num>=1)
	{
		done = 1;
		int rmm = 0;
		for(int j=0; j<=num-1; ++j)
		{
			double wl = mvec[j].lambda_r();
			if(wl < twi.x0 || wl > twi.x1) { rmm = j; done = 0; }
		}
		if(done == 0) remove(rmm);
	}
	return;
}

/* ------------------------------------------------------------------------
   WGM 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    WGMsP_NumSlices = 100;
// radial distance of interior boundary to pos. R, in vacuum wavelengths
double WGMsP_IntBd = 2.5;
// radial distance of exterior boundary to pos. R, in vacuum wavelengths
double WGMsP_ExtBd = 2.5;
// number of spectral terms examined
int WGMsP_NSpT = 99;
// minimum effective index (real part) considered
double WGMsP_MinNeff = 0.1;
// relative extension of effective index range towards lower values
double WGMsP_NRangeExt = 0.5;
// acceptable relative attenuation for prospective confined modes 
double WGMsP_RelAttLim = 0.1;
// range of angular mode indices around each trial number  
int WGMsP_mRange = 1;

/* complex root finding */

// initial wavelength interval, half-width in mum */
double WGMsP_IniDLambda = 0.01;
// initial attenuation, -log10(upper level) */
double WGMsP_IniMaxAtt = 2;
// initial attenuation, -log10(lower level) */
double WGMsP_IniMinAtt = 14;
// secant method, terminal distance of root prototypes
double WGMsP_RootAcc = 1.0e-12;
// number of trials for convergence with modified initial conditions
int WGMsP_LimNumIts = 30;
// maximum growth of trial value distance permitted per iteration step
double WGMsP_ConvGrowth = 2.0;
// mimimum distance in abs. frequency for resonances to be considered different
double WGMsP_MinOmDiff = 1.0e-8;
// maximum acceptable level of absolute TRC violation 
double WGMsP_MaxTRC = 1.0e-3;


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

/* for the cavity with radial layering w, radius r, centered at (xc, zc),  
   try to find pol-polarized WGMs of low radial order, 
   with eigenfrequencies next to val_omega(twl), for the target wavelength twl, 
   consider ntre * ntim trial values for complex eigenfrequencies,
   lorm=1: consider modes of lowest radial order only,  
   quiet = 0, 1, 2: all, less, no  log output  */
int wgmsolve(Waveguide w, double r, double xc, double zc, Polarization pol,
             double twl, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet)
{	
	Interval twi;
	twi.x0 = 0.0;
	twi.x1 = 1.0;
	return wgmsolve('t', w, r, xc, zc, pol, twl, twi, lorm, ntre, ntim, wgma, quiet);
}

/* for the cavity with radial layering w, radius r, centered at (xc, zc),  
   try to find pol-polarized WGMs of low radial order, 
   with eigenfrequencies corresponding to wavelengths in the interval twi, 
   consider ntre * ntim trial values at central resonance wavelength,  
   lorm=1: consider modes of lowest radial order only,  
   quiet = 0, 1, 2: all, less, no  log output  */
int wgmsolve(Waveguide w, double r, double xc, double zc, Polarization pol,
             Interval twi, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet)
{	
	double twl = 1.0;
	return wgmsolve('i', w, r, xc, zc, pol, twl, twi, lorm, ntre, ntim, wgma, quiet);
}

/* ... generic: 
   stype = 't': search resonances at target wavelength twl,
   stype = 'i': find resonances within wavelength interval twi */
int wgmsolve(char stype, Waveguide w, double r, double xc, double zc, Polarization pol,
             double twl, Interval twi, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet)
{	
	if(quiet < 2)
	{
fprintf(stderr, "\n------------- WGMsolve --------------------------------------------- '18 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
fprintf(stderr, "R: %g   (x0, z0): (%g, %g)\n", r, xc, zc);
if(stype == 't') fprintf(stderr, "target lambda_r: %g\n", twl);
else             fprintf(stderr, "target lambda_r: [%g, %g]\n", twi.x0, twi.x1);
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");
	}
	if(stype == 't') twi = Interval(twl-1.0, twl + 1.0);
	else             twl = 0.5*(twi.x0+twi.x1);

	Dvector tbv;

	if(ntre < 0)
	{

		ModeArray ma, mas;
		mas.clear();

		double ib, eb, ai, ae;
		ib = r+w.hx(0)-WGMsP_IntBd*twl; if(ib < r/10.0) ib = r/10.0; 
		eb = r+w.hx(w.nx)+WGMsP_ExtBd*twl; 
		ai = r*log(ib/r);
		ae = r*log(eb/r);
		int nx = WGMsP_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 = twl;
		
		modespectrum(eqwg, pol, Interval(ai, ae), LIMD, WGMsP_NSpT, ma, 1);
		double wnm;
		if(lorm == 1) wnm = w.n.min()*val_k0(twl);
		else          wnm = WGMsP_MinNeff*val_k0(twl);
		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();

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(twl);
                }
                else
                {
                        mib = WGMsP_MinNeff*val_k0(twl);
                }
		double mab = w.n.max()*val_k0(twl); 
		tbv = Dvector(ntre);
		for(int i=0; i<=ntre-1; ++i) 
			tbv(i) = mib + (mab-mib)/((double) (ntre-1))*((double) i);
	}

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

	if(ntim <= 0) ntim = 10;
	double astep = (WGMsP_IniMinAtt-WGMsP_IniMaxAtt)/((double) ntim);

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

	Dvector pom; pom.strip();
	for(int j=0; j<=tbv.nel-1; ++j)  
	{
		double neff = tbv(j)/val_k0(twl);
		double m_d = 2.0*PI/twl*r*neff;
		pom.append(m_d);
	}	
	for(int j=0; j<=pom.nel-1; ++j)  pom(j) = floor(pom(j)+0.5);
	int nm=pom.nel;
	for(int j=0; j<=nm-1; ++j)
	{
		int m = (int) pom(j);
		for(int i=1; i<=WGMsP_mRange; ++i)
		{
			if(m-i>=0) pom.append((double) (m-i));
			pom.append(m+i);
		}
	}
	pom.sort();
	int done = 0;
	while(done==0)
	{
		done = 1;
		int rmi=-1; 
		for(int i=0; i<=pom.nel-2; ++i)
			if(fabs(pom(i)-pom(i+1))<0.1) { rmi = i; done = 0; }
		if(done == 0) pom.remove(rmi);
	}
	
if(quiet < 1)
{
	fprintf(stderr, "-> prosp. m =");
	for(int j=0; j<=pom.nel-1; ++j)  
	{
		fprintf(stderr, " %g", pom(j));
	}
	fprintf(stderr, ".\n");
	fprintf(stderr, "----------------------------------------------------------------------------\n");
}
	wgma.clear();

	WGMode wgm;
	wgm.pol = pol;
	wgm.R = r;
	wgm.x0 = xc;
	wgm.z0 = zc;
	wgm.wg = w;
	wgm.cav = OvlStruct(w.n(w.nx+1), twl); 
	wgm.cav.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);
			wgm.cav.place(o);
		}
	}	
	wgm.dir = FORW;
	wgm.cRi = (wgm.wg.hx(0)+wgm.R)-WGMp_R_CROP_IN*twl; 
	if(wgm.cRi < wgm.R*WGMp_R_TINY) wgm.cRi = wgm.R*WGMp_R_TINY;
	wgm.cRe = (wgm.wg.hx(wgm.wg.nx)+wgm.R)+WGMp_R_CROP_OUT*twl; 

	wgm.cA = Cvector(w.nx+2);
	wgm.cB = Cvector(w.nx+2);

	for(int mi=0; mi<=pom.nel-1; ++mi)
	{
		wgm.m = (int) floor(pom(mi)+0.5);
		for(int oi=0; oi<=1; ++oi)
		{
			for(double ai=WGMsP_IniMaxAtt; ai<=WGMsP_IniMinAtt; ai += astep)
			{
				double maxatt = pow(10.0, -ai);
				double minatt = pow(10.0, -(ai+astep));

				Complex z1, z2;

				if(oi==0)
				{
					z1 = Complex(val_omega(twl-WGMsP_IniDLambda), maxatt); 
					z2 = Complex(val_omega(twl+WGMsP_IniDLambda), minatt); 
				}
				else 
				{
					z1 = Complex(val_omega(twl+WGMsP_IniDLambda), maxatt); 
					z2 = Complex(val_omega(twl-WGMsP_IniDLambda), minatt); 
				}

				int succ = wgm.converge(z1, z2, quiet);
				if(succ == 1) wgma.add(wgm);
			}
		}
	}
if(quiet < 1)
{
	fprintf(stderr, "----------------------------------------------------------------------------\n");
}
	wgma.prune();

	if(stype == 't')
	{
		wgma.prune(twl);
		wgma.sort();

if(quiet < 2)
{
		for(int i=0; i<=wgma.num-1; ++i)
		{
			fprintf(stderr, "[%d] %s, om=%g+i%g, lambda_r=%g, FWHM_lambda=%g, Q=%g\n", 
			i, wgma(i).ids, wgma(i).c_omega().re, wgma(i).c_omega().im, wgma(i).lambda_r(), wgma(i).fwhm_lambda(), wgma(i).Q());
		}
		fprintf(stderr, "Ok.\n\n");
}
		return wgma.num;
	}
	WGModeArray twgma;
	twgma = wgma;
	twgma.prune(twl);
	for(int j=0; j<=twgma.num-1; ++j)
	{
		double wl, dwl, att; 
		int mt;
		wl = val_lambda(twgma(j).omega);
		att = twgma(j).alpha;
		mt = twgma(j).m;
		dwl = wl/((double) mt);
		wl -= dwl;
		++mt;
		while(wl >= twi.x0)
		{
			wgm.m = mt;
			Complex z1 = Complex(val_omega(wl-WGMsP_IniDLambda), att); 
			Complex z2 = Complex(val_omega(wl+WGMsP_IniDLambda), att); 
			int succ = wgm.converge(z1, z2, quiet);
			if(succ == 1) wgma.add(wgm);
			wl -= dwl;
			++mt;
		}
		wl = val_lambda(twgma(j).omega);
		mt = twgma(j).m;
		wl += dwl;
		--mt;
		while(wl <= twi.x1)
		{
			wgm.m = mt;
			Complex z1 = Complex(val_omega(wl-WGMsP_IniDLambda), att); 
			Complex z2 = Complex(val_omega(wl+WGMsP_IniDLambda), att); 
			int succ = wgm.converge(z1, z2, quiet);
			if(succ == 1) wgma.add(wgm);
			wl += dwl;
			--mt;
		}
	}
if(quiet < 1)
{
	fprintf(stderr, "----------------------------------------------------------------------------\n");
}
	wgma.prune();
	wgma.prune(twi);
	wgma.sort();

if(quiet < 2)
{
	for(int i=0; i<=wgma.num-1; ++i)
	{
		fprintf(stderr, "[%d] %s, om=%g+i%g, lambda_r=%g, FWHM_lambda=%g, Q=%g\n", 
		i, wgma(i).ids, wgma(i).c_omega().re, wgma(i).c_omega().im, wgma(i).lambda_r(), wgma(i).fwhm_lambda(), wgma(i).Q());
	}
	fprintf(stderr, "Ok.\n\n");
}
	return wgma.num;
}