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

/*
 * slams.cpp
 * Mode solver routines,
 * guided mode analysis & rigorous spectral dicretization
 */

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

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

/* --- mode analysis, "tuning" parameters -------------------------------- */

/* tolerance for fixing the propagation constants */
double SLAMS_BQTOL = 1.0e-13;

/* absolute tolerance */ 
double bqTol(double bq)
{
	double abq = fabs(bq);
	if(abq < 10.0) return SLAMS_BQTOL;
	return abq*SLAMS_BQTOL;
}

/* minimum distance at which propagation constants are considered different */
double SLAMS_BQSEP = 1.0e-8;

/* absolute separation limit */ 
double bqSep(double bq)
{
	double abq = fabs(bq);
	if(abq < 10.0) return SLAMS_BQSEP;
	return abq*SLAMS_BQSEP;
}

/* tolerance for the maximum relative error in the continuity of the 
   basic field component for accepting a field as a valid mode */
double SLAMS_DISCERRTOL = 1.0e-6;

/* mode spectrum computation, BEP/QUEP field discretization:
   minimum difference between a (vaguely defined) "mode propagation angle"
   and the transverse axis. If modes close to the transition from propagating
   to evanescent type with almost vanishing propagation constant occur,
   the spectral discretization procedure issues a warning message
   ((!) in quiet mode), since this can lead to highly erroneous results.
   Workaround: a slight change of the computational window width */
double SLAMS_MODEANGLELIM = 1.0e-3;      // approx 0.05*PI/180.0

/* --- mode analysis, helper functions ----------------------------------- */

/* insert a new mode in the list of found ones, return one, if added */
int recognizemode(Mode& m, int lop, Ivector& found, Dvector& bqf, 
                  ModeArray& ma, int verb)
{
	if(m.special == 0)
	{
		if(m.ord > found.nel-1) 
		{
			fprintf(stderr, " %s, bq = %.15g ?\n", m.ids, m.betaq);
			m.plot('x', 'x');
			slamserror( "recognizemode: high order mode");
		}
		if(found(m.ord) == 1)
		{
			if(fabs(m.betaq-bqf(m.ord)) > bqSep(m.betaq))
			{
				fprintf(stderr, "\nFound: \n");
				for(int i=lop; i<=found.nel-1; ++i)
				if(found(i) != 0 && abs(i-m.ord)<5)
				fprintf(stderr, "bq_%d = %.15g\n", i, bqf(i));
				fprintf(stderr, " %s, bq = %.15g ?\n",
						m.ids, m.betaq);
				m.plot('x', 'x');
				slamserror( "recognizemode: classification");
			}
			return 0;
		}	
		if(found(m.ord) == 2)
		{
			double tbq = bqf(m.ord);
			int i = lop; 
			while(i<=ma.num-1 
			      && fabs(ma(i).betaq - tbq) > bqSep(tbq)) ++i;
			if(i >= ma.num) 
				slamserror("recognizemode: confused, 0");
			Mode tm = ma(i);
			ma.remove(i);
			ma.add(m);
			found(m.ord) = 1;
			bqf(m.ord) = m.betaq;
			if(verb == 2) fprintf(stderr, "%s ", m.ids);
			if(verb == 1) fprintf(stderr, ".");
			recognizemode(tm, lop, found, bqf, ma, verb); 
			return 1;
		}
		ma.add(m);	
		if(verb == 2) fprintf(stderr, "%s ", m.ids);
		if(verb == 1) fprintf(stderr, ".");
		bqf(m.ord) = m.betaq;
		found(m.ord) = 1;
		return 1;
	}
	// treat special mode
	for(int i=lop; i<=found.nel-1; ++i)
	{
		if(found(i) != 0 && fabs(m.betaq-bqf(i)) < bqSep(bqf(i))) 
			return 0;
	}
	int p = lop;
	while(p <= found.nel-1)
	{
		if(found(p) == 0)
		{
			int acc = 1;
			int ps = p-1;
			while(ps>=lop && acc==1)
			{ if(found(ps)==1 && bqf(ps)<m.betaq) acc = 0; --ps; }
			ps = p+1;
			while(ps<=found.nel-1 && acc==1)
			{ if(found(ps)==1 && bqf(ps)>m.betaq) acc = 0; ++ps; }
			if(acc) 
			{
				ma.add(m);
				if(verb == 2) fprintf(stderr, "%s ", m.ids);
				if(verb == 1) fprintf(stderr, "x");
				bqf(p) = m.betaq;
				found(p) = 2;
				return 1;
			}
		} 
		if(found(p) == 2)
		{
			int acc = 1;
			int ps = p-1;
			while(ps>=lop && acc==1)
			{ if(found(ps)==1 && bqf(ps)<m.betaq) acc = 0; --ps; }
			ps = p+1;
			while(ps<=found.nel-1 && acc==1)
			{ if(found(ps)==1 && bqf(ps)>m.betaq) acc = 0; ++ps; }
			if(acc == 1 && bqf(p) < m.betaq)
			{
				double tbq = bqf(p);
				int i = lop; 
				while(i<=ma.num-1 
				      && fabs(ma(i).betaq - tbq) > bqSep(tbq))
					 ++i;
				if(i >= ma.num) 
				      slamserror("recognizemode: confused, 1");
				Mode tm = ma(i);
				ma.remove(i);
				ma.add(m);
				found(p) = 2;
				bqf(p) = m.betaq;
				if(verb == 2) fprintf(stderr, "%s ", m.ids);
				if(verb == 1) fprintf(stderr, "x");
				recognizemode(tm, lop, found, bqf, ma, verb); 
				return 1;
			}
		}
		 ++p;
	}
	return 0;
}

