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

/*
 * wmm.cpp
 * WMM mode analysis
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmmat.h"
#include"wmmtrif.h"
#include"mode.h"
#include"wmmfld.h"
#include"wmmarr.h"
#include"wmm.h"
#include"wmmglob.h"
#include"wmmpol.h"
#include"wmmdisc.h"

/* WMM mode analysis: global variables */
Waveguide Wg;         // waveguide structure under investigation
Polarization Pol;     // polarization mode
Symmetry Sym;         // prescribed modal symmetry 
double Neffmin;       // interval for
double Neffmax;       // the propagation constant search
WMM_Parameters Par;   // less specific parameters
double K0;	      // vacuum wavenumber 
double Invomep0;      // 1/omega/epsilon_0  [V*um/A]  
double Invommu0;      // 1/omega/mu_0  [A*um/V]  

char Logname[13] = "___sev__";	
		      // logfile: values of the lsq field mismatch vs. Beta

double Beta;          // propagation constant, current trial value 
double SDBeta;        // trial value for spectral discretization

Imatrix Rtyp;         // rectangle types
Imatrix Styp;         // symmetry dependencies 
Dmatrix Swf;          // normalization weighting factors, symmet. dep.
int Crim;             // y-index of the center rectangles 

		      // edges connecting rectangles:
int NumVk;            // number	
Ivector Vktyp;        // orientation
Dvector Vkbeg;        // coordinates
Dvector Vkend;        //  
Dvector Vkcor;        //  
Imatrix Nrxi;         // indices of neighboring rectangles
Imatrix Nryi;         // indices of neighboring rectangles
Dvector Vkswf;        // symmetry weighting factors 

int AmDim;            // current dimension of the amplitude matrices
double* Amat;         // amplitude matrix  
double* AmVec;        // amplitude vector

Trifunfield Phi;      // mode field prototype


/* default analysis parameters */
WMM_Parameters::WMM_Parameters()
{
	vform = HXHY;
	vnorm = NRMMH;
	ccomp = CCALL;

	ini_d_alpha   = 0.05;
	ini_N_alpha   = 15;
	ini_alpha_max = 2.0;

	ini_steps = 30;
	ref_num   = 5;
	ref_exp   = 4.0;
	ref_sdf   = 0.5;

	fin_d_alpha   = 0.01;
	fin_N_alpha   = 30;
	fin_alpha_max = 3.0;

	btol = 1.0e-7;
	mshift = 1.0e-10;

	penfacTE.strip();
	penfacTM.strip();
	penfacEXT = 1.0;
}

/* free allocated memory in WMM_parameters */
void WMM_Parameters::strip()
{
	penfacTE.strip();
	penfacTM.strip();
	return;
}

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

/* ................................................................... */

/* setup global geometry parameters: helper function */
void addvkpar(int typ, double coord, double beg, double end, 
	      int xi1, int yi1, int xi2, int yi2)
{
	Vktyp(NumVk) = typ;
	Vkcor(NumVk) = coord;
	Vkbeg(NumVk) = beg;
	Vkend(NumVk) = end;
	Nrxi(0, NumVk) = xi1; Nryi(0, NumVk) = yi1;
	Nrxi(1, NumVk) = xi2; Nryi(1, NumVk) = yi2;
	++NumVk;
	return;
}