/* trap propagation constants */
#define MAXNUMSCF 100   // record this many suspicious configurations for the
                        // further propagation constant search       
int trap(ModeArray& ma, Mode& m, 
         double bqmin, double bqmax, double sf, int sm, 
         Dvector& bqf, Ivector& found, Dmatrix& susp, int& nsc, int lfsc, 
         int verb)
{
	int newv;
	double step;
	double bq0, bq1, bq2, f0, f1, f2, mbq; 
	double err;

	if(bqmax <= bqmin) return 0;
	if(bqmax-bqmin < bqSep(bqmax)) return 0;
	bq0 = bqmax-bqTol(bqmax);
	int zmin = m.nummodesabove(bq0);
	bq0 = bqmin+bqTol(bqmin);
	int zmax = m.nummodesabove(bq0);
	if(zmin >= zmax) return 0;

	newv = 0;
	double nst = 10.0*pow(10.0, sf);
	step = (bqmax-bqmin)/nst;
	if(step < bqTol(bqmax)) step = bqTol(bqmax); 
	if(step < bqTol(bqmin)) step = bqTol(bqmin); 

// fprintf(stderr, "trap(%.15g#%d <--T%c-- %.15g#%d <%.15g>, %d)\n", 
//                 bqmin, zmax, polCHR(m.pol), bqmax, zmin, step, sm);

	bq1 = bqmax;
	f1 = m.travres(bq1);
	bq0 = bqmax;
	f0 = f1; 
	for(double bq=bqmax-step; bq>=bqmin; bq-=step)
	{
		bq2 = bq;
		f2 = m.travres(bq2);
		if(f1*f2 <= 0.0) 
		{
			mbq = m.localize(bq2, bq1, bqTol(bq2));	
			m.polish(mbq);
			err = m.checkcontinuity();
// fprintf(stderr, "\n< %s >\n", m.ids);
// fprintf(stderr, "n_eff = %.16g\n", m.neff);
// fprintf(stderr, "betaq = %.16g\n", m.betaq);
// fprintf(stderr, "spec  = %d\n",    m.special);
// fprintf(stderr, "disc  = %.16g\n", m.checkcontinuity());
// fprintf(stderr, "dtol  = %.16g\n", SLAMS_DISCERRTOL);
			if(err < SLAMS_DISCERRTOL)
			{
				newv = recognizemode(m, zmin, found, bqf, 
				                     ma, verb);
				if(newv != 0 && sm == 1) return newv;
			}
			// else { m.plot('d', 'd'); }
		}
		if(lfsc != 0)
		{
			if(f1 > 0.0 && f0 > f1 && f1 < f2) 
			{
				if(nsc < MAXNUMSCF-1)
				{
					susp(nsc, 0) = bq2;
					susp(nsc, 1) = bq0;
					++nsc;
				}
			}
			if(f1 < 0.0 && f0 < f1 && f1 > f2) 
			{
				if(nsc < MAXNUMSCF-1)
				{
					susp(nsc, 0) = bq2;
					susp(nsc, 1) = bq0;
					++nsc;
				}
			}
		}
		bq0 = bq1; f0 = f1;
		bq1 = bq2; f1 = f2;
	}
	return newv;
}