/* setup global geometry parameters */
void initgeometry()
{
	int l, m, j;

	Rtyp = Imatrix(Wg.nx+2, Wg.ny+2);
	Styp = Imatrix(Wg.nx+2, Wg.ny+2);
	Swf  = Dmatrix(Wg.nx+2, Wg.ny+2);

	NumVk = (Wg.nx+1)*(Wg.ny+2)+(Wg.ny+1)*(Wg.nx+2);
	Vktyp = Ivector(NumVk);
	Vkbeg = Dvector(NumVk);
	Vkend = Dvector(NumVk);
	Vkcor = Dvector(NumVk);
	Nrxi = Imatrix(2, NumVk); 
	Nryi = Imatrix(2, NumVk); 
	Vkswf = Dvector(NumVk);
	NumVk = 0;

	/* rectangle types */
	for(m=0; m<=Wg.ny+1; ++m)
	{
		if(m == 0)          Rtyp(0,m) = 6; 
		if(m>0 && m<=Wg.ny) Rtyp(0,m) = 5;
		if(m == Wg.ny+1)    Rtyp(0,m) = 4;
		for(l=1; l<=Wg.nx; ++l)
		{
			if(m == 0)          Rtyp(l,m) = 7; 
			if(m>0 && m<=Wg.ny) Rtyp(l,m) = 0;
			if(m == Wg.ny+1)    Rtyp(l,m) = 3;
		}
		if(m == 0)          Rtyp(Wg.nx+1,m) = 8; 
		if(m>0 && m<=Wg.ny) Rtyp(Wg.nx+1,m) = 1;
		if(m == Wg.ny+1)    Rtyp(Wg.nx+1,m) = 2;
	}

	/* symmetry weighting factors */
	if(Sym != NOS && ((Wg.ny-1)/2)*2+1 != Wg.ny)
	{
		fprintf(stderr, "### Symmetry restriction not reasonable. ###\n");
		Sym = NOS;
	}
	switch(Sym)
	{
	case NOS:
		for(m=0; m<=Wg.ny+1; ++m) 
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				Styp(l,m) = CALC;
				Swf(l,m) = 1.0;
			}
		}
		break;
	case SYM:
		Crim = (Wg.ny-1)/2+1;
		for(m=0; m<=Crim-1; ++m) 
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				Styp(l,m) = CALC;
				Swf(l,m) = 2.0;
			}
		}
		for(l=0; l<=Wg.nx+1; ++l)
		{
			Styp(l,Crim) = VSYM;
			Swf(l,Crim) = 1.0;
		}
		for(m=Crim+1; m<=Wg.ny+1; ++m) 
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				Styp(l,m) = IDSY;
				Swf(l,m) = 0.0;
			}
		}
		break;
	case ASY:
		Crim = (Wg.ny-1)/2+1;
		for(m=0; m<=Crim-1; ++m) 
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				Styp(l,m) = CALC;
				Swf(l,m) = 2.0;
			}
		}
		for(l=0; l<=Wg.nx+1; ++l)
		{
			Styp(l,Crim) = VASY;
			Swf(l,Crim) = 1.0;
		}
		for(m=Crim+1; m<=Wg.ny+1; ++m) 
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				Styp(l,m) = IDAS;
				Swf(l,m) = 0.0;
			}
		}
		break;
	}

	/* edges */
	NumVk = 0;
	Rect r;
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx; ++l)
		{
			r = Wg.rectbounds(l,m);
			addvkpar(HOR, r.x1, r.y0, r.y1, l, m, l+1, m);
		}
	}
	for(l=0; l<=Wg.nx+1; ++l)
	{
		for(m=0; m<=Wg.ny; ++m)
		{
			r = Wg.rectbounds(l,m);
			addvkpar(VER, r.y1, r.x0, r.x1, l, m, l, m+1);
		}
	}
	/* ... their symmetryweightingfactors */
	if(Sym == NOS)
	{
		for(j=0; j<=NumVk-1; ++j) Vkswf(j) = 1.0;
		return;	
	}
	for(j=0; j<=NumVk-1; ++j)
	{
		
		switch(Vktyp(j))
		{
		case HOR: 
			m = Nryi(0,j);
			if(m < Crim) Vkswf(j) = 2.0;
			if(m == Crim) Vkswf(j) = 1.0;
			if(m > Crim) Vkswf(j) = 0.0;
			break;
		case VER:
			if(Nryi(0,j) <= Crim && Nryi(1,j) <= Crim)
				Vkswf(j) = 2.0;
			else Vkswf(j) = 0.0;
			break;
		}
	}
	return;
}

/* free global matrix memory */
void stripgeometry()
{
	Wg.strip();
	Rtyp.strip();
	Styp.strip();
	Swf.strip();
	Vktyp.strip();
	Vkbeg.strip();
	Vkend.strip();
	Vkcor.strip();
	Nrxi.strip();
	Nryi.strip();
	Vkswf.strip();
}

/* ................................................................... */

/* file for logging lsqerror evaluations */
void openlogfile(char ext0, char ext1)
{
	switch(Pol)
	{
		case QTE: Logname[0] = 't'; Logname[1] = 'e'; break;
		case QTM: Logname[0] = 't'; Logname[1] = 'm'; break;
		case VEC: Logname[0] = 'v'; Logname[1] = 'c'; break;
	}
	switch(Sym)
	{
		case SYM: Logname[2] = 's'; break;
		case ASY: Logname[2] = 'a'; break;
		case NOS: Logname[2] = 'n'; break;
	}
	Logname[6] = ext0;
	Logname[7] = ext1;
	FILE *dat;
	char name[20];
	int i;
	for(i=0; i<=7; ++i) name[i] = Logname[i];
	name[8] = '.'; name[9] = 'x'; name[10] = 'y'; name[11] = 'f';
	name[12] = 0;
	fprintf(stderr, "Log >> %s\n", name); 
	dat = fopen(name, "w+");
	fclose(dat);
	return;
}
void closelogfile() {;}

/* output parameters for the WMM mode analysis */
void writesetupinfo()
{
	int j, k, p;
	fprintf(stderr, 
		"\n------------- WMM -------------------------- '99 ---\n");
	switch(Pol)
	{
	case VEC: 
		fprintf(stderr, "VEC "); 
		switch(Par.vform)
		{
		case HXHY: fprintf(stderr, "HXHY "); break;
		case EXEY: fprintf(stderr, "EXEY "); break;
		case EZHZ: fprintf(stderr, "EZHZ "); break;
		}
		switch(Par.vnorm)
		{
		case NRMMH: fprintf(stderr, "NRMMH "); break;
		case NRMME: fprintf(stderr, "NRMME "); break;
		case NRMMZ: fprintf(stderr, "NRMMZ "); break;
		case NRMMS: fprintf(stderr, "NRMMS "); break;
		}
		switch(Par.ccomp)
		{
		case CCALL: fprintf(stderr, "CCALL "); break;
		case CCHXYZEZ: fprintf(stderr, "CCHXYZEZ "); break;
		case CCEXYZHZ: fprintf(stderr, "CCEXYZHZ "); break;
		}
		break;
	case QTE: fprintf(stderr, "QTE "); break;
	case QTM: fprintf(stderr, "QTM "); break;
	}
	switch(Sym)
	{
	case NOS: fprintf(stderr, "NOS "); break;
	case SYM: fprintf(stderr, "SYM "); break;
	case ASY: fprintf(stderr, "ASY "); break;
	}
	fprintf(stderr, "\n"); 
	fprintf(stderr, 
		"SD-ini    d: %5g,    N: %3d,    al: %5g\n", 
		Par.ini_d_alpha, Par.ini_N_alpha, Par.ini_alpha_max);
	fprintf(stderr, 
		"SD-fin    d: %5g,    N: %3d,    al: %5g\n", 
		Par.fin_d_alpha, Par.fin_N_alpha, Par.fin_alpha_max);
	fprintf(stderr, 
        	"I_steps: %d   Btol: %g\n", Par.ini_steps, Par.btol);
	fprintf(stderr, 
        	"Refine  Num: %d  Exp: %g  Sdf: %g\n", 
		Par.ref_num, Par.ref_exp, Par.ref_sdf);
	fprintf(stderr, 
		"----------------------------------------------------\n");
	fprintf(stderr, "    Nx: %d  Ny: %d\n", Wg.nx, Wg.ny);
	fprintf(stderr, "    Hx: ");
	for(j=0; j<=Wg.nx; ++j) fprintf(stderr, "%6.10g  ", Wg.hx(j));
	fprintf(stderr, "\n");
	fprintf(stderr, "    Hy: ");
	for(j=0; j<=Wg.ny; ++j) fprintf(stderr, "%6.10g  ", Wg.hy(j));
	fprintf(stderr, "\n");
	fprintf(stderr, "     N: ");
	for(j=Wg.nx+1; j>=0; --j) 
	{
		for(k=0; k<=Wg.ny+1; ++k) fprintf(stderr, "%6.10g  ", Wg.n(j,k));
		if(j>0) fprintf(stderr, "\n        ");
	}
	fprintf(stderr, "\n");
	fprintf(stderr, "Lambda: %.10g, K0: %g\n", Wg.lambda, K0);
	fprintf(stderr, " n_eff: ( %.10g, %.10g )\n", Neffmin, Neffmax);
	p = 0;
	for(j=Wg.nx+1; j>=0; --j) 
	{
		for(k=0; k<=Wg.ny+1; ++k) 
		if(Par.penfacTE(j,k) != 1.0 || Par.penfacTM(j,k) != 1.0) 
			p = 1;
	}
	if(p != 0)
	{
		fprintf(stderr, 
		"----------------------------------------------------\n");
		fprintf(stderr, "PenFac TE: ");
		for(j=Wg.nx+1; j>=0; --j) 
		{
			for(k=0; k<=Wg.ny+1; ++k) fprintf(stderr, "%5.10g  ", 
							Par.penfacTE(j,k));
			if(j>0) fprintf(stderr, "\n           ");
		}
		fprintf(stderr, "\n");
		fprintf(stderr, "       TM: ");
		for(j=Wg.nx+1; j>=0; --j) 
		{
			for(k=0; k<=Wg.ny+1; ++k) fprintf(stderr, "%5.10g  ", 
							Par.penfacTM(j,k));
			if(j>0) fprintf(stderr, "\n           ");
		}
		fprintf(stderr, "\n");
	}
	fprintf(stderr, 
		"----------------------------------------------------\n");
	fprintf(stderr, "\n");
	return;	
}	

/* ................................................................... */