int trap(ModeArray& ma, Mode& m, 
         double bqmin, double bqmax, double sf, int sm, int lfsc, 
         Dvector& bqf, Ivector& found, int verb)
{
	int newv, nsc = 0;
	Dmatrix susp(MAXNUMSCF, 2);
	newv = trap(ma, m, bqmin, bqmax, sf, sm, bqf, found, 
	            susp, nsc, 1, verb);
	if(newv > 0) return newv;
	if(lfsc == 0) return newv;
	double minbq, maxbq;
	while(nsc >= 1)
	{
		minbq = susp(nsc-1, 0);
		maxbq = susp(nsc-1, 1);
		--nsc;
		newv += trap(ma, m, minbq, maxbq, sf, sm, bqf, found, 
		             susp, nsc, 0, verb);
		if(sm == 1 && newv > 0) return newv;
	}
	return newv;
}


/* --- mode analysis --------------------------------------------------- */

/* determine all modes with squared propagation constants between 
   bqmin and bqmax, for the problem specified by wg, pol, 
   boundaries of type bdtb at bpb, of type bdtt at bpt; 
   verb: 0: no control output, 1: progress, 2: all, 
   nstd: 1: solve nonstandard  effective VEIMS problem,
   ec: VEIMS effective permittivity */
void findmodes(Waveguide wg, Polarization pol, 
               Boundary_type bdtb, double bpb,
               Boundary_type bdtt, double bpt,
               double bqmin, double bqmax, ModeArray& ma, int verb, 
               int nstd, Dvector ec)
{
	double k0 = val_k0(wg.lambda);
	double tbq, elim; 
	int cl, split, nsm;

// fprintf(stderr, "\nFM_I( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
// fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//                  bdtb, bpb, bdtt, bpt);
// fprintf(stderr, "      bqmin = %g, bqmax = %g )\n", bqmin, bqmax);
// verb = 2; 
// wg.write(stderr);
	
	if(bqmax < bqmin+bqSep(bqmax)) return;	
//	if(bqmax < bqmin+bqSep(bqmax)) slamserror("findmodes: bqmin, bqmax");	
	cl = wg.checksymmetry();
	split = 1;
	if(cl <= 1 || cl >= wg.nx) split = 0;
	if(bdtb != bdtt) split = 0;
	if(bdtb != OPEN && fabs((wg.hx(0)-bpb)-(bpt-wg.hx(wg.nx)))>COMPTOL_HX) 
		split = 0;
	if(split == 1)
	{
		Waveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		ModeArray sm; 
		double cp = (wg.hx(cl-1)+wg.hx(cl))/2.0;
		if(verb == 1) fprintf(stderr, "s");
		if(verb == 2) fprintf(stderr, "sym ");
		sm.clear();
		findmodes(wg0, pol, bdtb, bpb, LIMN, cp, 
		          bqmin, bqmax, sm, verb, nstd, ec); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				;
			else sm(j).expand(); 
		}
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "a");
		if(verb == 2) fprintf(stderr, "asy ");
		sm.clear();
		findmodes(wg0, pol, bdtb, bpb, LIMD, cp, 
		          bqmin, bqmax, sm, verb, nstd, ec); 
		nsm = sm.num;
		for(int j=0; j<=nsm-1; ++j) 
		{		
			if(sm(j).special==2 && sm(j).bdt_t==OPEN) 
				sm(j).mirror(cp);
			else sm(j).expand(); 
		}
		ma.merge(sm);
		ma.sort();
		return;
	}	
	elim = MAXREFIND*MAXREFIND; 
	cl = -1;
	for(int l=1; l<=wg.nx; ++l)
	{ 
		double te = wg.decoupepseff(l);
		if(te < elim) 
		{
			elim = te;
			cl = l;
		}
	}
	tbq = k0*k0*elim;
	if(tbq < bqmin) tbq = bqmin;
	if(bqmax > tbq+bqSep(tbq))
	{
		Waveguide wg0, wg1;
		wg.split(cl, wg0, wg1);
		ModeArray sm;
		if(verb == 1) fprintf(stderr, "d-");
		if(verb == 2) fprintf(stderr, "d^(0, %d) ", cl);
		sm.clear();
		findmodes(wg0, pol, bdtb, bpb, OPEN, AWAY, 
		          tbq, bqmax, sm, verb, nstd, ec); 
		for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
		ma.merge(sm);
		if(verb == 1) fprintf(stderr, "d+");
		if(verb == 2) fprintf(stderr, "d^(%d, %d) ", cl, wg.nx+1);
		sm.clear();
		findmodes(wg1, pol, OPEN, -AWAY, bdtt, bpt, 
		          tbq, bqmax, sm, verb, nstd, ec); 
		for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
		ma.merge(sm);
		if(tbq > bqmin+bqSep(bqmin))
		{
			if(verb == 1) fprintf(stderr, "dr");
			if(verb == 2) fprintf(stderr, "d_(0, %d) ", wg.nx+1);
			sm.clear();
			findmodes(wg, pol, bdtb, bpb, bdtt, bpt, 
				  bqmin, tbq, sm, verb, nstd, ec); 
			ma.merge(sm);
		}
		ma.sort();
		return;
	}
	if(bdtb != OPEN)
	{
		elim = wg.lbdecepseff(bpb); 
		tbq = k0*k0*elim;
		if(tbq < bqmin) tbq = bqmin;
		if(bqmax > tbq+bqSep(tbq))
		{
			ModeArray sm;
			if(verb == 1) fprintf(stderr, "l^");
			if(verb == 2) fprintf(stderr, "l^ ");
			sm.clear();
			findmodes(wg, pol, OPEN, -AWAY, bdtt, bpt, 
				  tbq, bqmax, sm, verb, nstd, ec); 
			for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
			ma.merge(sm);
			if(bqmin < tbq-bqSep(tbq))
			{
				if(verb == 1) fprintf(stderr, "l_");
				if(verb == 2) fprintf(stderr, "l_ ");
				sm.clear();
				findmodes(wg, pol, bdtb, bpb, bdtt, bpt, 
					  bqmin, tbq, sm, verb, nstd, ec); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt != OPEN)
	{
		elim = wg.ubdecepseff(bpt); 
		tbq = k0*k0*elim;
		if(tbq < bqmin) tbq = bqmin;
		if(bqmax > tbq+bqSep(tbq))
		{
			ModeArray sm;
			if(verb == 1) fprintf(stderr, "u^");
			if(verb == 2) fprintf(stderr, "u^ ");
			sm.clear();
			findmodes(wg, pol, bdtb, bpb, OPEN, AWAY, 
				  tbq, bqmax, sm, verb, nstd, ec); 
			for(int j=0; j<=sm.num-1; ++j) sm(j).special = 2;
			ma.merge(sm);
			if(bqmin < tbq-bqSep(tbq))
			{
				if(verb == 1) fprintf(stderr, "u_");
				if(verb == 2) fprintf(stderr, "u_ ");
				sm.clear();
				findmodes(wg, pol, bdtb, bpb, bdtt, bpt, 
					  bqmin, tbq, sm, verb, nstd, ec); 
				ma.merge(sm);
			}
			ma.sort();
			return;
		}
	}
	if(bdtt==OPEN && bdtb==OPEN)
	{
		double emin = wg.defaultepseffmin();	
		double emax = wg.defaultepseffmax();	
		tbq = k0*k0*emin;
		if(bqmin < tbq) bqmin = tbq;
		tbq = k0*k0*emax;
		if(bqmax > tbq) bqmax = tbq;
		if(bqmax < bqmin+bqSep(bqmin)) return;
	}
	double bq, dbq;
	int zmax, zmin;
	Mode m;
	if(nstd != 1) m = Mode(wg, bdtb, bpb, bdtt, bpt, pol);
	else          m = Mode(wg, ec, bdtb, bpb, bdtt, bpt);
	bq = bqmax-bqTol(bqmax);
	zmin = m.nummodesabove(bq);
	bq = bqmin+bqTol(bqmin);
	zmax = m.nummodesabove(bq);
	if(zmin >= zmax) return;

// fprintf(stderr, "\nFM_F( wg.nx = %d, pol = T%c\n", wg.nx, polCHR(pol));
// fprintf(stderr, "      bdtb = %d, bpb = %g, bdtt = %d, bpt = %g\n", 
//                 bdtb, bpb, bdtt, bpt);
// fprintf(stderr, "      bqmin = %.15g, zmax = %d, bqmax = %g, zmin= %d)\n", 
//                 bqmin, zmax, bqmax, zmin);
// m.travresinspect(bqmin+bqTol(bqmin), bqmax-bqTol(bqmax));

	Ivector found(zmax);
	found.init(0);	
	Dvector bqf(zmax);
	bqf.init(0.0);

	trap(ma, m, bqmin, bqmax, 2.0, 0, 1, bqf, found, verb);
	
	int fa, mn, i, newv;
	double bqup, bqdn, sf;
	fa = 0;
	while(fa == 0)
	{
		fa = 1; 
		mn = zmax-1;
		while(fa == 1 && mn>= zmin)
		{
			if(found(mn) != 0) --mn;
			else
			{
				fa = 0;
				i=mn+1; 
				while(i<zmax && found(i) == 0) ++i;
				if(i >= zmax) bqdn = bqmin;
				else          bqdn = bqf(i);
				i=mn-1;
				while(i>zmin-1 && found(i) == 0) --i;
				if(i <= zmin-1) bqup = bqmax;
				else 	  bqup = bqf(i);
				dbq = bqup-bqdn;
				sf = 1.0;
				while((newv = trap(ma, m, 
					      bqdn+bqTol(bqdn), 
				              bqup-bqTol(bqup), sf, 1, 1,
					      bqf, found, verb)) == 0 
				      && sf < 4.0) sf += 0.5;
				if(newv >= 1) goto fnd;
				if(mn < zmax-1 && found(mn+1) != 0)
				{
				     bqdn = bqf(mn+1);
				     sf = 1.0;
				     tbq = bqdn+dbq/100.0;
				     while((newv = trap(ma, m, bqdn, tbq, sf, 
				                        1, 1, bqf, found, verb)
					    ) == 0 && sf < 4.0) sf += 0.5;
				}
				if(newv >= 1) goto fnd;
				if(mn > zmin && found(mn-1) != 0)
				{
				     bqup = bqf(mn-1);
				     sf = 1.0;
				     tbq = bqup-dbq/100.0;
				     while((newv = trap(ma, m, tbq, bqup, sf,
				                        1, 1, bqf, found, verb)
					    ) == 0 && sf < 4.0) sf += 0.5;
				}
				if(newv >= 1) goto fnd;

				if(mn == zmax-1)
				{
				     bqdn = bqmin;
				     sf = 1.0;
				     tbq = bqdn+dbq/100.0;
				     while((newv = trap(ma, m, bqdn, tbq, sf, 
				                        1, 1, bqf, found, verb)
					    ) == 0 && sf < 4.0) sf += 0.5;
				}
				if(newv >= 1) goto fnd;
				if(mn == zmin)
				{
				     bqup = bqmax;
				     sf = 1.0;
				     tbq = bqup-dbq/100.0;
				     while((newv = trap(ma, m, tbq, bqup, sf,
				                        1, 1, bqf, found, verb)
					    ) == 0 && sf < 4.0) sf += 0.5;
				}
				if(newv >= 1) goto fnd;

				fprintf(stderr, " T%c%d ?\n", 
				                polCHR(m.pol), mn);
				m.travresinspect(bqmin, bqmax);
				m.wg.write(stderr);
				m.wg.plot('e','r');
				slamserror("findmodes: missing mode");
fnd: 					;;
			}
		}
	}
	ma.sort();
	return;
}
void findmodes(Waveguide wg, Polarization pol, 
               Boundary_type bdtb, double bpb,
               Boundary_type bdtt, double bpt,
               double bqmin, double bqmax, ModeArray& ma, int verb)
{
	Dvector dummy;
	findmodes(wg, pol, bdtb, bpb, bdtt, bpt, bqmin, bqmax, ma, verb, 
	          0, dummy);
	return;
}

/* homogeneous medium, constant refractive index: 
   determine mnm modes up to order mnm-1, for the problem specified by 
   wg, pol, with boundaries of type bdt at bpb and bpt; 
   verb: 0: no control output, 1: progress, 2: all */
void freespacemodes(Waveguide wg, Polarization pol, 
                    Boundary_type bdt, double bpb, double bpt,
                    int mnm, ModeArray& ma, int verb)
{
	double k0 = val_k0(wg.lambda);
	double e = wg.eps(0);
	if(wg.constmed() != 1) 
			slamserror("freespacemodes: eps not constant");
	if(fabs(bpt-bpb)< LIMBDSEP) 
			slamserror("freespacemodes: bpb, bpt");
	if(bdt != LIMD && bdt != LIMN) 
			slamserror("freespacemodes: bdt");
	if(verb == 2) fprintf(stderr, "fs ");
	if(verb == 1) fprintf(stderr, "fs");
	ma.clear();
	Mode m(wg, bdt, bpb, bdt, bpt, pol);
	double gamma, bq;
	for(int z=1; z<=mnm-1; ++z)
	{
		gamma = ((double) z)*PI/(bpt-bpb);
		bq = k0*k0*e-gamma*gamma;
		m.polish(bq);	
		ma.add(m);
		if(verb == 2) fprintf(stderr, "%s ", m.ids);
		if(verb == 1) fprintf(stderr, ".");
	}
	switch(bdt)
	{
	case LIMD:
		gamma = ((double) mnm)*PI/(bpt-bpb);
		bq = k0*k0*e-gamma*gamma;
		m.polish(bq);	
		ma.add(m);
		if(verb == 2) fprintf(stderr, "%s ", m.ids);
		if(verb == 1) fprintf(stderr, ".");
		break;
	case LIMN:
		bq = k0*k0*e;
		m.polish(bq);	
		ma.add(m);
		if(verb == 2) fprintf(stderr, "%s ", m.ids);
		if(verb == 1) fprintf(stderr, ".");
		break;
	default:
		slamserror("freespacemodes: ?");
		break;
	}
	ma.sort();
	return;	
}

/* --- mode analysis, simpler driver routines -------------------------- */

/* complete guided mode analysis, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   ma: the modes found (output)
   quiet == 1: suppress log output */
int modeanalysis(Waveguide wg, Polarization pol, ModeArray& ma, int quiet)         
{
	wg.consistency();
	wg.optimize();
	double emin = wg.defaultepseffmin();	
	double emax = wg.defaultepseffmax();	
	double k0 = val_k0(wg.lambda);
	double bqmin = k0*k0*emin;
	double bqmax = k0*k0*emax;

	if(quiet != 1)
	{
fprintf(stderr, "\n------------- Slams - MA ------------------------------------------- '04 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
if(wg.special) fprintf(stderr, "(!) eps_eff: (%.10g, %.10g)  ", emin, emax);
else fprintf(stderr, "n_eff: (%.10g, %.10g)  ", sqrt(emin), sqrt(emax));
fprintf(stderr, "Tol: %g  ", SLAMS_BQTOL);
fprintf(stderr, "Sep: %g\n", SLAMS_BQSEP);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", wg.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=wg.nx; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
fprintf(stderr, "\n");
if(wg.special) fprintf(stderr, "   eps: ");
else           fprintf(stderr, "     n: ");
for(int j=0; j<=wg.nx+1; ++j) fprintf(stderr, "%6.10g  ", wg.n(j));
fprintf(stderr, "\n");
fprintf(stderr, "lambda: %.10g  k0: %g\n", wg.lambda, k0);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	bqmin += SLAMS_BQSEP;
	int verb = 2; if(quiet == 1) verb = 0;
	ma.clear();
	if(bqmax > bqmin+bqSep(bqmax)) 
		findmodes(wg, pol, OPEN, -AWAY, OPEN, AWAY, 
	                  bqmin, bqmax, ma, verb);
	ma.sort();
	for(int j=0; j<=ma.num-1; ++j)
	{
		if(ma(j).special==0 && ma(j).ord!=j) 
			slamserror("modeanalysis: mode set corrupted");
		ma(j).setids(j);
	}
	if(quiet != 1)
	{
		if(ma.num<=0) fprintf(stderr, "No guided modes.\n");
		else
fprintf(stderr, "\n----------------------------------------------------------------------------\n");
		for(int j=0; j<=ma.num-1; ++j)
		{
fprintf(stderr, "%s: n_eff = %.8g, beta = %.8g, npcB = %.8g",
                 		ma(j).ids, ma(j).neff, ma(j).beta, ma(j).npcB);
			if(wg.nx==1 && wg.special==0) 
fprintf(stderr, ", theta = %.4g\n", acos(ma(j).neff/wg.n(1))*180.0/PI);
			else 
fprintf(stderr, "\n");
		}
		if(ma.num == 2 && 
                   fabs(ma(0).betaq-ma(1).betaq)>bqSep(ma(1).betaq))
fprintf(stderr, "Lc = %.6g\n", PI/(ma(0).beta-ma(1).beta));
fprintf(stderr, "\n");
	}
	return ma.num;
}

int modeanalysis(Waveguide wg, Polarization pol, ModeArray& ma)         
{
	return modeanalysis(wg, pol, ma, 0);
}

/* compute the discretized mode spectrum on an interval with artificial
   boundaries, returns the number of found modes
   wg: the structure under consideration
   pol: polarization type, TE or TM
   cw: the computational window
   bt: the type of boundary, one of LIMD or LIMN
   mnm: highest order mode (+1) to be included, determine up to mnm eigenmodes 
   ma: the modes found (output)
   quiet: == 1: suppress log output */
int modespectrum(Waveguide wg, Polarization pol, Interval cw, Boundary_type bt,
                 int mnm, ModeArray& ma, int quiet)
{
	wg.consistency();
	wg.optimize();
	if(cw.x0 >= wg.hx(0)-LIMBDSEP) 
		slamserror("modespectrum: cw.x0");
	if(cw.x1 <= wg.hx(wg.nx)+LIMBDSEP) 
		slamserror("modespectrum: cw.x1");
	if(bt != LIMD && bt != LIMN) 
		slamserror("modespectrum: bdtype");
	double k0 = val_k0(wg.lambda);
	double emax;
	if(wg.special) emax = wg.n.max();
	else         { emax = wg.n.max(); emax = emax*emax; }
	double bqmax = k0*k0*emax;

	if(quiet != 1)
	{
fprintf(stderr, "\n------------- Slams - SP ------------------------------------------- '04 ---\n");
fprintf(stderr, "T%c  ", polCHR(pol));
if(wg.special) fprintf(stderr, "(!) eps_eff: (... , %.10g)  ", emax);
else fprintf(stderr, "n_eff: (... , %.10g)  ", sqrt(emax));

fprintf(stderr, "#: %d  ", mnm);
fprintf(stderr, "x in [%g, %g]_", cw.x0, cw.x1);
if(bt == LIMD) fprintf(stderr, "D  ");
else           fprintf(stderr, "N  ");
fprintf(stderr, "Tol: %g  ", SLAMS_BQTOL);
fprintf(stderr, "Sep: %g\n", SLAMS_BQSEP);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", wg.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=wg.nx; ++j) fprintf(stderr, "%6.10g  ", wg.hx(j));
fprintf(stderr, "\n");
if(wg.special) fprintf(stderr, "   eps: ");
else           fprintf(stderr, "     n: ");
for(int j=0; j<=wg.nx+1; ++j) fprintf(stderr, "%6.10g  ", wg.n(j));
fprintf(stderr, "\n");
fprintf(stderr, "lambda: %.10g  k0: %g\n", wg.lambda, k0);
fprintf(stderr, "----------------------------------------------------------------------------\n");
	}
	int verb = 1; if(quiet == 1) verb = 0;
	if(wg.constmed())
	{		
		ma.clear();
		freespacemodes(wg, pol, bt, cw.x0, cw.x1, mnm, ma, verb);
	}
	else
	{
		Mode m(wg, bt, cw.x0, bt, cw.x1, pol);
		double bq, bqmin;
		int z;
		bqmin = bqmax;
		bq = bqmax-bqTol(bqmax);
		// z = m.nummodesabove(bq);
		// if(mnm<z) return 0;
		double elow; 
		elow = wg.eps(0);
		for(int l=1; l<=wg.nx+1; ++l) 
			{ if(wg.eps(l) < elow) elow = wg.eps(l); }
		bq = m.k0*m.k0*elow;
		bq += bqTol(bq);
		z = m.nummodesabove(bq);
		if(z >= mnm) bqmin = bq;
		else 
		{	
			double bql, bqu;
			int zl, zu;
			bqu = bq; zu = z;
			bql = m.k0*m.k0*elow
			      -((double) (mnm*mnm))*PI*PI
			      /(cw.len())/(cw.len());
			zl = m.nummodesabove(bql);
			while(zl < mnm)
			{
				double d = fabs(bql)*0.01;
				if(d < 10.0) d = 10.0;
				bql -= d;
				zl = m.nummodesabove(bql);
			}
			if(zl == mnm) bqmin = bql;
			else
			{
				while(zl > zu)
				{
					if(zl < mnm || zu > mnm) 
					slamserror("modespectrum: confused");
					bq = (bqu+bql)*0.5;	
					z = m.nummodesabove(bq);
					if(z == mnm) { bqmin=bq; zl=zu=z; }
					if(z  < mnm) { zu = z; bqu = bq; }
					if(z  > mnm) { zl = z; bql = bq; }
				}	
			}
		}
		ma.clear();
		findmodes(wg, pol, bt, cw.x0, bt, cw.x1, 
		          bqmin, bqmax, ma, verb);
	}
	if(ma.num > mnm) 
	{
		ma.sort();
		while(ma.num > mnm) ma.remove(ma.num-1);
	}
	double nmax;
	if(wg.special) 
	{
		emax = wg.n.max(); if(emax<1.0) emax = 1.0; nmax = sqrt(emax);
	}
	else nmax = wg.n.max();
	for(int j=0; j<=ma.num-1; ++j)
	{
		if(ma(j).special==0 && ma(j).ord!=j) 
		{
			fprintf(stderr, "\nFound: \n");
			for(int k=0; k<=ma.num-1; ++k)
			{
				fprintf(stderr, "%d: %s, bq = %g\n", 
				                k, ma(k).ids, ma(k).betaq);
				// ma(k).plot(dig10(k), dig1(k));
			}
			slamserror("modespectrum: mode set corrupted (1)");
		}
		ma(j).setids(j);
		if(sqrt(fabs(ma(j).betaq)) < k0*nmax*SLAMS_MODEANGLELIM)
		{
			if(quiet == 0) 
		fprintf(stderr, 
		"\n(!)  %s: small |beta^2| = %g, change computational window slightly.\n", 
		ma(j).ids, fabs(ma(j).betaq));
			else fprintf(stderr, "(!)");
		}
	}
	if(quiet != 1) fprintf(stderr, " Ok.\n\n");
	return ma.num;
}

int modespectrum(Waveguide wg, Polarization pol, Interval cw, Boundary_type bt,
                 int mnm, ModeArray& ma)
{
	return modespectrum(wg, pol, cw, bt, mnm, ma, 0);
}

/* compute the discretized mode spectrum on an interval with Dirichlet
   boundaries, returns the number of found modes
   wg: the structure under consideration
   pol: polarization type, TE or TM
   cw: the computational window
   mnm: highest order mode (+1) to be included, determine up to mnm eigenmodes 
   ma: the modes found (output)
   quiet: == 1: suppress log output */
int modespectrum(Waveguide wg, Polarization pol, 
                 Interval cw, int mnm, ModeArray& ma, int quiet)
{
	return modespectrum(wg, pol, cw, LIMD, mnm, ma, quiet);
}

int modespectrum(Waveguide wg, Polarization pol, 
                 Interval cw, int mnm, ModeArray& ma)
{
	return modespectrum(wg, pol, cw, LIMD, mnm, ma, 0);
}