/* output setup information */
void diagsetup()
{
	int l, m, vk;
	double beg, end, j;
	FILE *dat;
	dat = fopen("diag", "w+");
	fprintf(stderr, ">> diag\n");
	fprintf(dat,"WMM:\n\n");
	fprintf(dat,"        Pol: %d\n", Pol); 
	fprintf(dat,"        Sym: %d\n\n", Sym); 
	fprintf(dat,"    VecForm: %d\n", Par.vform); 
	fprintf(dat,"VecNormMode: %d\n", Par.vnorm); 
	fprintf(dat,"      Ccomp: %d\n\n", Par.ccomp); 


	fprintf(dat,"* Structure:\n");
	fprintf(dat,"  Nx: %d\n", Wg.nx);
	fprintf(dat,"  Ny: %d\n", Wg.ny);
	fprintf(dat,"  Horizontal boundaries:\n");
	for(l=0; l<=Wg.nx; ++l)
		fprintf(dat,"  %d - %g\n", l, Wg.hx(l));
	fprintf(dat,"  Vertical boundaries:\n");
	for(m=0; m<=Wg.ny; ++m)
		fprintf(dat,"  %d - %g\n", m, Wg.hy(m));
	fprintf(dat,"  refractive index values:\n");
	for(l=0; l<=Wg.nx+1; ++l)
	{
		for(m=0; m<=Wg.ny+1; ++m)
			fprintf(dat, "  [%d,%d] - %g\n", l, m, Wg.n(l,m)); 
	}

	fprintf(dat,"* Wavelength: %g\n", Wg.lambda);
	fprintf(dat,"  K0: %g\n", K0);
	fprintf(dat,"* rectangle types:\n");
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
			fprintf(dat, "  [%d,%d] - %d\n", l, m, Rtyp(l,m)); 
	}
	fprintf(dat,"* edges: # %d\n", NumVk);
fprintf(dat,"  [Nr.] - type  (coord.)  start   end  [Nb1]  [Nb2]  *Swf\n");
	for(vk=0; vk<=NumVk-1; ++vk)
	fprintf(dat, 
		"  [%d] - %d  (%g)  %g->%g   [%d,%d]  [%d,%d]  *%g\n", 
		vk, Vktyp(vk), Vkcor(vk), Vkbeg(vk), Vkend(vk),
		Nrxi(0,vk), Nryi(0,vk), Nrxi(1,vk), Nryi(1,vk), Vkswf(vk)); 
	fclose(dat);
	dat = fopen("vk.ppf", "w+");
	fprintf(stderr, ">> vk.ppf\n");
	for(vk=0; vk<=NumVk-1; ++vk)
	{
		switch(Vktyp(vk))
		{
			case HOR: 
			   beg = Vkbeg(vk); 
			   if(beg < Wg.hy(0)) beg = Wg.hy(0)-Wg.lambda;
			   end = Vkend(vk); 
			   if(end > Wg.hy(Wg.ny)) end = Wg.hy(Wg.ny)+Wg.lambda;
			   for(j=0.0; j<=1.0; j+=0.05)
		           	fprintf(dat, "%g  %g\n", 
			   		beg+j*(end-beg), Vkcor(vk));
			   break;
			case VER: 
			   beg = Vkbeg(vk); 
			   if(beg < Wg.hx(0)) beg = Wg.hx(0)-Wg.lambda;
			   end = Vkend(vk); 
			   if(end > Wg.hx(Wg.nx)) end = Wg.hx(Wg.nx)+Wg.lambda;
			   for(j=0.0; j<=1.0; j+=0.05)
		           	fprintf(dat, "%g  %g\n", 
			   		Vkcor(vk), beg+j*(end-beg));
			   break;
		}
	}
	fclose(dat);
	return;
}

/* ................................................................... */

/* determine the least squares error mju_beta */
double minAIA(double bet, int ev, char mo, int diag, int nolog)
{
	int l, m, d, k, j;
	int td;
	int *dim, *il, *im, *si;
	double* *lmat;
	dim = new int[(Wg.nx+2)*(Wg.ny+2)];
	il  = new int[(Wg.nx+2)*(Wg.ny+2)];
	im  = new int[(Wg.nx+2)*(Wg.ny+2)];
	si  = new int[(Wg.nx+2)*(Wg.ny+2)];
	lmat = new double*[(Wg.nx+2)*(Wg.ny+2)];
	int jl, jm, kl, km;
	int maxdim;
	int numbk;
	int con;
	double *dmat, *ymat;
	double minaia;
	double *xvec;

	if(bet < K0*Neffmin || bet > K0*Neffmax)
		wmmerror("minAIA: wrong beta !");

	Beta = bet;

	fprintf(stderr, " [%c] ", mo);

	fillampmat();

	fprintf(stderr, "(%.10g)", Beta);

	if(diag == 1)
	{
		fprintf(stderr, "\n minAIA:");
		fprintf(stderr, "Beta: %.16g\n", Beta);
		diagsetup();
		diagspecdisc();
		savemat(Amat, AmDim, AmDim, "amat");
		testmatsym(Amat, AmDim, 1.0e-14);
	}

	numbk = 0;
	maxdim = 0;
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
			if(Swf(l,m) != 0.0)
			{
				il[numbk] = l;
				im[numbk] = m;
				si[numbk] = (*Phi(0,l,m,0)).Amsi;
				td = Phi.ntf(0,l,m);
				if(Pol == VEC) td += Phi.ntf(1,l,m);
				lmat[numbk] = new double[td*td];
				matto0(lmat[numbk], td, td);
				d = fillsnmat(l, m, lmat[numbk], 1); 
				if(diag == 1) 
				{
					char nam[5] = "nm00";
					nam[2] = '0'+l;
					nam[3] = '0'+m;
					savemat(lmat[numbk], d, d, nam);
				}

				choldc(lmat[numbk], d, 0);
				if(d > maxdim) maxdim = d;
				dim[numbk] = d;
				++numbk;
			}
		}
	}

	fprintf(stderr, " -");

	dmat = new double[maxdim*maxdim];
	ymat = new double[maxdim*maxdim];

	for(j=0; j<=numbk-1; ++j)
	{
		jl = il[j];
		jm = im[j];
		for(k=0; k<=j; ++k)
		{
			kl = il[k];
			km = im[k];
			if(j == k) con = 1;
			else con = Wg.testconnect(jl, jm, kl, km);
			if(con != 0)
			{
				matto0(dmat, dim[j], dim[k]);
				extractbloc(Amat, AmDim, AmDim, si[j], si[k],
					    dmat, dim[j], dim[k]);
				bwdsubst(lmat[j], dim[j], ymat, dim[k], dmat); 
				transpose(ymat, dim[j], dim[k], dmat);
				bwdsubst(lmat[k], dim[k], ymat, dim[j], dmat); 
				insertbloc(ymat, dim[k], dim[j],
					   Amat, AmDim, AmDim, si[k], si[j]);
			}
		}
	}

	if(ev == 1 || diag == 1) AmVec = new double[AmDim];

	if(diag == 1) 
	{
		for(k=1; k<=AmDim-1; ++k)
		{
			for(j=0; j<=k-1; ++j) Amat[k+j*AmDim] = Amat[j+k*AmDim];
		}

		savemat(Amat, AmDim, AmDim, "rmat");
		xvec = new double[AmDim];
		minaia = smallesteigval(Amat, AmDim, 1, xvec, 0, 1);
		savemat(xvec, AmDim, 1, "xvec");
		fprintf(stderr, "sev = %g\n", minaia);
		for(j=0; j<=numbk-1; ++j)
		{
			vecbwdsubs(lmat[j], dim[j], 
				   &(AmVec[si[j]]), &(xvec[si[j]]));
		}
		delete[] xvec;
		savemat(AmVec, AmDim, 1, "lvec");
		closeampmat();
		delete[] AmVec;
		exit(1);
	}

	fprintf(stderr, "-");

	if(ev != 1) minaia = smallesteigval(Amat, AmDim, 0, NULL, 0, 1);
	else
	{
		xvec = new double[AmDim];
		minaia = smallesteigval(Amat, AmDim, 1, xvec, 0, 1);
		for(j=0; j<=numbk-1; ++j)
		{
			vecbwdsubs(lmat[j], dim[j], 
				   &(AmVec[si[j]]), &(xvec[si[j]]));
		}
		delete[] xvec;
	}
	if(nolog != 1) apptoxyf(Logname, Beta, minaia);
	fprintf(stderr, "> %.10g\n", minaia);

	for(j=0; j<=numbk-1; ++j) delete[] lmat[j];
	delete[] lmat;
	delete[] dmat;
	delete[] ymat;
	delete[] dim;
	delete[] il;
	delete[] im;
	delete[] si;

	closeampmat();
	if(minaia < 0.0) { minaia = minAIA(Beta, 1, 'x', 1, nolog); }

	return minaia;
}

double minAIA(double bet, int ev, char mo, int diag)
{
	return minAIA(bet, ev, mo, diag, 0);
}

/* ................................................................... */

/* setup WMM_Mode for given propagation constant b */ 
WMM_Mode setupmode(double b, int nolog)
{
	double nu, nd;
	int mcp, cp, l, m, j;
	Trifun* fun;
	WMM_Mode mode;
	specdisc(SD_FIN, b, b, b);
	minAIA(b, 1, '*', 0, nolog);
	mcp = 0; if(Pol == VEC) mcp = 1;
	for(cp=0; cp<=mcp; ++cp)
	{
		for(m=0; m<=Wg.ny+1; ++m)
		{
			for(l=0; l<=Wg.nx+1; ++l)
			{
				for(j=0; j<=Phi.ntf(cp,l,m)-1; ++j)
				{
					fun = Phi(cp,l,m,j);
					(*fun).Amplit = AmVec[(*fun).Amsi];
				}
			}
		}
	}
	delete[] AmVec;

	mode.wg = Wg;
	mode.pol = Pol;
	mode.sym = Sym; 
	mode.beta = b;
	mode.k0 = K0;
	mode.neff = b/K0;
	nu = Wg.defaultneffmax();
	nd = Wg.defaultneffmin();
	mode.npcB = (b*b/K0/K0-nd*nd)/(nu*nu-nd*nd);
	mode.vform = Par.vform;
	mode.phi = Phi;
	mode.ampfac = 1.0;
	mode.normalize(1.0);
	mode.setfieldmax();

	return mode;
}
WMM_Mode setupmode(double b)
{
	return setupmode(b, 0);
}

/* ................................................................... */

/* fix propagation constant --- after a minimum is trapped 
 * partly based on 
 * brent.c (C) Copr. 1986-92 Numerical Recipes Software 5.13. */
double trapbeta(double b0, double v0,
	        double b1, double v1,
	        double b2, double v2)
{
	int i;
	double a,b,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
	double e, d;
	double xmin;

//	double b0s, b2s;
	double dx, dxi, dxf;
	double sf, n;		

//	b0s = b0;
//	b2s = b2;
	dxi = (b2-b0);
	dxf = Par.btol*b1;

	w = dxi;
	n = (double) Par.ref_num;

	dx = (b2-b0);
	b1 = 0.5*(b0*b0*(v2-v1)-b1*b1*(v2-v0)+b2*b2*(v1-v0))
                /(b0*(v2-v1)-b1*(v2-v0)+b2*(v1-v0));
	b0 = b1-0.5*dx; 
	b2 = b1+0.5*dx; 
	
	for(i=1; i<=Par.ref_num-1; ++i)
	{
		sf = (n+1.0-Par.ref_sdf*n)*i/(n+1-Par.ref_sdf*i)/n;
		specdisc(sf, b0, b1, b2);
		v0 = minAIA(b0, 0, '-', 0);
		v1 = minAIA(b1, 0, '0', 0);
		if(v1>v0)
		{
			while(v1>v0)
			{
				b2 = b1;
				v2 = v1;
				b1 = b0;
				v1 = v0;
				b0 = b0 - dx;
				if(b0<=K0*(Neffmin+(Neffmax-Neffmin)/100.0))
				{
					return -1.0;
				}
				v0 = minAIA(b0, 0, '-', 0);
			}
		}
		else
		{
			v2 = minAIA(b2, 0, '+', 0);
			while(v1>v2)
			{
				b0 = b1;
				v0 = v1;
				b1 = b2;
				v1 = v2;
				b2 = b2 + dx;
				if(b2>=K0*(Neffmax-(Neffmax-Neffmin)/100.0))
				{
					return -1.0;
				}
				v2 = minAIA(b2, 0, '+', 0);
			}
		}

		w = dxf + (dxi-dxf)*exp(Par.ref_exp*log((n-i)/n));
/*
		while(b2-b0 > w)
		{
			if(b1-b0 > b2-b1)
			{
				b = b1-0.3819660*(b1-b0);
				v = minAIA(b, 0, 's', 0);
				if(v > v1)
				{
					b0 = b;
					v0 = v;
				}
				else
				{
					b2 = b1;
					v2 = v1;
					b1 = b;
					v1 = v;
				}
			}
			else
			{
				b = b1+0.3819660*(b2-b1);
				v = minAIA(b, 0, 's', 0);
				if(v > v1)
				{
					b2 = b;
					v2 = v;
				}
				else
				{
					b0 = b1;
					v0 = v1;
					b1 = b;
					v1 = v;
				}
			}
		} 
		dx = (b2-b0);
		b1 = 0.5*(b0*b0*(v2-v1)-b1*b1*(v2-v0)+b2*b2*(v1-v0))
			/(b0*(v2-v1)-b1*(v2-v0)+b2*(v1-v0));
		b0 = b1-0.5*dx; 
		b2 = b1+0.5*dx; 
*/
		b1 = 0.5*(b0*b0*(v2-v1)-b1*b1*(v2-v0)+b2*b2*(v1-v0))
			/(b0*(v2-v1)-b1*(v2-v0)+b2*(v1-v0));
		b0 = b1-0.5*w; 
		b2 = b1+0.5*w; 
	}

//	if(b0<b0s) b0 = b0s;
//	if(b2>b2s) b2 = b2s;
//	b1 = (b0+b2)/2.0;

	dx = (b2-b0);
	specdisc(SD_FIN, b0, b1, b2);
	v0 = minAIA(b0, 0, '-', 0);
	v1 = minAIA(b1, 0, '0', 0);
	if(v1>v0)
	{
		while(v1>v0)
		{
			b2 = b1;
			v2 = v1;
			b1 = b0;
			v1 = v0;
			b0 = b0 - dx;
			if(b0<=K0*(Neffmin+(Neffmax-Neffmin)/100.0))
			{
				return -1.0;
			}
			v0 = minAIA(b0, 0, '-', 0);
		}
	}
	else
	{
		v2 = minAIA(b2, 0, '+', 0);
		while(v1>v2)
		{
			b0 = b1;
			v0 = v1;
			b1 = b2;
			v1 = v2;
			b2 = b2 + dx;
			if(b0>=K0*(Neffmax-(Neffmax-Neffmin)/100.0))
			{
				return -1.0;
			}
			v2 = minAIA(b2, 0, '+', 0);
		}
	}

	a=b0;
	b=b2;
	x=w=v=b1;
	fw=fv=fx=v1;
	e=0.0; d=0.0;
	for (i=1;i<=100;i++) {
		xm=0.5*(a+b);
		tol2=2.0*(tol1=Par.btol*fabs(x)+1.0e-10);
		if (fabs(x-xm) <= (tol2-0.5*(b-a))) 
		{
			xmin=x;
			return xmin;
		}
		if (fabs(e) > tol1) {
			r=(x-w)*(fx-fv);
			q=(x-v)*(fx-fw);
			p=(x-v)*q-(x-w)*r;
			q=2.0*(q-r);
			if (q > 0.0) p = -p;
			q=fabs(q);
			etemp=e;
			e=d;
			if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
				d=0.3819660*(e=(x >= xm ? a-x : b-x));
			else {
				d=p/q;
				u=x+d;
				if (u-a < tol2 || b-u < tol2)
					d=((xm-x) >= 0.0 
						? fabs(tol1) : -fabs(tol1));
			}
		} else {
			d=0.3819660*(e=(x >= xm ? a-x : b-x));
		}
		u=(fabs(d) >= tol1 ? x+d : 
		x+(d >= 0.0 ? fabs(tol1) : -fabs(tol1)));
		fu=minAIA(u, 0, 'f', 0);
		if (fu <= fx) {
			if (u >= x) a=x; else b=x;
			v=w; w=x; x=u;
			fv=fw; fw=fx; fx=fu;
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				v=w;
				w=u;
				fv=fw;
				fw=fu;
			} else if (fu <= fv || v == x || v == w) {
				v=u;
				fv=fu;
			}
		}
	}
	fprintf(stderr, "### trapbeta: more than 100 iterations ! ###\n");
	xmin=x;
	return xmin;
}

/* ................................................................... */

/* analysis procedures  */
#define MAXNFM 100        // maximum number of modes recognized 
int     wmmprocedure(Waveguide wg,       // the structure under consideration
                     Polarization pol,   // polarization type: QTE, QTM, VEC
                     Symmetry sym,       // enforced symmetry: SYM, ASY, NOS
		     double nmin,        // interval investigated for
		     double nmax,        //   effective mode indices 
		     WMM_Parameters par, // further WMM parameters
                     char ext0,          // extensions for the lsqerror
		     char ext1,          //         log file
                     WMM_ModeArray& mode,   // found modes 
		     int opm)            // 0: complete analysis
		                         // 1: find only fundamental mode 
		                         // 2: write survey to sevlogfile   
{
	double b, s;
	double b0, v0, b1, v1, b2, v2;


	Wg = wg;		
	Pol = pol;
	Sym = sym;
	if(nmin == nmax) { nmax=nmin=0.0; }
	if(nmin > nmax) { b=nmin; nmin=nmax; nmax=b; }
	if(nmin >= 1.0) Neffmin = nmin;
	else Neffmin = Wg.defaultneffmin();
	if(nmax >= 1.0) Neffmax = nmax;
	else Neffmax = Wg.defaultneffmax();
	Par = par;

	if(Par.penfacTE.r == 0)
	{
		Par.penfacTE = Dmatrix(Wg.nx+2, Wg.ny+2);
		Par.penfacTE.init(1.0);
		for(int m=0; m<=Wg.ny+1; ++m)
		{
			Par.penfacTE(      0,m) = Par.penfacEXT; 
			Par.penfacTE(Wg.nx+1,m) = Par.penfacEXT; 
		}
		for(int l=1; l<=Wg.nx; ++l)
		{
			Par.penfacTE(l,       0) = Par.penfacEXT;
			Par.penfacTE(l, Wg.ny+1) = Par.penfacEXT; 
		}
	}
	if(Par.penfacTM.r == 0)
	{
		Par.penfacTM = Dmatrix(Wg.nx+2, Wg.ny+2);
		Par.penfacTM.init(1.0);
		for(int m=0; m<=Wg.ny+1; ++m)
		{
			Par.penfacTM(      0,m) = Par.penfacEXT; 
			Par.penfacTM(Wg.nx+1,m) = Par.penfacEXT; 
		}
		for(int l=1; l<=Wg.nx; ++l)
		{
			Par.penfacTM(l,       0) = Par.penfacEXT;
			Par.penfacTM(l, Wg.ny+1) = Par.penfacEXT; 
		}
	}
	K0 = val_k0(Wg.lambda);
	Invomep0 = val_invomep0(Wg.lambda);
	Invommu0 = val_invommu0(Wg.lambda);

	writesetupinfo();

	initgeometry();
	openlogfile(ext0, ext1);

	Phi = Trifunfield(Wg.nx, Wg.ny);

	mode.clear();	

	s = K0*(Neffmax-Neffmin)/Par.ini_steps;		

/* ------------------------------------------ 
	Beta = K0*(Neffmax+Neffmin)/2.0; 
	diagsetup();
	specdisc(SD_INI, Beta, Beta, Beta);
	diagspecdisc();
	v1 = minAIA(Beta, 0, 'x', 1);
	closelogfile();
	stripgeometry();
	Par.strip();
	exit(1);	
 ------------------------------------------ */

/* ------------------------------------------ 
	Beta = K0*(Neffmax+Neffmin)/2.0; 
//	Beta = 5.887478;
	mode.add(setupmode(Beta));
	closelogfile();
	stripgeometry();
	Par.strip();
	return mode.num;
 ------------------------------------------ */

	if(opm == 2)
	{
		b1 = K0*Neffmin+s;
		while(b1<=K0*Neffmax-s)
		{
			specdisc(SD_FIN, b1, b1, b1);
			v1 = minAIA(b1, 0, 'o', 0);
			b1 += s;
		}
		closelogfile();
		stripgeometry();
		Par.strip();
		return 0;
	}

	b2 = K0*Neffmax-s/2.0;
	b1 = K0*Neffmax-s;
	b0 = K0*Neffmax-2.0*s;
	while(b0>=K0*Neffmin+s/2.0)
	{
		specdisc(SD_INI, b0, b1, b2);
		v0 = minAIA(b0, 0, 'l', 0);
		v1 = minAIA(b1, 0, 'c', 0);
		v2 = minAIA(b2, 0, 'r', 0);
		if(v1<=v0 && v1<=v2)
		{
			b = trapbeta(b0, v0, b1, v1, b2, v2);
			if(b > 0.0)
			{
				if(mode.num < MAXNUMMODES-1)
				{
					mode.add(setupmode(b));
					if(opm == 1)
					{
						closelogfile();
						stripgeometry();
						Par.strip();
						return mode.num;
					}
				}
				else
				{
					fprintf(stderr, "WMM: maximum number %d of modes reached, analysis stopped !\n", MAXNUMMODES);
					closelogfile();
					stripgeometry();
					Par.strip();
					return mode.num;
				}
			}
		}
		b2 = b1; 
		b1 = b0; 
		b0 -= s;
	}

	closelogfile();
	stripgeometry();
	Par.strip();
	return mode.num;
}

/* complete WMM mode analysis, returns number of found modes */
int WMM_modeanalysis(Waveguide wg,       // the structure under consideration
                     Polarization pol,   // polarization type: QTE, QTM, VEC
                     Symmetry sym,       // enforced symmetry: SYM, ASY, NOS
		     double nmin,        // interval investigated for
		     double nmax,        //   effective mode indices 
		     WMM_Parameters par, // further WMM parameters
                     char ext0,          // extensions for the lsqerror
		     char ext1,          //                    log file
                     WMM_ModeArray& mode)   // found modes 
{
	return wmmprocedure(wg, pol, sym, nmin, nmax, par, ext0, ext1, mode, 0);
}

/* WMM mode analysis, find highest (fundamental) mode in [k0*nmin, k0*nmax], 
   returns number of found modes */
int WMM_findfundmode(Waveguide wg,       // the structure under consideration
                     Polarization pol,   // polarization type: QTE, QTM, VEC
                     Symmetry sym,       // enforced symmetry: SYM, ASY, NOS
		     double nmin,        // interval investigated for
		     double nmax,        //   effective mode index 
		     WMM_Parameters par, // further WMM parameters
                     char ext0,          // extensions for the lsqerror
		     char ext1,          //                    log file
                     WMM_ModeArray& mode)   // found mode 
{
	return wmmprocedure(wg, pol, sym, nmin, nmax, par, ext0, ext1, mode, 1);
}

/* WMM mode analysis, survey mju_beta int [k0*nmin, k0*nmax] */
void      WMM_survey(Waveguide wg,       // the structure under consideration
                     Polarization pol,   // polarization type: QTE, QTM, VEC
                     Symmetry sym,       // enforced symmetry: SYM, ASY, NOS
		     double nmin,        // interval investigated for
		     double nmax,        //   effective mode indices 
		     WMM_Parameters par, // further WMM parameters
                     char ext0,          // extensions for the lsqerror
		     char ext1)          //                    log file
{
        WMM_ModeArray mode;  
	wmmprocedure(wg, pol, sym, nmin, nmax, par, ext0, ext1, mode, 2);
	return;
}